\section{Introduction}
\label{sec:intro}

Central to radiative modelling of pulsar wind nebulae (PWN) is the numerically solution the diffusion-loss equation for electron-positron pairs,
\begin{equation}
\label{eq:difflos}
\frac{\partial N(\gamma,t)}{\partial t}=-\frac{\partial}{\partial \gamma}\left[\dot{\gamma}(\gamma)N(\gamma,t) \right]-\frac{N(\gamma,t)}{\tau(\gamma,t)}+Q(\gamma,t),
\end{equation}
where $N(\gamma,t)$ is the electron-positron distribution function with energy $\gamma$ at time $t$.
The full version of the equation is found in \citet{ginzburg1964}, where the radial coordinate and other second order terms are taken into account.
Whereas the lack of radial coordinate $r$ in the equation usually adopted hampers the implementation of a detailed accounting of diffusion of particles and of the magnetic field distribution, it makes the problem more tractable.
Despite this caveat, this kind of simplified models have been very successful in order to fit the spectra of most of the known PWNe and has taught us as a lot about the underlying PWN physics. Among the most important pieces of information gathered we now know (see e.g., \cite{Aharonian1997,zhang2008, gelfand2009, dejager2009, tanaka2011, martin2012, Bucciantini2011, martin2012, torres2013, vorster2013, Torres2014,torres2018,zhu2021}:

\begin{itemize}

\item PWNe are particle dominated: magnetization of a few percent are common, although there might be some observational bias embedded into this fact.

\item PWNe are found having very similar injection parameters (high energy slope: 2.2-2.8; low-energy slope: 1.0-1.6)  and break energies (Lorenz factor around $10^{5 \div 6}$), but there is hardly any correlation among parameters.

\item PWNe have high multiplicity (order $10^{4 \div 6}$).

\item Comparing SEDs of the PWNe as observed today mixes pulsars of different spin-down power and age; generates a variety of distributions. A normalized comparison of the SEDs reduces the dispersion leaving no significant outliers.

\item From the known population of PWNe, only Crab is self-synchrotron Compton (SSC) dominated at high energies (which is understood since we can derive the conditions when this can plausibly happen and prove that essentially only Crab fulfills them). Most others are dominated by the comptonization of the far infrared photon density, with the
SSC being less relevant.

\item  The usually adopted GALPROP code may tend to underpredict local values of FIR and NIR needed for PWNe to shine up in TeV as they do.
GALPROP bins could be coarse for particular PWNe studies.

\item We need to use a more complex dynamics in aged PWNe. Reverberation is very important, but not clearly understood. Reverberation can produce stages of super-efficiency, where, as a result of contraction, e.g., the luminosity in X-rays $L_x > L_{sd}$

\item One key point related to Eq. \ref{eq:difflos} is that, as demonstrated in \citet{martin2012},  all terms of the equation are relevant when studying time evolution, since their relative significance change along time.
Relative differences in flux can reach 100\% between time dependent models once we let them evolve, meaning that the parameters obtained for even a young (worse for a middle-aged PWN) can differ significantly from one model to another depending on the approximations.

\end{itemize}

The recent 3DMHD numerical simulations of the Crab PWN \citep{Porth:2014,Olmi:2016}, despite they are yet limited in time,
showed that the onset of kink instabilities leads to a large amount of dissipation with respect to that observed at lower dimensionality.
Such dissipation renders the total pressure nearly constant across the whole PWN, as well as to a relatively limited
magnetic field gradient.
Given that, by definition, one-zone models have zero gradients, this may explain the success of radiative models to explain the spectral energy distribution of PWNe \citep{Gelfand2017}.

In this context, TIDE is a radiative model prepared to compute the spectral and dynamical evolution of young and middle-aged pulsar wind nebulae (PWNe). It solves the diffusion-loss equation for the electron-positron pairs injected into the nebula from the termination shock taking into account synchrotron, inverse Compton (IC), adiabatic, and Bremmstrahlung energy losses, evolving in time the magnetic field, the radius of the PWN
and other dynamical variables, as
the shock trajectories of the SNR, allowing to follow the evolution further than the beginning of the reverberation phase.
The TIDE code has been used in several occasions to model PWNe observations and study its physics (e.g., see the latest in
\cite{torres2018,bandiera2020} and references therein.
Here, we are presenting the latest significant upgrades, including parallelization of some of its components and automatic fitting in a fully time-dependent scenario, and using these upgrades to tackle the issue of possible degeneracy of solutions.
Real fitting of PWNe multi-frequency data is not usually done (despite it is common to find the word {\it fit} in the literature (see the formerly quoted papers, for instance), it rather refer to visually acceptable descriptions when time-dependent, complex models are involved.
It is thus relevant to briefly clarify what is incorporated in the fitting model and where is room for future improvements.

Using this automatic fitting tool we can more properly address whether degeneracy exists.
Given a relatively complete set of observational data, how unique are PWN models?
Can a significantly different set of model parameters lead to similar spectral energy distributions?
We also use such tool to consider, for the first time to our knowledge, whether a significant set of PWN multi-frequency spectral energy distribution (SED) data can actually imposes a constraint over the age of the system that emits it.
The rest of this work goes directly into presenting the fitting technique and the results obtained,
leaving a more detailed accounting of the physical components to an Appendix.

\section{Fitting method}
\label{sec:min}

\subsection{Identifying free parameters}

Time-dependent PWNe models, despite all simplifying assumptions,
are still subject to a significant degree of freedom.
Table \ref{tab:free-para} identifies all (in principle) free parameters that enter into our model.
The meaning and formulae associated with this set of parameters can be read in Appendix~\ref{sec:appa}.
The parameters in Table \ref{tab:free-para} --unless in obvious cases such as the age of the PWN-- are additionally assumed to be
constant along time, which is in itself an approximation.
There could be occasions (e.g., spin-down rate transitions, glitches, bursts, etc.) where the
injection or the energy distribution or the magnetization could be subject to sudden changes.

\begin{table*}
\scriptsize
\centering
\begin{threeparttable}
\caption{Identification of the whole set of possible free parameters for a time-dependent automatic fitting. All of these parameters can be changed directly in the code, and searched in the minimization procedure when sufficient data warrants it.}
\label{tab:free-para}
\begin{tabular}{lll}
\hline
 & Name & Comment\\
\hline
\hline
$\delta$ & Data systematic uncertainty & Accounts for the systematic uncertainty between data points obtained from different detectors \\
$t_\text{age}$ & Age of the PWN & Relevant for the particle accumulation, dynamics, energetics. \\
$v_\text{psr}$ & Proper motion of the pulsar & Activates/deactivates the injection of new particles depending on the pulsar position with respect to the PWN shell \\
$n$\tnote{$^\dagger$} & Braking index & Affects the evolution of the spin-down luminosity \\
$L_0$\tnote{$^\dagger$} & Initial spin-down luminosity & Affects the evolution of the spin-down luminosity \\
$\tau_0$\tnote{$^\dagger$} & Initial spin-down age & Affects the evolution of the spin-down luminosity \\
$d$ & Distance & Relevant for the flux normalization \\
$\gamma_\text{min}$ & Minimum energy at injection & Relevant for the spectrum normalization and can affect to the luminosity of the spectrum at radio frequencies \\
$\gamma_\text{b}$ & Energy break & Relevant for the spectrum population and the position of the spectrum break \\
$\alpha_l$ & low-energy index & Relevant for the spectral definition at low energies.  \\
$\alpha_h$ & high-energy index  & Relevant for the spectral definition at high energies. \\
$\varepsilon$ & Containment factor & Relevant for fixing the synchrotron cut-off \\
$E_\text{sn}$ & Supernova explosion energy & Relevant for the evolution of the PWN radius \\
$M_\text{ej}$ & Ejected mass & Relevant for the evolution of the PWN radius \\
$\rho_\text{ism}$ & Interstellar medium density & Relevant for the evolution of the PWN radiu \\
$\delta_\text{core}$ & SNR core density index & Relevant for the evolution of the PWN radius \\
$\omega$ & SNR envelope density index & Relevant for the evolution of the PWN radius \\
$\eta_\text{B}$ & magnetic fraction & Controls synchrotron radiation and magnetic field. \\
$\eta_\text{other}$ & Other energy releases fraction & Controls spectrum normalization and magnetic field. \\
$w_\text{i}$ & Photon background energy density & Affects inverse Compton luminosity. Can come from observations or environmental modelling. \\
$T_\text{i}$ & Photon background energy density & Affects inverse Compton luminosity. Can come from observations or environmental modelling. \\
$n_\text{He}/n_\text{H}$ & He/H density ratio & Relevant for Bremsstrahlung radiation spectrum. \\
\hline
\end{tabular}
\begin{tablenotes}
\item [$\dagger$] Their values are taken consistently
with the age of the PWN is the latter is let to vary.
\end{tablenotes}
\end{threeparttable}
\end{table*}

\subsection{Fitting }

Our model now incorporates a new routine package in order to perform fits of multi-wavelength PWNe
spectra from radio to VHE, in a multi-dimensional space, along time. We refer to this as TIDEfit.
It uses the Downhill Simplex Method (also called Nelder-Mead method, see \cite{numerical_recipes}) where the initial set of parameters is represented in a polygon (called simplex) in a $N$-dimensional parameter space, being $N$ the number of free parameters, with $N+1$ vertices.
The Nelder-Mead method prompts a sort of transformations to this polygon, which provide new solutions that are tested, and accepted or rejected according to its nearness to the data, in such a way that the polygon surrounds the final solution up to the desired accuracy for final convergence.
It is known that such a method suffer of a great loss of efficiency when the number of parameters is high enough \citep{gao2012}.
However, these authors proposed a modified version of the method to improve the efficiency in these situations, which is already implemented in the {\it minimize} function of the {\it Scipy} library for Python.

Getting a good fit via methods such as the Monte Carlo Markov Chain implies realizing a high number of iterations  (of the order of tens of thousands or even more) in which the diffusion-loss equation must be solved at every time.
Due to the complexity that implies solving this equation in detail and the potentially high number of free parameters, other typical methods (e.g., minimum least squared) which require the calculations of first partial derivatives (and probably second partial derivatives too) with respect of all the free parameters are not preferred either, due to their time-consumption.
Instead, the simplex method avoid the use of first partial derivatives and can typically converge in $\sim 1000$ iterations for 6-8 free parameters if the initial choice of the simplex is not extremely far from the fitting one.
See Table 1 for typical search ranges.

\subsection{Data treatment}

The multi-wavelength data for a generic PWN come from different detectors, and from different observation times (sometimes decades apart). Thus, in addition of the observational uncertainty of each detector per se and the complexities arising from treating data with vastly different relative error bars, there exist a cross-calibration uncertainty that must be taken into account.
This issue --as well as some misconceptions arising when fitting data in astrophysics-- are discussed by \cite{hogg2010}.
We pose that the most convenient prescription is to introduce an extra free parameter in our fit such that the real uncertainty for each point is systematically enlarged in proportion to the flux calculated from the model.
That is (see \cite{hogg2010}),
\begin{equation}
s_i^2 = \sigma_i^2 + \delta^2 F_i^2,
\end{equation}
where $\sigma_i$ is the observational uncertainty, $\delta$ is the free parameter that accounts for the systematic uncertainty and $F_i$ is the spectral flux obtained from the model.
As discussed in \cite{hogg2010}, the $\chi^2$ function typically used in astrophysics is only a part of the full formula for the likelihood $\chi^2$ function, because the coefficient depending on the uncertainty is neglected. This approximation commonly done in the literature is right if we do not consider systematic errors in our data, so they are constant independently of the values of the model parameters. Here, instead of minimizing the typical $\chi^2$ function, we maximize the likelihood function \citep{hogg2010}
\begin{eqnarray}
\label{eq:likefunc}
\ln p(y | \nu, \sigma, \vec{x}, \delta) = - \frac{1}{2}
\sum_i \left \{ \frac{[y_i(\nu) - F(\nu,\vec{x})]^2}{s_i^2(\nu,\sigma_i,\delta)} - \right.  \nonumber \\  \left.
\ln 2 \pi s_i^2(\nu,\sigma_i,\delta) \right \},
\end{eqnarray}
where $\nu$ is the observed frequency and $\vec{x}$ is our set of free parameters.
Note that Eq.~\ref{eq:likefunc} uses the logarithm of $p$
to have smoother control over the function's gradient and the variation tolerance before consider convergence of the numerical scheme.

If the likelihood function has a complicated shape in the parameter phase-space, it is possible --in the context of the Nelder-Mead method-- to fall into local minima of the likelihood function
depending on the first guess of the parameters.
To avoid this, we noticed that a good practise is to fix first the systematic uncertainty at an a-priori seemingly reasonable value (10 - 20\%),
then fit the rest of parameters,
in order to use the results as an initial guess for a second refined fit now leaving the systematic uncertainty as a free parameter too.
All of this is automatically implemented in our code (see Appendix B for some further notes on implementation and numerical cost).

\section{Applications to two canonical PWNe with varying data quality: How
unique are their fits?}
\label{sec:results}

We now fit the spectra of two canonical PWNe with TIDE: the Crab Nebula and 3C 58. Fortunately, most of the parameters can be fixed for observations in both nebulae, as we can see in Table \ref{tab:fixed}. However, the number of data points and their uncertainty varies.
For instance, whereas the age of Crab is certain, the age of 3C 58 is assumed.
A total of 8 parameters has been left free to fit both PWNe.
Apart the systematic uncertainty, these are the energy break and the spectral slopes of the injection function, the ejected mass of the SNR, the magnetic fraction and the FIR and NIR energy densities.

A special note goes for the containment factor.
The value taken for Crab comes from a previous fit where this parameter is free, but since it only affects the synchrotron cut-off where the number of data points is small, the uncertainty in this parameter is high because the $\chi^2$ changes little when we look for the confidence level.
Thus we have taken it fixed to minimize the number of parameters.
We come back to this in Appendix~\ref{sec:appa}, and show in detail how different containment factors affects the spectra.
The fitting ranges and the results obtained for each nebula are compiled in Table \ref{tab:results}.

\begin{table}
\scriptsize
\centering
\caption{Data from observations -first panel-, computed from the first panel -second panel-, or assumed -third panel- in the models of the Crab nebula and 3C 58.}
\label{tab:fixed}
\begin{tabular}{llcc}
\hline
\hline
Period (ms) & $P$ & 33.4 & 65.7\\
Period derivative (s s$^{-1}$) & $\dot{P}$ & $4.2 \times 10^{-13}$ & $1.93 \times 10^{-13}$\\
Braking index & $n$ & 2.509 & 3\\
Observed radius (pc) & $R_\text{pwn}$ & 1.8 & 2.3\\
Distance (kpc) & $d$ & 2 & 2\\
Pulsed radiation fraction & $\eta_\text{other}$ & 0.0142 & 0.0018\\
\hline
Characteristic age (yr) & $\tau_\text{c}$ & 1296 & 5398\\
Age (yr) & $t_\text{age}$ & 968 & 2500\\
Spin-down  (erg s$^{-1}$) & $L$ & $4.5 \times 10^{38}$ & $2.7 \times 10^{37}$\\
Initial spin-down  (erg s$^{-1}$) & $L_0$ & $3.1 \times 10^{39}$ & $9.3 \times 10^{37}$\\
Initial spin-down age (yr) & $\tau_0$ & 750 & 2878\\
\hline
SN energy explosion (erg) & $E_\text{sn}$ & $10^{51}$ & $10^{51}$\\
ISM density (cm$^{-3})$ & $\rho_\text{ism}$ & 0.5 & 0.1\\
FIR temperature (K) & $T_\text{cmb}$ & 70 & 25\\
NIR temperature (K) & $T_\text{cmb}$ & 5000 & 2900\\
Containment factor & $\varepsilon$ & 0.3 & 0.5\\
\hline
\end{tabular}
\end{table}

\subsection{Crab Nebula (with LHASSO data)}
\label{sec:crab}

\begin{figure*}
\centering
\includegraphics[width=0.45\textwidth]{crab_sed-eps-converted-to.pdf}
\includegraphics[width=0.45\textwidth]{crab_freeage_sed-eps-converted-to.pdf}
\\
\includegraphics[width=0.45\textwidth]{3c58_sed-eps-converted-to.pdf}
\includegraphics[width=0.45\textwidth]{3c58_freeage_sed-eps-converted-to.pdf}
\caption{Models fitted for the Crab Nebula (top) and 3C 58 (bottom). On the left column, the age is fixed and on the right, it is a free parameter (see results in Table \ref{tab:results}) The data points for Crab are taken from \cite{baldwin1971,maciasperez2010}, (compilations of radio observations), \cite{buehler2014} (from IR to $\gamma$-rays and VHE) and the new points offered by the LHAASO collaboration \citep{lhaaso2021}. For 3C 58, they come from \cite{green1986,morsi1987,salter1989} (radio), \cite{torii2000} (infrared), \cite{green1994,slane2008} (X-rays), \cite{abdo2013,ackermann2013,li2018} (gamma-rays) and \cite{aleksic2014} (VHE).}
\label{fig:sed}
\end{figure*}

The Crab Nebula is the most studied PWN, and has precise measurements of its flux at all frequencies from radio to VHE.
Regarding the flux data, we included the compilation shown in \cite{buehler2014} and the recently VHE data points by LHAASO \citep{lhaaso2021}.
The resulting spectral energy distribution (SED) from the automatic fitting process is shown in Figure \ref{fig:sed}.
The automatic fit has reduced $\chi^2$ of 1.013, and the 1$\sigma$ band (depicted in grey in the plot) is almost invisible at all energies in this scale (compare with the companion panel for 3C 58, where this band is more prominently seen).
The fit reproduces the data, particularly in critical regions such as in the synchrotron decay or in the gamma-ray peak.
The Crab PWN is SSC dominated at the gamma-ray peak, with a significant contribution of comptonized CMB photons at larger energies.
Despite the generous range of the parameters search, there are no surprises in the fit, neither in the radiative process dominance nor in the free parameters values.
There is no alternative set of solutions outside the range quoted for the free parameters providing a similarly good fit.
We have tested this thoroughly, starting from dozens of very different initial conditions in the Nelder-Mead simplex:
We performed $\sim 150$ fits with randomly selected initial simplex
obtaining convergence to the solution quoted in more than 73\% of the cases.
In the remaining ones, visual inspection or the value of the reduced $\chi^2$ would immediately qualify the fit as non-compliant.
These non-convergent cases happen when the values of the energy break $\gamma_b$ and/or the magnetic fraction $\eta_\text{B}$, and/or the initial systematic uncertainty
are close to one of the extreme values inside its range of variation.
The Crab nebula spectra is thus uniquely fit in the surrounding of the parameters described, and in the surrounding of no others.
There seems to be little room for degeneracy in our PWNe models when such significant set of data is considered.

It is interesting to remark that the best model actually slightly overpredicts the flux at TeV energies, especially around 100 TeV (see Figure \ref{fig:sedzoom} to zoom in the details).
\begin{figure}
\centering
\includegraphics[width=0.5\textwidth]{crab_sed_zoom-eps-converted-to.pdf}
\caption{Zoom of the Crab Nebula spectrum shown in Figure \ref{fig:sed} between 100 GeV and 1 PeV.}
\label{fig:sedzoom}
\end{figure}

Given that the PWN radius in the model is just an average value representing a more complex physical situation, one can think that the inverse Compton yield at high energies could be reduced by if the radius of the nebula is larger, at least within the corresponding uncertainty.
If so, the magnetic energy density would be lower and, in consequence the magnetic field.
In turn, the magnetic fraction could be larger in order to power up the magnetic field sufficiently, decreasing the normalization of the particle distribution function to compensate.
Despite of the reduction of the normalization, the synchrotron luminosity could be fitted because it is proportional to $\sim B^2$.
A lower number of pairs would decrease the IC flux helping to fit the IC flux in both components (IC-CMB and SSC).
This strategy may however be considered disadvantageous since it would require an ejected mass of about $6M_{\odot}$, which seems too low for current models of supernovae \cite{Smartt:2009}.
\cite{Owen2015} state
that the total mass of gas plus dust in the Crab Nebula is $7.2 \pm 0.5 M_\odot$, consistent with a progenitor star mass of $\sim
9 M_\odot$.
However,
we note on the side that disregarding this fact, and allowing for smaller ejecta masses than the range we adopt in Table
\ref{tab:fixed}, we obtain even better fits than the ones shown.

Another option to decrease the inverse Compton yield at high energies is that part of the energy invested in the acceleration of electron-positron pairs, $(1-\eta_\text{B})$, is used to accelerate an hadronic component, or to emit pulsed emission, or other multi-messenger emission such as neutrinos or gravitational waves.
That is, that part of the energy is simply not available due to it being spent in other processes that are or may be, in fact, naturally occuring.
We have taken this into account. We take the fraction of pulsed emission integrating the data from \cite{Torres2019}, which corresponds to 1.42\% of the spin-down luminosity, and we have actually discounted such a quantity from our contribution to particles (as explained in the Appendix~\ref{sec:appa}). However,  Figure \ref{fig:sed} shows that this imply little change.

Figure \ref{fig:crabchi} shows the shape of the reduced chi-squared function
around the solution given by TIDE for each parameter pair.
Most of the uncertainty in the model is introduced by the subdominant FIR and NIR energy densities: the 1, 2 and 3$\sigma$ ellipsoids are elongated and cover the entire possible value parameter range. This occurs as a result of the fact that none of them is actually relevant for the fit. Interestingly, however, the automatically selected FIR inverse Compton contribution is larger than the NIR one.
The rest of the parameters, those related with the injection function (energy break and low and high energy injection slopes) together with the magnetic fraction and the ejected mass are well determined in the $\chi^2_r$ maps, confirming that the degeneracy is minimal.

\subsubsection{The automatic fit in the context of other models}

The automatically chosen parameters are compatible with earlier studies using previous incarnations of the code (e.g., \cite{Torres2014, aleksic2015}).
Differences introduced by the improved treatment of the dynamics explain 10-20\% changes in the resulting radius and ejected mass (we now discount magnetospheric pulsed emission, use the new SNR/PWN shock positions as in \cite{bandiera2021}, time-dependence of the spin-down luminosity in the radius evolution equations as in \cite{martin2016}, and an improved treatment for the energy conservation in the dynamical evolution as in \cite{bandiera2020}).
The small variation in radius also reflects in the magnetic fraction, while the magnetic field is similar.
Regarding the FIR and NIR energy densities, the fitting algorithm confirms that the IC-FIR and IC-NIR contributions are negligible in order to explain the IC flux.

Apart from TIDE, other models use the new data from LHAASO to update their corresponding fit, but it is difficult to compare them directly due to the features of each model being different.
Here, we will compare our results with three different models. One of them is previous to the LHAASO data, and based on the model of \cite{meyer2010} and published with the MAGIC data points released in \cite{aleksic2015}.
The other two models include the new data: the paper of the LHAASO collaboration itself \citep{lhaaso2021} and the one proposed in \cite{nie2022}.
\cite{meyer2010} (hereinafter MEY) fits the Crab spectrum with a steady one zone model composed by two pair populations: one which accounts for the radio emitting pairs and another one for the wind accelerated particles. These two populations are not calculated by evolving the diffusion-loss equation but  as a fit of the particle spectrum observed today.
The model uses 15 parameters to fit the spectrum.
These are thus significant differences, which makes a direct comparison unfeasible.
The model used in \citep{lhaaso2021} (hereinafter LHA) is also a steady one zone model, but more simplified in comparison with MEY, needing only three parameters, in particular,  the spectral index at high energies, the magnetic field, and the super exponential energy cut-off.
Finally, TIDE and the model of \cite{nie2022} (hereinafter NIE) assume a broken power-law injection of particles and evolve them in time taking into account adiabatic and only synchrotron radiative energy losses.
From the right-hand side of Equation \ref{eq:difflos}, NIE neglects the second term and radiative losses are included in the escape term (third term of Equation \ref{eq:difflos}).
As it was shown in \citealt{martin2012}, this kind of approximations has a significant impact in the results along the time evolution, obtaining differences in the calculated flux of more than a factor 2 by using the same parameters.
Qualitatively bearing in mind account all said differences, selected model values are not far from one another.
For instance, the maximum energy at injection for the models MEY, LHA and NIE is, respectively, 3.7, 2.15 (we take the super-exponential energy cut-off as a reference) and 5.649 PeV (this is the initial maximum energy at injection, so we should take it as an upper limit).
In the case of TIDE, the maximum energy gives 4.7 PeV (see Table \ref{tab:results}) which is an intermediate value of those obtained by these models.
We have seen in \cite{lhaaso2021,nie2022} some trials to include an hadronic component into the spectrum to explain the enhance of the VHE flux at $~$1 PeV in the Crab Nebula.
Introducing this component is certainly possible, but for the current set of data, there is not enough data to really favour a lepto-hadronic model over a leptonic one: it would introduce more free parameters (at least two) to fit only one additional data point.
However, the appearance of such enhancement may require this inclusion in the future.

\subsubsection{A lower-energy explosion}
\label{crab-low}

It has been earlier suggested that SN 1054 was an atypical, low-energy explosion (e.g., \cite{Frail1995,Yang2015}).
Here, we shall also explore this possibility by considering, on the one hand, a model with $E_{sn} = 10^{50}$ erg, and on the other
by additionally consider $E_{sn}$ as a free parameter within a range $10^{49}$ --  $10^{51}$ erg.

For the $E_\text{sn}$ case, the model clearly worsens and requires an ejected mass lower than 7M$_\odot$ to be competitive.
If we prevent this, then we cannot fit the radius of the nebula correctly.
For the latter case, when we allow $E_\text{sn}$ to vary, we obtain a very similar result as show in Table~\ref{tab:results}: $E_\text{sn}=8.3^{+0.5}_{-0.4} \times 10^{50}$ erg, which is close to the a priori canonical value of $10^{51}$ erg as initially assumed.

\subsection{3C 58}
\label{sec:3c58}

3C 58 is a PWN associated with the pulsar PSR J0205+6449 \citep{murray2002}.
Regarding the morphology, this PWN shows a faint jet and a toroidal structure perpendicular to it.
3C 58 has been detected in radio, IR, X-rays, and gamma-rays (see the references in the caption of Figure \ref{fig:sed}).
This makes 3C 58 another good example to test the automatic fitting code, albeit
there is still an ongoing discussion on its age and distance in comparison to the more certain values of the Crab nebula.
Regarding the age, \cite{stephenson2002} associate the nebula to the supernova SN 1181, resulting a current age of 841 yr, but other estimates give an age of up to 7 kyr (see e.g., \citealt{fesen2008} for a compilation).
For our model, we use the dynamical age of 2.5 kyr given by \cite{chevalier2004,chevalier2005} and come back to this discussion below, with a new approach.
We also note that the required ejected mass in order to fit the radius using an age of 841 yr would also be low ($< 8M_\odot$).
On the other hand, the distance ranges from 2 kpc \citep{kothes2010} to 3.2 kpc \citep{roberts1993}.
We take the newest determination of 2 kpc, where new H I data is used.

We obtain a fit with a reduced chi-squared of 0.74 and a systematic uncertainty of 0.01, which is on the limit considered for this parameter.
In terms of the reduced $\chi^2$, the fit is a bit farther from ideal, but we understand that this may be strongly influenced by the number of data points considered in X-rays (189 of 232 in total) and radio to a lesser extent (24 of the remaining 43 taking out the X-ray points), where the fit is really good. The data outside these ranges have errors larger than an order of magnitude, then their contribution to the reduced $\chi^2$ becomes really low.
Regarding the solution convergence (see Figure \ref{fig:3c58chi}), if we compare with the solution for the Crab Nebula, the ellipsoids in the $\chi^2$ maps seem more elongated. Therefore the relative errors of the parameters are higher, but in terms of degeneracy, the ellipsoids are still closed and the solution is not degenerate.

As in the case of the Crab Nebula,  we also discount from the energy input for particles the pulsed radiation emitted ($1.8 \times 10^{-3}$, see Table~\ref{tab:fixed} and \citealt{li2018}), with no noticeable effect on the spectrum and the parameters. The modelled flux surpasses the data points at energies higher than 100 GeV. For Crab, this situation gives a chance to speculate on possible corrections to the model due to the high quality of the data, but in this case, we have less data points and their dispersion makes this difference no-significant in comparison with the variations of the model in other parts of the IC spectrum.

Finally, regarding convergence of the solution, and as in the previous case, we have also performed $\sim$150 simulations with random sets of initial parameters converging to the final solution in more than 70\% of the cases. The initial sets of parameters for the non-convergent cases have in general values close to the limit range for the energy break $\gamma_\text{b}$ and magnetic fraction $\eta_\text{B}$. Thus, a very similar scenario and results as in the case of the Crab Nebula appears.

\subsubsection{The fit in the context of other models}

There are not many models of 3C 58 apart from \cite{tanaka2013} and
the studies we did earlier with TIDE in the framework of different observations
\cite{torres2013,li2018}.
As expected from the outcome of Crab, the differences with the previous fit published in \cite{li2018} are also few, and are mainly explained by the improvement on the expansion model and the mathematical fitting used in the version presented in this paper, as commented before for the case of Crab.

Regarding the model of \cite{tanaka2013}, we recall that the nebula was detected one year later at gamma-rays and VHE \cite{abdo2013,ackermann2013,aleksic2014} and thus these data has not been used then. Despite of this, there is a general good agreement with our obtained slopes for the low and high energy injection power-laws (1.0 \& 3.0, respectively) and the location of the energy break in the two alternative models that they present ($5-9 \times 10^4$).
The expansion model is a ballistic expansion of the nebula where the parameter to be fitted is the expansion velocity.
In this case, the two models are indeed different, from 780 to 2000 km s$^{-1}$.
In our model, the expansion velocity at 2.5 kyr is 1056 km s$^{-1}$, which is roughly in agreement with mentioned models too, but in our case we have to take into account that it is not constant and changes with time like $\sim t^{1/5}$.
Finally, the magnetic field obtained by us is 15.6 $\mu$G, which is close but lower than the minimum value given by \cite{tanaka2013} where the magnetic field lies between 17 and 40 $\mu$G. This small deviation in the magnetic field can be produced by the disagreement of the models on the considered radius of the PWN. While \cite{tanaka2013} consider a radius of 2 pc for 3C 58, our model predicts a slightly larger radius of 2.31, which is closer to what is observed. The magnetic fraction in \cite{tanaka2013} is also a factor 3 lower than what we obtain from the fit.
However it is very difficult to compare both models directly since the evolution in time is computed differently.
In fact, we remark that comparing models based on different theoretical priors provides not more than a qualitative feeling of the stability of the solution with respect to the assumptions.

\begin{table*}
\begin{threeparttable}
\center
\scriptsize
\caption{Parameters obtained for the Crab Nebula and 3C 58 from TIDEfit. The parameters where a searching range is defined are those used to do the fit. The quantities inside the parenthesis are the confidence interval at 95\%.}
\label{tab:results}
\begin{tabular}{llccccc}
\hline
Parameter & Symbol & Range & \multicolumn{2}{c}{Crab Nebula} & \multicolumn{2}{c}{3C 58}\\
 & & & Age fixed & Age free & Age fixed & Age free\\
\hline
\hline
Free parameters \\
\hline
Age (yr) & $t_\text{age}$ & [300,1500 / 2000,5000] & 968 & 914 (901,931) & 2500 & 2108 (2075,2140)\\
Magnetic fraction ($\times 10^{-2}$) & $\eta_\text{B}$ & [0.01,50] & 2.048 (2.027,2.053) & 2.2 (2.1,2.4) & 1.06 (1.03,1.09) & 1.77 (1.70,1.84)\\
Ejected mass ($M_\odot$) & $M_\text{ej}$ & [7,12 / 12,25] & 7.9 (7.6,8.1) & 7.0 (7,7.2) & 17.2 (16.8,17.6) & 10.8 (10.5,11.0)\\
Energy break ($\times 10^5$) & $\gamma_{b}$ & [0.1,100] & 4.9 (4.7,5.3) & 5.8 (5.4,6.2) & 0.88 (0.86,0.90) & 0.76 (0.74,0.78)\\
Low energy index & $\alpha_{l}$ & [1,4] & 1.46 (1.44,1.48) & 1.51 (1.49,1.53) & 1.00 (1,1.03) & 1.00 (1,1.04)\\
High energy index & $\alpha_{h}$ & [1,4] & 2.470 (2.468,2.473) & 2.478 (2.476,2.481) & 3.012 (3.009,3.015) & 3.012 (3.009,3.015)\\
FIR energy density (eV cm$^{-3}$) & $w_\text{fir}$ & [0.1,5] & 0.10 (0.1,0.15) & 0.10 (0.1,0.18) & 0.22 (0.1,0.42) & 0.37 (0.1,0.66)\\
NIR energy density (eV cm$^{-3}$) & $w_\text{nir}$ & [0.1,5] & 0.1 (0.1,1.38) & 0.3 (0.1,2.4) & 0.41 (0.1,0.99) & 0.75 (0.1,1.56)\\
Systematic uncertainty & f & [0.01,0.5] & 0.106 & 0.103 & 0.010 & 0.010\\
\hline
Derived values\\
\hline
Magnetic field ($\mu$G) & $B$ & \ldots &  86.7 (84.4,89.1) & 88.5 (86.6,91.4) & 15.6 (15.3,15.9) & 18.0 (17.6,18.4)\\
Maximum energy at injection ($\times 10^9$) & $\gamma_{max}(t)$ & \ldots & 9.25 (9.21,9.26) & 9.6 (9.4,10.1) & 2.70 (2.66,2.74) & 3.51 (3.44,3.58)\\
PWN Radius (pc) & $R_\text{pwn}$ & \ldots & 1.82 (1.79,1.85) & 1.78 (1.73,1.83) & 2.31 (2.28,2.33) & 2.30 (2.25,2.35)\\
\hline
Statistics\\
\hline
Degrees of freedom\tnote{$^\dagger$} & DOF & \ldots & 580 & 579 & 456 & 455\\
Reduced chi-squared & $\chi^2_\text{r}$ & \ldots & 1.013 & 0.987 & 0.740 & 0.718\\
\hline
\end{tabular}
\begin{tablenotes}
\item [$\dagger$] Age free models have one more free parameter, but also one more data point, the spin-down luminosity. It is necessary to fix correctly the initial spin-down luminosity and the initial spin-down age.
\end{tablenotes}
\end{threeparttable}
\end{table*}

\begin{figure*}
\raggedleft
\includegraphics[width=0.16\textwidth]{crab_gridmap_8_9.jpg}
\includegraphics[width=0.16\textwidth]{crab_gridmap_8_10.jpg}
\includegraphics[width=0.16\textwidth]{crab_gridmap_8_13.jpg}
\includegraphics[width=0.16\textwidth]{crab_gridmap_8_17.jpg}
\includegraphics[width=0.16\textwidth]{crab_gridmap_8_22.jpg}
\includegraphics[width=0.16\textwidth]{crab_gridmap_8_24.jpg}
\\
\includegraphics[width=0.16\textwidth]{crab_gridmap_9_10.jpg}
\includegraphics[width=0.16\textwidth]{crab_gridmap_9_13.jpg}
\includegraphics[width=0.16\textwidth]{crab_gridmap_9_17.jpg}
\includegraphics[width=0.16\textwidth]{crab_gridmap_9_22.jpg}
\includegraphics[width=0.16\textwidth]{crab_gridmap_9_24.jpg}
\\
\includegraphics[width=0.16\textwidth]{crab_gridmap_10_13.jpg}
\includegraphics[width=0.16\textwidth]{crab_gridmap_10_17.jpg}
\includegraphics[width=0.16\textwidth]{crab_gridmap_10_22.jpg}
\includegraphics[width=0.16\textwidth]{crab_gridmap_10_24.jpg}
\\
\includegraphics[width=0.16\textwidth]{crab_gridmap_13_17.jpg}
\includegraphics[width=0.16\textwidth]{crab_gridmap_13_22.jpg}
\includegraphics[width=0.16\textwidth]{crab_gridmap_13_24.jpg}
\\
\includegraphics[width=0.16\textwidth]{crab_gridmap_17_22.jpg}
\includegraphics[width=0.16\textwidth]{crab_gridmap_17_24.jpg}
\\
\includegraphics[width=0.16\textwidth]{crab_gridmap_22_24.jpg}
\caption{Reduced $\chi^2$ maps for Crab Nebula. The dark blue, turquoise and yellow contours correspond to the 1, 2 and 3$\sigma$ confidence levels.}
\label{fig:crabchi}
\end{figure*}

\begin{figure*}
\raggedleft
\includegraphics[width=0.16\textwidth]{3c58_gridmap_8_9.jpg}
\includegraphics[width=0.16\textwidth]{3c58_gridmap_8_10.jpg}
\includegraphics[width=0.16\textwidth]{3c58_gridmap_8_13.jpg}
\includegraphics[width=0.16\textwidth]{3c58_gridmap_8_17.jpg}
\includegraphics[width=0.16\textwidth]{3c58_gridmap_8_21.jpg}
\includegraphics[width=0.16\textwidth]{3c58_gridmap_8_23.jpg}
\\
\includegraphics[width=0.16\textwidth]{3c58_gridmap_9_10.jpg}
\includegraphics[width=0.16\textwidth]{3c58_gridmap_9_13.jpg}
\includegraphics[width=0.16\textwidth]{3c58_gridmap_9_17.jpg}
\includegraphics[width=0.16\textwidth]{3c58_gridmap_9_21.jpg}
\includegraphics[width=0.16\textwidth]{3c58_gridmap_9_23.jpg}
\\
\includegraphics[width=0.16\textwidth]{3c58_gridmap_10_13.jpg}
\includegraphics[width=0.16\textwidth]{3c58_gridmap_10_17.jpg}
\includegraphics[width=0.16\textwidth]{3c58_gridmap_10_21.jpg}
\includegraphics[width=0.16\textwidth]{3c58_gridmap_10_23.jpg}
\\
\includegraphics[width=0.16\textwidth]{3c58_gridmap_13_17.jpg}
\includegraphics[width=0.16\textwidth]{3c58_gridmap_13_21.jpg}
\includegraphics[width=0.16\textwidth]{3c58_gridmap_13_23.jpg}
\\
\includegraphics[width=0.16\textwidth]{3c58_gridmap_17_21.jpg}
\includegraphics[width=0.16\textwidth]{3c58_gridmap_17_23.jpg}
\\
\includegraphics[width=0.16\textwidth]{3c58_gridmap_21_23.jpg}
\caption{Reduced $\chi^2$ maps for 3C 58. The contours have the same correspondence as in Figure \ref{fig:crabchi}.}
\label{fig:3c58chi}
\end{figure*}

\section{Using detailed observations and PWN models to determine the system's age}

The quantity of data available and its completeness in the entire electromagnetic spectrum in both cases allow us to explore using
TIDEfit to estimate the system's age directly from the fit, adding it to the rest of the free parameters.
In Table \ref{tab:results}, we show the results obtained for both nebulae when such approach is adopted.
In terms of reduced $\chi^2$, the quality of the fits are quite similar when comparing with those where the age is fixed, and the age obtained is very close to the expected.
In the case of Crab, the true age falls inside the 2$\sigma$ range obtained from the model.
For 3C 58, for which a real age is not known, we find an age of 2059 yr, which implies a $\sim$440 yr younger PWN than assumed in the previous model.
The change in the age has concurrently changed the values of the magnetic fraction, the ejected mass, and the FIR and NIR energy densities.
The rest of parameters remain vary similar, included derived parameters as the magnetic field, the radius and the maximum energy at injection.
A difference of only 440 years is in fact already pointing towards the larger values in the age dichotomoy stated earlier (where the
discussion was between 871 and 2500 years).
This is the first time that a radiative code is used to fit the age of PWNe.
Of course, it is crucial that a complete set of data, from radio to VHE exists,  to avoid
possible degeneracies when the age is fitted.

\section{Concluding remarks}
\label{sec:conc}

\begin{itemize}

\item We have considered how to take into account that the multi-wavelength data for a generic PWN come from different
detectors, and from different observation times (sometimes separated decades
apart). Using the prescriptions given in \citet{hogg2010}, we proposed that a cross-calibration
uncertainty must be taken into account when fitting models, and posed that introducing an extra parameter to deal with it is the best option.

\item We have successfully incorporated a fitting technique into a fully-time dependent numerical model for PWN.
The appendix~\ref{sec:appa} accounts for the physical components of the model.
They are treated in much the same way as in earlier incarnations of TIDE, prior to automatic fitting,
including some improvements like those related to the dynamics, or the accounting of the energy reservoir spent in non-PWN related emission, and others mentioned.
The difference with former versions of the code, and in general with other codes, is
that all of these components are now included into an automatic fitting technique. This is a result of an improved single-model computational cost (see Appendix~\ref{sec:imple})
as well as the introduction of a relatively fast algorithm technique that avoids the use of derivatives.
The algorithm introduced in the code TIDEfit is based on the Nelder-Mead search for solutions, which was also used recently in
the search of best models for magnetospheric emission from pulsars (see \cite{Iniguez2022}).
TIDEfit is able to cope with variations and minimization of a significant number of parameters (see Table \ref{tab:free-para}).

\item Our model maximizes the likelihood  $\chi^2$ function to search for a solution  instead of minimizing the typical expression for $\chi^2$ function. This is done to take into account systematic uncertainties that appear due to the diverse origin of the data, which is generally ignored in the literature regarding radiative models for PWNe. Our model fits the multi-frequency data once the free parameters of the model
are selected. This is done in a modular way, such that a different number of parameters (and which are those) can be chosen according to the problem at hand
and the number and quality of the data points at our disposal.

\item This represents an improvement over other PWNe models, where the `fits' are in fact qualitatively obtained without a quality assessment.
In addition, it allows us to tackle whether there is any significant degeneracy where very different sets of parameters would lead to equally good fits.
We performed convergence tests with random sets of initial conditions for the cases studied in this paper (Crab Nebula and 3C 58). The model converges to the final solution in more than 70\% of the cases for both nebulae independently of the initial set of parameters.
When not, parameters selected are on the verge of the search ranges and the fits are non-compliant.
There is no significant degeneracy since we find that neither Crab nor 3C58 data can be encompasses by any other set of parameters leading to an equally good solution.

\item Another application we introduced was the determination of the system's age by assuming it as a free parameter of the model.
We find that when a significant dataset is considered, the age can be fixed from the spectrum at a relatively high accuracy, for instance, $<50$ years in the case of Crab.
We have used this possibility to actually provide a new estimation of the age of 3C58, for which a discussion is yet ongoing. We obtain an estimation
of its age at $2108_{-33}^{+32}$ years, within 95\% confidence level.

\end{itemize}

\section*{Data Availability}

No new observational data is herein presented. The dataset used in all figures is publicly available, and references are given correspondingly.
Any additional theoretical detail required is available from the authors on reasonable request.

\section*{Acknowledgements}

We thank
Elena Amato,
Rino Bandiera,
Niccolo Bucciantini,
Agnibha de Sarkar,
Barbara Olmi,
and
Wei Zhang
for comments, discussions, encouragement, and tests.
This work has been supported by the grants INAF grants MAINSTREAM 2018, SKA-CTA, PRIN-INAF 2019, by ASI-INAF n.2017-14-H.O, and the Spanish grant
PID2021-124581OB-I00 and
by the program Unidad de Excelencia María de
Maeztu CEX2020-001058-M.

\bibliographystyle{elsarticle-harv}


\section{Introduction}
\IEEEPARstart{R}{}ECONFIGURABLE Intelligent Surfaces (RISs) are attracting much consideration as a leading technology to achieve intelligent wireless channels environment for the next generation networks \cite{basar2019wireless}. RISs are planar surfaces of electromagnetic (EM) material comprising a large number of cheap passive reflecting elements. A microcontroller controls each element to alter the amplitude and phase of the reflected signal. The RIS technology has many advantages, including the ability to change transmission environments into intelligent ones, enhancing the quality of the received signals at the destination,  reducing the power consumption compared with other technologies, increasing the physical layer security (PLS), and alleviating the undesired interference  \cite{wu2019towards, wu2021intelligent, RISKhoshafa, 9530403}. {Passive RISs prototypes were assembled in \cite{fara2022prototype,  trichopoulos2022design,  dai2020reconfigurable} to acquire more practical and precise results regarding the actual performance of RISs-aided systems by taking experimental measurements.}

With the envision new technologies and utilizing higher frequency bands, secure communications are significant in the sixth generation (6G) wireless networks, where new security challenges arise \cite{chorti2022context}. {Present research contributions have established RISs as cutting-edge technology, with promising research directions toward the 6G. To take things further, integrating RISs with emerging communication technologies results in higher performance gains that can be achieved \cite{basharat2021reconfigurable}}.  The PLS, initially investigated by Wyner \cite{wyner1975wire},  has evolved as an attractive technique for improving the cellular network's secrecy performance against signal leakage. In this respect,  PLS utilizes the natural properties and characteristics of wireless communication channels and noise to secure data transmission by limiting the amount of data that can be leaked at the bit level by eavesdroppers. Thanks to their distinctive characteristics, which enable them to control the transmission environment, RISs can be utilized to eliminate interference and improve the received signal without using active elements. In this respect, the RIS technology has been recently utilized to improve the PLS of wireless communication system \cite{9652031,guan2020intelligent,RISKhoshafa}. To guarantee a secure transmission, the RIS was deployed near the eavesdropper to cancel out the eavesdropping signal received by the eavesdropper \cite{zhang2021improving}, which can actually decrease the information leakage to enhance the PLS of the wireless network.

On the other hand, mobile  wireless communication has experienced rapid development in data traffic due to the dramatic growth of smart devices. According to Cisco, the average number of mobiles per capita will be 3.6 by 2023 \cite{cisco}, leading to an enormous demand for radio spectrum resources, including bandwidth and energy. Consequently, spectral and energy efficiency are two crucial principles for designing future wireless networks. Cognitive radio (CR) has been introduced as an efficient technique to improve spectral efficiency. In CR networks, the spectrum can be shared by two different networks, the primary network (PN) and the secondary network (SN), provided that the interference produced by the SN to the PN is controlled by interference constraint. {The authors in \cite{khoshafa2021performance} studied the PLS of energy harvesting for CR networks using the cooperative relaying technique.
} A common technique is to employ beamforming to improve the performance of the SN while guaranteeing that the interference power received by the PN users is below the predefined interference limit. Nevertheless, the beamforming gain is restricted when the link between the SN transmitter and SN receiver is weak due to severe attenuation. To address such a problem, an RIS can be deployed to improve the performance of the SN while enhancing the secrecy rate of the PN \cite{pan2021reconfigurable}. In \cite{zhang2020intelligent}, the RIS technology has been employed to aid data transmission in CR networks. The authors in \cite{yuan2020intelligent} proposed an RIS-assisted CR network to enhance the SN's achievable rate. In this work, we propose a secure RISs-aided underlay CR network. To the best of our knowledge, there is no previous work studying the advantage of the RIS technology to secure CR network. Furthermore, the influence of the RIS technology on the PN secrecy capacity is investigated. The main contributions of this paper can be summarized as follows:
\begin{figure}[t]
\centering
\includegraphics[width=3.2in]{Figure01.eps}
\caption{{System Model.}}
\label{Sys}
\end{figure}
\begin{itemize}
\item The RIS technology is introduced to improve the reliability of the SN, while concurrently increasing the physical layer security of the PN.

\item To compensate for the spectrum sharing, the RIS technology is utilized as a friendly jammer to ensure a high-secrecy performance for the PN, consequently enabling a win-win situation between the two networks, i.e., security provisioning for the PN and high reliability and robustness for the SN.
\item The SN outage probability is studied, and a novel analytical expression is derived. Besides, closed-form expressions for the secrecy outage probability (SOP) of the PN are also derived, considering two combining techniques, namely maximal ratio combining (MRC) and selection combining (SC), which are employed by an eavesdropper.
\item Asymptotic analysis is provided for the SOP of the PN. Moreover, the benefits of the proposed system model are confirmed through numerical and simulation results.
\end{itemize}

\begin{table}
\centering
\caption{{Table of Symbols}}
\label{table}
\setlength{\tabcolsep}{3pt}
\begin{tabular}{|p{25pt}|p{180pt}|}
\hline
{Symbol}&
{Description}\vspace{2mm} \\
\hline
S&
Secondary transmitter \\
D&
Secondary receiver\\
PT&
Primary transmitter\\
PR&
Primary receiver\\
Eav&
Eavesdropper\\
$N_P$&
Number of antennas at PT\\
$N_E$&
Number of antennas at Eav\\
$N$&
Number of reflecting elements at RIS\\
$h_{ab}$&
Channel coefficient of $ab$ link\\
$d_{ab}$&
Euclidean distance between $ab$\\
$n_{a}$&
AWGN node at $a$\\
$\sigma^2_a$&
AWGN variance at node $a$\\
$y_a$&
Received signal at node $a$\\
$x_s$&
SN transmitted signal \\
$x_p$&
PN transmitted signal \\
$P_\textup{S}$&
SN transmitted power\\
$P_\textup{P}$&
PN transmitted power\\
$d_o$&
Reference distance\\
$\phi_{i}$&
Phase coefficient of the $i^{th}$ element of the RIS\\
$\theta _{i}$&
Residual phase errors affecting the PR\\
$\psi_{i}$&
Residual phase errors affecting the Eav\\
$\gamma_a$&
SINR at node $a$\\
$\bar{\gamma}_a$ &
Average SNR at node $a$\\
$\eta$&
Path loss exponent\\
$f_{X}(\cdot)$&
PDF of random variable $X$\\
$F_{X}(\cdot)$&
CDF of random variable $X$\\
$F^{ \infty }_{X}(\cdot)$&
Asymptotic CDF of random variable $X$\\
$\mathcal{Q}$&
Threshold of interference temperature\\
$\mathcal{C}_\textup{S}$&
PN secrecy capacity\\
$\mathcal{C}_\textup{P}$&
PN capacity\\
$\mathcal{C}_\textup{E}$&
Eav capacity\\
$\mathcal{R}_d$&
SN achievable data rate\\
$\mathcal{R}_s$&
PN target secrecy rate\\
$\textrm{Pr}(\cdot)$&
Probability of an event\\
SOP&
Secrecy outage probability\\
$P_{out}$&
SN outage probability\\
$\textup{SOP}^{ \infty }$&
Asymptotic SOP\\
$\mathcal{G}_a$&
Secrecy diversity order\\
$\mathcal{G}_d$&
Secrecy array gain\\
$\mathcal{O}(\cdot)$&
Higher order term\\
$\beta$&
$2^{\mathcal{R}_s}$\\
$\alpha$&
$\beta-1$ \\
\hline
\end{tabular}
\label{tab1}
\end{table}

\section{System Model}
\label{S01}
A proposed RISs-aided underlay CR system, including a license-holding PN and an unlicensed SN, is considered as shown in Fig. \ref{Sys}. Specifically, the SN comprises a secondary transmitter (S) and a secondary receiver (D), each equipped with a single antenna, while the PN comprises a primary transmitter (PT) and a single-antenna primary receiver (PR). The PT is equipped with $N_P\geq1$ antennas. In addition, an eavesdropper (Eav), equipped with $N_E\geq1$ antennas, intends to overhear the PN's data streams. Therefore, an RIS, made of $N$ reflecting elements, is utilized to enhance the achievable secrecy rate of the PN by interfering with the eavesdropping signals at Eav while improving the transmission conditions of the SN.  {It is worth mentioning that a field-programmable gate array (FPGA) can be utilized as a controller to achieve adjustable control of the RIS in practice, which often communicates and coordinates with other network elements (e.g., BS and users) via dedicated connections  \cite{shi2022wireless}.} It is assumed that the channel state information (CSI) of all channels employed in the system is known\footnote{The traditional pilot signaling techniques can be utilized to estimate the CSI of the legitimate transmission links \cite{basar2019wireless,  zhang2020intelligent, guan2020intelligent, 76543, 9291402} }. As we consider a passive Eav, its CSI is unknown at both the RIS and the PT.

The channel coefficients for the
$\textup{PT} \to \textup{PR} ,\:
\textup{PT} \to \textup{Eav} ,\:
\textup{RIS} \to
\textup{Eav},\:
\textup{RIS} \to
\textup{PR},\:
\textup{S} \to \textup{Eav},\: \textup{S} \to \textup{PR},\: \textup{S} \to \textup{D} ,\:\textup{S} \to \textup{RIS},$ and $\textup{RIS} \to \textup{D}$ links are expressed as $h_{pp},\: h_{pe},\: \textbf{\textrm{h}}_\textbf{\textrm{E}}\footnote{ The bold font is used to indicate vectors.},\: \textbf{\textrm{h}}_\textbf{\textrm{P}},\:h_{se}, \:h_{sp},\: h_{sd},\: \textbf{\textrm{h}}_\textbf{\textrm{S}},$ and $\textbf{\textrm{h}}_\textbf{\textrm{D}},$ respectively. The above channel coefficients are assumed to undergo Rayleigh fading\footnote{Since the transmission links experience blockages and the RIS's location cannot be optimized to guarantee reliable line-of-sight links, similar to \cite{ guan2020intelligent, 76543, 9291402}, Rayleigh fading environment is assumed in this work.
}. {To elaborate, $\textbf{\textrm{h}}_\textbf{\textrm{S}}\in \mathbb{C}^{N\times 1},\: \textbf{\textrm{h}}_\textbf{\textrm{D}}\in \mathbb{C}^{1\times N},\: \textbf{\textrm{h}}_\textbf{\textrm{E}}\in \mathbb{C}^{1\times N},$ and $
\textbf{\textrm{h}}_\textbf{\textrm{P}}\in \mathbb{C}^{1\times N}$ denote the channel vector between the SN transmitter and RIS, RIS and SN receiver, RIS and eavesdropper, and RIS and PN receiver, respectively.}  Moreover, the Euclidean
distances between $\textup{S} \to \textup{RIS} ,\: \textup{RIS} \to \textup{D},\:
\textup{S} \to \textup{D} ,\:
\textup{PT} \to \textup{PR} ,\:
\textup{PT} \to \textup{Eav} ,\:
\textup{RIS} \to
\textup{Eav},\:
\textup{RIS} \to
\textup{PR},\:
\textup{S} \to \textup{Eav},$ and $ \textup{S} \to \textup{PR}$ links are denoted as $d_{sr},\:d_{rd},\: d_{sd},\: d_{pp},\: d_{pe},\: d_{re},\: d_{rp},\:d_{se},$ and $\:d_{sp},$ respectively. In addition, $n_{\varrho  }$ is the additive white
Gaussian noise (AWGN) at D, PR, and E, respectively, where $\varrho  \in \left \{ d, p, e \right \}$, with zero mean and variance  $\sigma^{2}_{\varrho  }$. Consequently, the received signal at D can be written as
\begin{equation}\begin{split}    y_{\textup{D}} =\sqrt{P_\textup{S}}\,x_s&\left [ \left ( \frac{d_{sr}\,  d_{rd}}{d_{o}^{2}} \right )^{-\frac{\eta }{2}}\sum_{i=1}^{N} h_{s_i}\, h_{d_i}\, e^{j\phi_{i}}\right.\\&\left.+h_{sd}\left ( \frac{d_{sd}}{d_{o}} \right )^{-\frac{\eta }{2}}  \right ]+n_d,  \end{split}\end{equation}
where $x_s$ is the SN transmitted signal, $P_\textup{S}$ denotes the SN transmitted power, $d_{o}$ is a reference distance, and $\eta$ is the path loss exponent. In addition, $h_{s_i}$ and $h_{d_i}$ are complex Gaussian random variables (RV) with a zero mean and unit variance, and $\phi_{i}$ is the alterable phase coefficient of the $i^{th}$ element of the RIS. Moreover, it is assumed that the PT position is distant from the RIS and D and, therefore, does not impose any real interference. Consequently, the interference at the RIS and D from PT is negligible; this is a well-known assumption that is widely used in the literature  \cite{lee2010outage, lei2017secrecy}. Furthermore,  the phases of the channels $h_{s_i}$ and $h_{d_i}$ are assumed to be perfectly known at the RIS. Hence, the optimal phase shift is selected to maximize the instantaneous signal-to-noise ratio (SNR) at D \cite{basar2019wireless, bjornson2019intelligent}.
Besides, the reflected gain of the $i^{th}$ reflecting element is assumed to equal to one \cite{guan2020intelligent, basar2019wireless}. Thus, the received signal at the $\textup{PR}$ can be written as
\begin{equation}\begin{split}   y_{\textup{P}} =&\sqrt{P_\textup{P}}\left ( \frac{d_{pp}}{d_{o}} \right )^{-\frac{\eta }{2}} \,h_{pp}\,x_p+\sqrt{P_\textup{S}}\,\, x_{s}\left [ \left ( \frac{d_{sr}\,  d_{rp}}{d_{o}^{2}} \right )^{-\frac{\eta }{2}} \right.\\&\left.\times\sum_{i=1}^{N} h_{s_i}\, h_{p_i}\, e^{j\theta_{i}}+h_{sp}\left ( \frac{d_{sp}}{d_{o}} \right )^{-\frac{\eta }{2}} \right ]+n_p,   \label{yp}    \end{split}\end{equation}
where $x_p$ is the PN transmitted signal, $P_\textup{P}$ denotes the PN transmitted power, and $\theta _{i}$ is the residual phase errors affecting the PR. In a similar way, the wiretapped signal at Eav can be written as
\begin{equation}\begin{split}    y_{\textup{E}} =&\sqrt{P_\textup{P}}\left ( \frac{d_{pe}}{d_{o}} \right )^{-\frac{\eta }{2}}\,h_{pe}\,x_p+\sqrt{P_\textup{S}}\,x_{s} \left [\left ( \frac{d_{sr}\,  d_{re}}{d_{o}^{2}} \right )^{-\frac{\eta }{2}}\right.\\&\left.\times\sum_{i=1}^{N} h_{s_i}\, h_{e_i}\, e^{j\psi_{i}} +h_{se}\left ( \frac{d_{se}}{d_{o}} \right )^{-\frac{\eta }{2}}\right ]+n_e,\end{split}\end{equation}
where $h_{e_i}$ is the channel coefficient between Eav and the $i^{th}$ reflecting element of the RIS and $\psi_{i}$ is the residual phase errors affecting the Eav. The instantaneous SNR at D, $\gamma_{D}$, is given by
\begin{equation}
    \gamma_{\textup{D}}
=\frac{P_\textup{S}}{\sigma_d^2}\left |\left ( \frac{d_{sr}  d_{re}}{d_{o}^{2}} \right )^{-\frac{\eta}{2}}\sum_{i=1}^{N} h_{s_i} h_{d_i} e^{j\phi_{i}}+  h_{sd}\left ( \frac{d_{sd}}{d_{o}} \right )^{-\frac{\eta}{2}} \right |^{2}.
    \label{S11}
\end{equation}

Moreover, the instantaneous signal-to-interference-and-noise ratio (SINR) at PR and Eav, denoted as $\gamma_P$ and $\gamma_E$, respectively, are given by

\begin{equation}
    \begin{split}
        \hspace {-0.5pc}\gamma_{\textup{P}} \hspace {-0.2pc}=\hspace {-0.2pc}\frac{P_\textup{P}\left ( \frac{d_{pp}}{d_{o}^{2}} \right )^{-\eta}\left | h_{pp}\right |^{2}}{P_\textup{S}\hspace {-0.2pc}\left |\left ( \frac{d_{sr}  d_{rp}}{d_{o}^{2}} \right )^{-\frac{\eta}{2}}\sum_{i=1}^{N} h_{s_i} h_{p_i}e^{j\theta_i }\hspace {-0.2pc}+\hspace {-0.2pc}  \left ( \frac{d_{sp}}{d_{o}} \right )^{-\frac{\eta}{2}}h_{sp} \right |^{2}\hspace {-0.3pc}+\sigma_p^2}\hspace {-0.2pc},
        \label{YPP25}
            \end{split}
\end{equation}

and

\begin{equation}
    \begin{split}
        \hspace {-0.2pc}\gamma_{\textup{E}} &\hspace {-0.2pc}=\hspace {-0.3pc}\frac{P_\textup{P}\left ( \frac{d_{sd}}{d_{o}} \right )^{-\eta}\left | h_{pe}\right |^{2}}{P_\textup{S}\left |\left ( \frac{d_{sr}  d_{re}}{d_{o}^{2}} \right )^{-\frac{\eta}{2}}\hspace {-0.2pc}\sum_{i=1}^{N} h_{s_i} h_{e_i}e^{j\psi_i }+\left ( \frac{d_{se}}{d_{o}} \right )^{-\frac{\eta}{2}}h_{se}  \right |^{2}\hspace{-0.3pc}+\hspace {-0.3pc}\sigma_e^{2}},
        \label{YYE}
    \end{split}
\end{equation}
which can be rewritten as
\begin{align}   \begin{split}       \gamma_{\textup{E}} =\frac{\Psi_\textup{PE}}{\Psi_\textup{E}+1},\label{c756}    \end{split}\end{align}
where
\begin{align*}   \begin{split}       \Psi_\textup{E}=\left |\Lambda_{1} \sum_{i=1}^{N} h_{s_i}\, h_{e_i}\,e^{j\psi_i }+\Lambda_{2}\,h_{se}  \right |^{2},\label{c7569}    \end{split}\end{align*}
$\Lambda_{1}=\sqrt{\bar{\gamma }_{se}}\left ( \frac{d_{sr}\,  d_{re}}{d_{o}^{2}} \right )^{-\frac{\eta}{2}}
,$ $\Lambda_{2}=\sqrt{\bar{\gamma }_{se}}\left ( \frac{d_{se}}{d_{o}^{2}} \right )^{-\frac{\eta}{2}}
,$ $\bar{\gamma }_{se}=\frac{P_\textup{S}}{\sigma_e^{2}},$ $\Psi_\textup{PE} =\omega_{e}\left | h_{pe}\right |^{2},$ $\omega_{e} =\bar{\gamma }_{e}\,\left ( \frac{d_{sd}}{d_{o}} \right )^{-\eta},$ and $\bar{\gamma }_{e}=\frac{P_\textup{P}}{\sigma_e^{2}}.$ As the phase shifts of the RIS elements are designed based on the legitimate SN link, the resulting phase distributions for each  RIS$\rightarrow $E link $\psi_{i}$ are i.i.d. and uniformly distributed RVs by virtue of \cite{76543}. Thus, $\Psi_\textup{E}$ can be approximated by an exponential RV according to  \cite[Corollary 2]{badiu2019communication} with a parameter $\Lambda_{{E}}=N\, \Lambda_{1}^{2}+\Lambda_{2}^{2}$. Therefore, the  PDF of $\Psi_\textup{E}$ is given by
\begin{equation}
    f_{\Psi_\textup{E}}\left ( \gamma  \right )=\frac{1}{\Lambda_{E}}\,\exp{\left (-\frac{\gamma }{\Lambda_{E} }  \right )}.\label{c76}
\end{equation}

It is worth mentioning that the SN transmitted power, $P_\textup{S}$, must be under a certain level to limit the interference. Consequently, the interference signal power towards the PR should be constrained as $P_\textup{S}\,\Psi_\textup{P}\leq \mathcal{Q}$, where $\Psi_\textup{P}$ is the summation of the channel power gains from the RIS and S, and $\mathcal{Q}$ is the threshold of interference temperature. Precisely, $\Psi_\textup{P}$ consists of the reflected link $\textup{RIS} \to\textup{PR}$ and the $\textup{S} \to\textup{PR}$ link, while $\mathcal{Q}$ expresses the maximum tolerant interference imposed on the PR. From (\ref{YPP25}), $\Psi_\textup{P}$ can be expressed as
\begin{equation}\begin{split}    \Psi_\textup{P}&=\left |\left ( \frac{d_{sr}  d_{rp}}{d_{o}^{2}} \right )^{-\frac{\eta}{2}}\sum_{i=1}^{N} h_{s_i} h_{p_i}e^{j\theta_i }+  \left ( \frac{d_{sp}}{d_{o}} \right )^{-\frac{\eta}{2}}h_{sp} \right |^{2}.\end{split}    \end{equation}
Similar to $\Psi_\textup{E}$, $\Psi_\textup{P}$ also can be approximated as an exponential RV with a parameter $\lambda_{\textup{P}}=\left ( \frac{d_{sr}  d_{rp}}{d_{o}^{2}} \right )^{-\eta}N+\left ( \frac{d_{sp}}{d_{o}} \right )^{-\eta}$, where the PDF of $\Psi_\textup{P}$ is given by
\begin{equation}\begin{split}    f_{\Psi_\textup{P}}\left ( \gamma  \right )=\frac{1}{\lambda_{P}}\,\exp{\left (-\frac{\gamma }{\lambda_{P} }  \right )}. \label{ccp}\end{split}    \end{equation}

\section{Performance Analysis}
\subsection{PN Secrecy Outage Probability}
In this subsection, the  $\textup{SOP}$ of the PN is investigated. The SOP can be expressed as
\begin{align}   \begin{split}       \textup{SOP}&=\textrm{Pr}\left ( \mathcal{C}_{\textup{S}} < \mathcal{R}_s \right ),\label{c6}    \end{split}\end{align}
where $\mathcal{C}_{\textup{S}}$ is the PN secrecy capacity and $\mathcal{R}_s$ is the PN target secrecy rate. In this regard, $\mathcal{C}_{\textup{S}}$ can be obtained by
\begin{align}   \begin{split}       \mathcal{C}_{\textup{S}}= \big [\mathcal{C}_{\textup{P}}-\mathcal{C}_{\textup{E}} ,0 \big ]^{+},\label{c612}    \end{split}\end{align}
where $\mathcal{C}_{\textup{P}}$ and $\mathcal{C}_{\textup{E}}$ are the PN and the Eav capacities, respectively, and $\left [x,0 \right ]^{+}=\max\left ( x,0 \right )$. Accordingly, $\mathcal{C}_{\textup{P}}$  is given by
\begin{align}   \begin{split}       \mathcal{C}_{\textup{P}}&=\textrm{log}_2\Big ( 1+\gamma_{\textup{P}}   \Big ),\label{c61}    \end{split}\end{align}
where $\gamma_{\textup{P}} $ is given by

\begin{equation}
    \begin{split}
        \gamma_{\textup{P}} =\frac{P_\textup{P}\left ( \frac{d_{pp}}{d_{o}^{2}} \right )^{-\eta}\left | h_{pp}\right |^{2}}{P_\textup{S}\Psi_\textup{P}+\sigma_p^2}=\frac{P_\textup{P}\left ( \frac{d_{pp}}{d_{o}^{2}} \right )^{-\eta}\left | h_{pp}\right |^{2}}{\mathcal{Q}+\sigma_p^2}=\Phi\left | h_{pp}\right |^{2},
        \label{YPP2}
            \end{split}
\end{equation}
where $\mathcal{Q}=P_\textup{S}\,\Psi_\textup{P},$ $\Phi=\omega_{p}\,\vartheta,$ $\omega_{p}=\frac{P_\textup{P}}{\sigma_p^2}\left ( \frac{d_{pp}}{d_{o}^{2}} \right )^{-\eta},$ and $\vartheta  =\left (\frac{\mathcal{Q}}{\sigma_p^2}+1   \right )^{-1}$. {Antenna
selection approach is employed at the PT to avoid the high hardware complexity while maintaining the diversity and reliability advantages of multiple antennas. More specifically, the importance
of using the antenna selection approach lies in the fact that the power consumption and the
complexity of signal processing overhead are low as compared with other techniques such as
beamforming techniques.} Antenna selection strategy is applied at the PT to maintain multiple antennas' diversity and reliability benefits, while avoiding high hardware complexity. Therefore, the best antenna at PT is selected according to the following criterion
\begin{equation}
    \left | h_{pp} \right |^{2}=\:\underset{n \in \left\{1, ..., N_P \right\} }{\textrm{max}}\left | h_{p_{n}p} \right |^{2}.
\end{equation}
The CDF of ${\gamma_{\textup{P}} }$ is given by

\begin{equation}
     F_{\gamma_{\textup{P}} }(\gamma)=\sum_{n=0}^{N_P-1}\frac{N_P\binom{N_P-1}{n}}{(-1)^{-n}\left ( n+1 \right )}\left (1-\textup{exp}\left (  \frac{-\gamma\:(n+1)}{\Phi} \right ) \right ).
     \label{s1021}
 \end{equation}
Moreover, $\mathcal{C}_{\textup{E}}$ is given by $\mathcal{C}_{\textup{E}}=\textrm{log}_2\left ( 1+\gamma_{\textup{E}}   \right )$, where $\gamma_{\textup{E}}$ is given in (\ref{YYE}). Now, the SOP can be derived as

\begin{equation}
    \textrm{SOP}_{{\varsigma  }_{}}=\int_{0}^{\infty }F_{{\gamma}_P}(\beta \gamma+\alpha )f_{{\gamma}_{E_{}} }^{\varsigma  }(\gamma )\:d\gamma,
    \label{M018}
\end{equation}
where $\alpha=\beta-1$, $\beta=2^{\mathcal{R}_s}$, and $\varsigma \in\left \{ \textrm{SC},\textrm{MRC} \right \}$. For the SC technique, the PDF of $\gamma_{E_{}}$, $f_{\gamma_{E_{}}}^{\textrm{SC}}(\gamma)$, is given by \cite[eq. (31)]{khoshafa2020improving}
\begin{equation}
\begin{split}
       f_{\gamma_{\textup{E}} }^{\textrm{SC}}(\gamma)=&\sum_{k=0}^{N_{E}-1}\frac{N_{E}(-1)^k\binom{N_{E}-1}{k}}{\omega_{e}\Lambda_{E}}\exp\left (- \frac{\gamma(k+1)}{\omega_{e}}  \right )\\&\times\left (\frac{ 1+\frac{\gamma\left ( k+1 \right )}{\omega_{e}}+\frac{1}{\Lambda_{E}} }{\left (\frac{\gamma\left ( k+1 \right )}{\omega_{e}}+\frac{1}{\Lambda_{E}}  \right )^2}\right ).
        \label{o1042}
        \end{split}
\end{equation}
By plugging (\ref{s1021}) and (\ref{o1042}) into (\ref{M018}), and after simple algebraic manipulations, then with the help of \cite[eq. (3.383.9)]{gradshteyn2014table}, the SOP for SC, $\textrm{SOP}_{\textrm{SC}}$, can be derived as
\begin{equation}
\begin{split}
       \textrm{SOP}_{\textrm{SC}}=&N_P\sum_{n=0}^{N_P-1}\:\frac{(-1)^n\:\binom{N_P-1}{n}}{\left ( n+1 \right )}\left [ 1-\frac{N_{E}}{\omega_{e}\:\Lambda_{E}}\right.\\&\left.\times\sum_{k=0}^{N_{E}-1}\:\frac{(-1)^k\:\binom{N_{E}-1}{k}}{\,\mathcal{H}_{1}\exp\left ( -\left (\mathcal{H}_{3}-\mathcal{H}_{2}  \right ) \right )}\right.\\&\left.\left(\frac{\Lambda_{E}}{\exp\left (- \mathcal{H}_{3} \right )} +\frac{\left ( \mathcal{H}_{4}-\mathcal{H}_{1} \right )}{\mathcal{H}_{1}}\Gamma \left ( 0, \mathcal{H}_{3} \right ) \right ) \right ],
        \label{o10420}
        \end{split}
\end{equation}
where $\mathcal{H}_{1}=\frac{\left ( k+1 \right )}{\omega_{e}}
,$ $\mathcal{H}_{2}=\frac{\alpha \,(n+1)}{\Phi},$ $\mathcal{H}_{3}=\frac{\mathcal{H}_{4}}{\Lambda_{E}\,\mathcal{H}_{1}},$  $\mathcal{H}_{4}=\frac{\beta\left ( n+1 \right )}{\Phi}+\frac{\left ( k+1 \right )}{\omega_{e}}
,$ and $\Gamma (\cdot ,\cdot )$ denotes the upper incomplete gamma function \cite[eq. (8.350.2)]{gradshteyn2014table}. For the MRC technique, the PDF of $\gamma_{E_{}}$, $f_{\gamma_{E_{}}}^{\textrm{MRC}}(\gamma)$, is given by \cite[eq. (37)]{khoshafa2020improving}
\begin{equation}
\begin{split}
       f_{\gamma_{\textup{E}} }^{\textrm{MRC}}(\gamma)=&\frac{\gamma ^{N_E-1}\:\exp\left (\frac{-\gamma }{\omega_{e}}  \right )}{\Gamma(N_E)\:\omega_{e}^{N_E}\Lambda_{E}}\:\sum_{k=0}^{N_E}\frac{\binom{N_E}{k}\Gamma(k+1) }{\left (\frac{\gamma}{\omega_{e}}+\frac{1}{\Lambda_{E}}  \right )^{k+1}}.
       \label{ss0119}
       \end{split}
\end{equation}
By plugging (\ref{s1021}) and (\ref{ss0119}) into (\ref{M018}), and after simple algebraic manipulations, then with the help of \cite[eq. (3.383.4)]{gradshteyn2014table}, the SOP for MRC, $\textrm{SOP}_{\textrm{MRC}}$, can be derived as
\begin{equation}
\begin{split}
       \textrm{SOP}_{\textrm{MRC}}=&\sum_{n=0}^{N_P-1}\frac{N_P\binom{N_P-1}{n}}{(-1)^{-n}\left ( n+1 \right )}\left [ 1-\sum_{k=0}^{N_E}\frac{\binom{N_E}{k}\Gamma(k+1)}{\Lambda_{E}^{{N_E-k}}}\right.\\&\left.\frac{\mathcal{H}_{5}^{-\left ( \frac{N_E-k}{2} \right )}W_{\frac{-N_E-k}{2}, \frac{-N_E+k+1}{2}}\big ( \mathcal{H}_{5} \big )}{\exp\big ( -\left ( 0.5\,\mathcal{H}_{5}-\mathcal{H}_{2} \right ) \big )} \right ],
       \label{ss011901}
       \end{split}
\end{equation}
where $\mathcal{H}_{5}=\left (\frac{\Phi+\beta\,\omega_{e}\,(n+1)}{\Phi\,\Lambda_{E}}  \right )$, and ${W}_{a,b}(\cdot )$ denotes the Whittaker function \cite[eq. (9.220.4)]{gradshteyn2014table}.

\subsection{Asymptotic SOP Analysis}
The PN asymptotic SOP, $\textup{SOP}^{ \infty }$, is studied {when $\omega_{p}\rightarrow \infty$}.  In this scenario, we consider that $ \omega_{p}   >  >    \omega_{e}$. $\textup{SOP}^{ \infty }$ is given by
\begin{equation} \textup{SOP}^{ \infty }=\left (\mathcal{G}_{a}^{} \overline{ \gamma }_{d} \right )^{-\mathcal{G}_{d}}+\mathcal{O}(\overline{ \gamma }_{d}^{\:-\mathcal{G}_{d}} ),\end{equation}
\noindent
where $\mathcal{G}_d$ is the secrecy diversity order, $\mathcal{G}_a$ is the secrecy array gain, and $\mathcal{O}(.)$ is the higher order terms. In this respect, the $\textup{SOP}^{ \infty }$ can be derived by first obtaining the asymptotic CDF, $F^{ \infty }_{\gamma_{\textup{P}} }(\gamma)$, \cite[eq. (42)]{khoshafa2020physical}. Then, by plugging $F^{ \infty }_{\gamma_{\textup{P}} }(\gamma)$ into (\ref{M018}), and using \cite[eq. (3.382.4)]{gradshteyn2014table} the $\textup{SOP}^{ \infty }$ can be derived. For the SC technique, $\mathcal{G}_{d}^{\textrm{SC}}$ = $N_P$ and $\mathcal{G}_{a
}^{\textrm{SC}}$ is given by
\begin{equation}
    \begin{split}
       \mathcal{G}_{a}^{\textrm{SC}}=&\left [\sum_{k=0}^{N_{E}-1}\sum_{n=0}^{N_{P}}\frac{\binom{N_{E}-1}{k}\binom{N_{P}}{n}\mathcal{Z}_{1}}{(-1)^{-k}\left ( k+1 \right )^{n}}\left (W_{\frac{-n-1}{2}, \frac{-n}{2}}\left ( \frac{1}{\Lambda_{E}} \right ) \right.\right.\\&\left.\left.\times\frac{1}{\sqrt{\Lambda_{E}}}+W_{\frac{-n-2}{2}, \frac{-i+1}{2}}\left (\frac{1}{\Lambda_{E}} \right ) \right )  \right ]^{\frac{-1}{N_P}},
    \end{split}
    \end{equation}
where $\mathcal{Z}_{1}=\frac{N_{E}\,\beta^{n}\,\alpha^{N_P-n}\Gamma \left ( n+1 \right )\,\omega_{e}^{n}   }{\vartheta ^{N_P}\exp\left ( \frac{-1}{2\omega_{e} } \right )}$. For the MRC technique, $\mathcal{G}_{d}^{\textrm{MRC}}$ = $N_P$ and $\mathcal{G}_{a
}^{\textrm{MRC}}$ is given by
\begin{equation}
    \begin{split}
       \mathcal{G}_{a}^{\textrm{MRC}}=&\left [\sum_{k=0}^{N_{E}}\sum_{n=0}^{N_{P}}\:\binom{N_{E}}{k}\binom{N_{P}}{n}\mathcal{Z}_{2}\right.\\&\left.\times W_{\frac{-N_E-k-n}{2}, \frac{-N_E+k-n+1}{2}}\left ( \frac{1}{\omega_{e}} \right )   \right ]^{\frac{-1}{N_P}},
    \end{split}
\end{equation}
where $\mathcal{Z}_{2}=\frac{\beta^{n}\,\alpha^{N_P-n}\Gamma \left ( N_E+n \right )\Gamma \left ( k+1 \right )\,\omega_{e}^{n}   }{\Gamma \left ( N_{E} \right )\vartheta ^{N_P}\exp\left ( \frac{-1}{2\omega_{e} } \right )\,\Lambda_{E}^{\frac{N_E+n-k}{2}} }$.

\subsection{Probability of Non-zero Secrecy Capacity}
In this subsection, the requirement for the presence of non-zero secrecy capacity is investigated. It is worth noting that the non-zero secrecy capacity is achieved when $\gamma_{C}>\gamma_E$. From (\ref{c6}), the PNSC is given by
\begin{equation}
\begin{split}
    \textrm{PNSC}_{{\varsigma  }}=\textrm{Pr}(\mathcal{C}_{\textup{S}}>0)&=\textrm{Pr}\left ( \frac{1+\gamma_{\textup{P}}}{1+\gamma_{\textup{E}}}> 1 \right )
\\&=1-\int_{0}^{\infty }F_{\gamma _{\textup{P}}}(\gamma)f_{\gamma _\textup{E}}^{\varsigma  }(\gamma)\:d\gamma.
    \label{o1}
\end{split}
  \end{equation}
\subsubsection{Eavesdropper's Channel with SC}
By plugging (\ref{s1021}) and (\ref{o1042}) into (\ref{o1}), and after simple algebraic manipulations, then with the help of \cite[eq. (3.383.9)]{gradshteyn2014table}, the PNSC for SC, $\textrm{PNSC}_{\textrm{SC}}$, can be derived as
\begin{equation}
\begin{split}
       \textrm{PNSC}_{\textrm{SC}}=&1-N_P\sum_{n=0}^{N_P-1}\:\frac{(-1)^n\:\binom{N_P-1}{n}}{\left ( n+1 \right )}\left [ 1-\frac{N_{E}}{\omega_{e}\:\Lambda_{E}}\right.\\&\left.\times\sum_{k=0}^{N_{E}-1}\:\frac{(-1)^k\:\binom{N_{E}-1}{k}}{\,\mathcal{H}_{1}\exp\left ( -\mathcal{H}_{6}  \right )}\right.\\&\left.\times\left(\frac{\Lambda_{E}}{\exp\left (- \mathcal{H}_{6} \right )} +\frac{\left ( \mathcal{H}_{7}-\mathcal{H}_{1} \right )}{\mathcal{H}_{1}}\Gamma \left ( 0, \mathcal{H}_{6} \right ) \right ) \right ],
        \label{o104201}
        \end{split}
\end{equation}
where  $\mathcal{H}_{6}=\frac{\mathcal{H}_{7}}{\Lambda_{E}\,\mathcal{H}_{1}},$  $\mathcal{H}_{7}=\frac{\left ( n+1 \right )}{\Phi}+\frac{\left ( k+1 \right )}{\omega_{e}}
.$

\subsubsection{Eavesdropper's Channel with MRC}
By plugging (\ref{s1021}) and (\ref{ss0119}) into (\ref{o1}), and after simple algebraic manipulations, then with the help of \cite[eq. (3.383.4)]{gradshteyn2014table}, the PNSC for MRC, $\textrm{PNSC}_{\textrm{MRC}}$, can be derived as
\begin{equation}
\begin{split}
       \textrm{PNSC}_{\textrm{MRC}}=&1-\sum_{n=0}^{N_P-1}\frac{N_P\binom{N_P-1}{n}}{(-1)^{-n}\left ( n+1 \right )}\left [ 1-\sum_{k=0}^{N_E}\frac{\binom{N_E}{k}}{\Lambda_{E}^{{N_E-k}}}\right.\\&\left.\times\frac{\Gamma(k+1)\mathcal{H}_{8}^{-\left ( \frac{N_E-k}{2} \right )}W_{\frac{-N_E-k}{2}, \frac{-N_E+k+1}{2}}\big ( \mathcal{H}_{8} \big )}{\exp\big ( - 0.5\,\mathcal{H}_{8} \big )} \right ],
       \label{ss01190}
       \end{split}
\end{equation}
where $\mathcal{H}_{8}=\left (\frac{\Phi+\omega_{e}\,(n+1)}{\Phi\,\Lambda_{E}}  \right )$, and ${W}_{a,b}(\cdot )$ denotes the Whittaker function \cite[eq. (9.220.4)]{gradshteyn2014table}.

\subsection{SN Outage Probability}
For the SN, the outage probability, $P_{out}$, can be expressed by

\begin{equation}    P_{out}=\textup{Pr}\left ( \gamma_{\textup{D}}\leq 2^{\mathcal{R}_d}-1   \right )=F_{\gamma_{\textup{D}}}(2^{\mathcal{R}_d}-1),    \label{R01}\end{equation}
where $\mathcal{R}_d$ is the SN achievable data rate, and $\gamma_{D}$ is the instantaneous SNR of the CR link.
Now, by replacing $P_\textup{S}$ with $\frac{\mathcal{Q}}{\Psi_\textup{P}}$ in (\ref{S11}), $\gamma_{\textup{D}} $ can be written as
\begin{equation}    \begin{split}       \gamma_{\textup{D}} =\frac{\mathcal{Q}}{\Psi_\textup{P}\,\sigma_d^2}\left |\left ( \frac{d_{sr}  d_{rd}}{d_{o}^{2}} \right )^{\hspace{-0.2pc}-\frac{\eta}{2}}\hspace{-0.2pc}\sum_{i=1}^{N} h_{s_i} h_{d_i} e^{j\phi_{i}}+  h_{sd}\left ( \frac{d_{sd}}{d_{o}} \right )^{-\frac{\eta}{2}}  \right |^{2},        \label{S1}    \end{split}\end{equation}
which can be rewritten as
$\gamma_{\textup{D}} =\frac{\Psi_\textup{D}}{\Psi_\textup{P}},
 $ where $\Psi_\textup{D}= \left (\Omega_{1} \sum_{i=1}^{N} \left |h_{s_i}  \right |\, \left |h_{d_i}  \right |+\Omega_{2} \left |h_{sd}  \right | \right )^{2},$ $\Omega_{1} = \frac{{\sqrt{\mathcal{Q}}}}{\sigma_d}\left ( \frac{d_{sr}\,  d_{rd}}{d_{o}^{2}} \right )^{-\frac{\eta}{2}},$ $ \Omega_{2} = \frac{{\sqrt{\mathcal{Q}}}}{\sigma_d}\left ( \frac{d_{sd}}{d_{o}^{2}} \right )^{-\frac{\eta}{2}}
.$ According to the central limit theorem, $\chi _{1}=\sum_{i=1}^{N} \left |h_{s_i}  \right |\, \left |h_{d_i}  \right |$ can be approximated as a Gaussian RV with a mean value $\varepsilon= \frac{N\pi}{4}$ and variance $\sigma^2=N\left ( 1-\frac{\pi}{16} \right )$ \cite{basar2019wireless}. Moreover, $\chi _{2}=\left |h_{sd}  \right |$ is a Rayleigh-distributed RV with a parameter $\delta  $. Thus, the pdfs of $\chi_{1}$ and $\chi_{2}$ are given by
\begin{equation}
    f_{\chi_{1}}(\gamma)=\frac{1}{\sqrt{2\pi\sigma^{2} }}\exp\left ( \frac{-\left ( \gamma -\mu  \right )^{2}}{2\sigma^{2}} \right ),
\end{equation}
and
\begin{equation}
    f_{\chi_{2}}\left (\gamma  \right )=\frac{\gamma }{\delta }\exp\left ( -\frac{\gamma ^{2}}{2\delta } \right ),
\end{equation}
respectively. Thus, $\Psi_\textup{D}$ can be expressed as $\Psi_\textup{D} =\Big (\Omega_{1}  \chi _{1}+\Omega_{2}\,\chi _{2} \Big )^{2},$ leading to the cumulative distribution function (CDF) given by
\cite{yang2020secrecy}
\begin{equation}
    \begin{split}
        &F_{\Psi_\textup{D} }(\gamma )=0.5\left [ \textup{erf}\left ( \frac{(\varphi_1/\varphi_2 )\sqrt{\gamma }-\varepsilon}{\sqrt{2\sigma ^{2}}} \right )+\textup{erf}\left ( \frac{\varepsilon }{\sqrt{2\sigma ^{2}}} \right ) \right ]\\&-\frac{\sqrt{\delta   }}{2\xi_1 }\exp\left ( \frac{-\left (\varphi_1\sqrt{\gamma }-\varphi_2\varepsilon   \right )^{2}}{2\xi _{1}^{2}} \right )\textup{erf}\left ( \frac{\xi _{4}\sqrt{\gamma }-\xi _{5}}{\xi _{1}\xi _{2}} \right )\\&-\frac{\sqrt{\delta   }}{2\xi_1 }\exp\left ( \frac{-\left (\varphi_1\sqrt{\gamma }-\varphi_2\varepsilon   \right )^{2}}{2\xi _{1}^{2}} \right )\textup{erf}\left ( \frac{\xi _{3}\sqrt{\gamma }+\xi _{5}}{\xi _{1}\xi _{2}} \right ),
        \label{Fs}
    \end{split}
\end{equation}
where $\delta$ is a Rayleigh-distributed RV parameter,  $\varepsilon= \frac{N\pi}{4}$, $\sigma^2=N\left ( 1-\frac{\pi}{16} \right )$,   $\varphi_1=\frac{1}{\Omega_{2}}$, $\varphi_2=\frac{\Omega_{1}}{\Omega_{2}}$, $\xi_1=\sqrt{\sigma^{2} \varphi_{2}+\delta    }
$, $\xi_2=\sqrt{2\sigma^{2}\delta  }
$, $\xi_3=\sigma^{2}\varphi_{1} \varphi_{2}
$, $\xi_4=\frac{\delta   \varphi_{1}}{\varphi_{2}}$, $\xi_5=\delta\:   \varepsilon $, and $\textup{erf}(\cdot)$ is the error
function \cite[eq. (8.250.1)]{gradshteyn2014table}. $P_{out}$ can be further written mathematically as \cite{papoulis2002probability}
\begin{equation}
    P_{out}=\int_{0}^{\infty }F_{\Psi_\textup{D}}(\gamma\,x)\,f_{\Psi_\textup{P}}(x)\:dx.
    \label{R02}
\end{equation}
However, utilizing (\ref{Fs}) to derive $P_{out}$  is not mathematically tractable. Therefore, the below approximation of the $\textup{erf}(\cdot)$ function is utilized \cite{sadhwani2017tighter}
\begin{equation}
    \textup{erf}\left ( x \right )\approx \left\{\begin{matrix}
1-\sum_{m=1}^{4}\Upsilon_m \exp\left ( -\Theta_m\, x^{2} \right ) & x\geq 0\\
 -1+\sum_{m=1}^{4}\Upsilon_m \exp\left ( -\Theta _m\, x^{2} \right )& x< 0,
\end{matrix}\right.
\label{App}
\end{equation}
where $\Theta  = [1, 2, 20/ 3 , 20 /17]$, and $\Upsilon = [1/ 8, 1/4, 1/ 4, 1/ 4]$.

\newpage
\begin{strip}
\begin{equation}
    \begin{split}
      \mathcal{A}_{1}=&\exp\left ( \frac{-b_{2}^{2}}{b_{1}^{2}\lambda_{P} } \right )+\frac{1}{2}\left ( \textup{erf}\left ( \frac{\varepsilon }
{\sqrt{2\sigma^{2} }} \right )-1 \right )+
\sum_{m=1}^{4}\frac{\Upsilon_m\exp\left ( -\frac{a_{2}}{2}  \right )}{4b_{1}^{2}\sqrt{(a_{2}/2)^{3}}\lambda_{P} }\left [  2\sqrt{a_{1}}\left ( \exp\left ( b_{2} \left ( a_{2}-a_{1} b_{2}\right )\right )-2 \right )+\exp\left ( \frac{a_{2}^{2}}{4a_{1}} \right )\right.\\&\left. \sqrt{\pi}\left (  a_{2}-2a_{1} b_{2}\right )\left ( \textup{erfc}\left ( \frac{a_{2}}{2\sqrt{a_{1}}} \right )-\textup{erf}\left ( \frac{a_{2}}{2\sqrt{a_{1}}} \right ) +\textup{erf}\left ( \frac{a_{2}-2a_{1}b_{2}}{2\sqrt{a_{1}}} \right ) \right )\right ],
    \end{split}
\end{equation}
\begin{align}
    \begin{split}
      \mathcal{A}_{2}=&\frac{\exp\left ( -c_{3} \right )\Xi_1}{4}\left [\frac{1}{\sqrt{c_{1}^{3}}}\left ( 2\sqrt{c_{1}}\left ( 2-\exp\left ( \left ( \frac{\xi_{5} }{\xi_{1} \xi_{2} } \right )\left ( c_{2}-c_{1}\left ( \frac{\xi_{5} }{\xi_{1} \xi_{2} } \right ) \right )\right ) \right )-\Xi_2\sqrt{\pi}\exp\left ( \frac{c_{2}^{2}}{4c_{1}} \right ) \left (\textup{erf}\left ( \frac{\Xi_2}{2\sqrt{c_{1}}} \right ) \right.\right.\right.\\&\left.\left.\left.-\textup{erf}\left ( \frac{c_{2}}{2\sqrt{c_{1}}} \right )+\textup{erfc}\left ( \frac{c_{2}}{2\sqrt{c_{1}}} \right )  \right )\right )+
\sum_{m=1}^{4}\frac{\Upsilon_m}{\sqrt{c_{4}^{3}}}\left ( 2\sqrt{c_{4}}\left ( \exp\left ( \left ( \frac{\xi_{5} }{\xi_{1} \xi_{2} } \right )\left ( c_{2}-c_{4}\left ( \frac{\xi_{5} }{\xi_{1} \xi_{2} } \right ) \right )\right )-2\right )\right.\right.\\&\left.\left. +\Xi_3\sqrt{\pi }\exp\left ( \frac{c_{2}^{2}}{4c_{4}} \right )\left ( \textup{erf}\left ( \frac{\Xi_3}{2\sqrt{c_{4}}} \right ) -\textup{erf}\left ( \frac{c_{2}}{2\sqrt{c_{4}}} \right )+\textup{erfc}\left ( \frac{c_{2}}{2\sqrt{c_{4}}} \right ) \right ) \right )  \right ],
\label{Pout}
    \end{split}
\end{align}
\begin{align}
    \begin{split}
      \mathcal{A}_{3}=&\frac{\Xi_{4}}{4}\exp\left ( -\left ( \upsilon _3-\left ( \frac{\xi_{5} }{\xi_{1} \xi_{2} } \right ) \upsilon _2 \right ) \right )\left [\upsilon _{1}^{-\frac{3}{2}} \exp\left ( -\left ( \frac{\xi_{5} }{\xi_{1} \xi_{2} } \right )^2 \upsilon _1 \right )
\left ( 2\sqrt{\upsilon _{1}}-\Xi_{5}\sqrt{\pi}\;\exp\left ( \frac{\Xi_{5}^{2}}{4\upsilon _1} \right )\textup{erf}\left ( \frac{\Xi_{5}}{2\sqrt{\upsilon _1}} \right ) \right )\right.\\&\left.+
\sum_{m=1}^{4}\frac{\Upsilon_m}{\sqrt{\upsilon _{4}^{3}}}\exp\left ( -\left ( \frac{\xi_{5} }{\xi_{1} \xi_{2} } \right )^2 \upsilon _4 \right )\left (-2\sqrt{\upsilon _{4}}+\Xi_{6}\sqrt{\pi}\,\exp\left ( \frac{\Xi_{6}^{2}}{4\upsilon _4} \right )\textup{erf}\left ( \frac{\Xi_{6} }{2\sqrt{\upsilon _4}} \right )   \right )  \right ],
\label{Poutt}
    \end{split}
\end{align}
where
\begin{equation*}
  b_{1}=\frac{ \varphi_
1   \sqrt{\gamma }}{\sqrt{2\sigma ^{2}}},\:\:\:\:  b_{2}=\frac{\varepsilon }{\sqrt{2\sigma ^{2}}},\:\:\:\: a_1=\frac{1}{\lambda _{p}\,b_{1}^{2}}+\Theta _m,\:\:\:\:  a_{2}=\frac{2b_{2}}{\lambda _{p}\,b_{1}^{2}},\:\:\:\:\Xi_{1}=\frac{\xi_{1}\xi_{2}^{2} \sqrt{\delta   } }{\lambda_{P}\xi_{4}^{2}  },\:\:\:\:\Xi_{2}=c_{2}-2c_{1}\left ( \frac{\xi_{5} }{\xi_{1}\xi_{2}} \right ),\:\:\:\:
\Xi_{3}=c_{2}-2c_{4}\left ( \frac{\xi_{5} }{\xi_{1}\xi_{2}} \right ),\end{equation*}
\begin{equation*}
c_{1}=\frac{\varphi_{1}^{2}\xi_{2}^{2}  }{2\xi_{4}^{2}}+\frac{\xi_{1}^{2} \xi_{2}^{2} }{\lambda _{P}\xi_{4}^{2}\gamma },\:\:\:\:
c_{2}=\frac{\xi_{2}}{\xi_{1}\xi_{4}^{2}\gamma }\left ( \varphi_{1}\left ( \varphi _{1}\xi_{5} \gamma -\varepsilon \xi_{4}\gamma \right ) +\frac{2\xi_{1}^{2}\xi_{5}}{\lambda _{P}}\right ),\:\:\:\: c_{3}=\frac{\xi_{5}^{2}}{\lambda _{P}\xi_{4}^{2}}+\frac{1 }{2\xi_{1}^{2}}\left ( \frac{\varphi_{1}\xi_{5}}{\xi_{4}} -\varepsilon\right )^{2},\:\:\:\:c_{4}=c_{1}+\Theta _{m},
\end{equation*}

\begin{equation*}
\upsilon _{1}=\frac{\varphi_{1}^{2}\xi_{2}^{2} \gamma  }{2\xi_{3}^{2}}+\frac{\xi_{1}^{2} \xi_{2}^{2} }{\lambda _{P}\xi_{3}^{2} },\:\:\:\:
\upsilon _{2}=\frac{\xi_{2}}{\xi_{1}\xi_{3}^{2} }\left ( \varphi_{1}\left ( \varphi _{1}\xi_{5} \gamma +\varepsilon \xi_{3}\sqrt{\gamma} \right ) +\frac{2\xi_{1}^{2}\xi_{5}}{\lambda _{P}}\right ),\:\:\:\: \upsilon _{3}=\frac{\xi_{5}^{2}}{\lambda _{P}\xi_{3}^{2}}+\frac{1 }{2\xi_{1}^{2}}\left ( \frac{\varphi_{1}\xi_{5}\sqrt{\gamma }}{\xi_{3}} +\varepsilon\right )^{2},\:\:\:\:\upsilon _{4}=\upsilon _{1}+\Theta _{m},
\end{equation*}
\begin{equation*}
\Xi_{4}=\frac{\xi_{1}\xi_{2}^{2} \sqrt{\delta   } }{\lambda_{P}\xi_{3}^{2}},\:\:\:\:\Xi_{5}=2\upsilon _{1}\left ( \frac{\xi_{5} }{\xi_{1}\xi_{2}} \right )-\upsilon _{2},\:\:\:\:\textup{and}\:\:\Xi_{6}=2\upsilon _{4}\left ( \frac{\xi_{5} }{\xi_{1}\xi_{2}} \right )-\upsilon _{2}.
\end{equation*}
\noindent\rule{\textwidth}{0.4pt}
\end{strip}

\noindent

\begin{figure}[t]
\centering
\includegraphics[width=\linewidth]{Figure02.eps}
  \caption{The PN's SOP vs. ${\omega}_{p}$, for different values of the number of reflecting elements, $N$, where $N_P=N_E=3$.}
 \label{fig122}
\end{figure}

\begin{figure}[t]
\centering
\includegraphics[width=\linewidth]{Figure03.eps}
  \caption{The PN's PNSC vs. ${\omega}_{p}$, for different values of the number of reflecting elements, $N$, where $N_P=N_E=3$.}
 \label{fig13}
\end{figure}

\begin{figure}[t]
\centering
\includegraphics[width=\linewidth]{Figure04.eps}
  \caption{The SN's $P_{out}$ vs. $\mathcal{Q}$, for different values of the number of reflecting elements, $N$, where $\mathcal{R}_d = 1$ b/s/Hz.}
  \label{fig11}
\end{figure}

\begin{figure}[t]
\centering
\includegraphics[width=\linewidth]{Figure05.eps}
  \caption{The SN's $P_{out}$ vs. $\mathcal{Q}$, for different
    scenarios, where $N= 30$, $\mathcal{R}_d = 1$ b/s/Hz. }
  \label{fig12}
\end{figure}
\noindent
By substituting (\ref{ccp}) and (\ref{Fs}) to (\ref{R02}) using (\ref{App}), then with the help of \cite[eq. (2.33.1)]{gradshteyn2014table}, $P_{out}$ can be obtained as
\begin{equation}   P_{out}= \mathcal{A}_{1}-\mathcal{A}_{2}-\mathcal{A}_{3},  \label{I1}\end{equation}
where $\mathcal{A}_{1},$ $\mathcal{A}_{2},$ and $\mathcal{A}_{3}$ are given at the top of this page.
\section{Results and Discussions}
\label{s3}
This section provides numerical and simulation results to confirm
the benefits of applying the RIS technology in the proposed system model. Unless otherwise stated, we set $\mathcal{R}_b$  = 1 b/s/Hz, $\mathcal{Q}$ = 10 dBW,  $\bar{\gamma}_{se}$  = 5 dB,  $\delta$ = 2, and $\mathcal{R}_s$  = 1 b/s/Hz.

In Fig. \ref{fig122}, we investigate the PN secrecy enhancement due to the deployment of the RIS technology. In this respect, the SOP of the PN is evaluated for the SC and MRC techniques versus  $\omega_p$, at Eav, for different values of $N$, where $\omega_e$ = 10 dB. The PN secrecy performance is improved as $N$ increases, showing the effect of the RIS's jamming signals toward Eav. Consequently, the PLS of the PN is increased. Moreover, the $\textrm{SOP}$ is enhanced as $\omega_p$ increases.
As revealed in our analysis and simulation, improved secrecy performance can be achieved using RIS as a friendly jammer. This is because the wiretapped  signal is degraded at Eav due to the jamming signals generated by the RIS, resulting in a more secure PN transmission. Since the MRC technique produces a higher SNR gain at Eav over the SC technique, the PN secrecy performance is degraded when Eav utilizes the MRC technique, as illustrated in Fig. \ref{fig122}. The asymptotic analyses are included, and perfect agreement with the theoretical results can be seen when $ \omega_p \rightarrow \infty $, confirming the preciseness of the asymptotic expressions. Finally, it is evident that theoretical and simulation results have an excellent match, confirming the exactness of the derived expressions.

Figure \ref{fig13} plots the PNSC versus $\bar{\gamma}_{p}$. It can be noted that the PNSC improves as $\bar{\gamma}_{p}$ increases for a fixed $\bar{\gamma}_e$. Moreover, the PNSC improves with decreasing $\bar{\gamma}_e$. Further, it is also remarkable that the PNSC increases as $N$ increases. Interestingly, secure transmission is guaranteed as $N$ increases. As expected, for the SC technique, the PNSC is lower than that of the MRC technique. Analytical results are also found to match simulation results, validating the accuracy of our analysis.

The SN outage probability, $P_{out},$ is presented in Fig. \ref{fig11}, where the numerical results are provided and compared with the simulated ones. Towards this end, the effect of $N$ on the RIS is evaluated. As shown in this figure, $P_{out}$ of the SN transmission decreases dramatically when $\mathcal{Q}$ increases. With this in mind, the reliability of SN communication increases as $N$ increases. As an illustration, $\mathcal{Q}$ decreases by nearly 4 dB, deploying an RIS technology with $N$ = 20 compared with $N$ = 50 to reach $P_{out}$ = $10^{-2}$.

In Fig. \ref{fig12}, the reliability of the proposed system model is studied and compared with different scenarios. Towards this end, the relay-aided transmission \cite{9060897}, phase shift error, and unavailability of the line of sight between S and D  scenarios are introduced and the results obtained through Monte-Carlo simulations. To evaluate the influence of the discrete phase shifts, simulation results where the phase error is uniformly distributed in $[-\frac{\pi}{4}, \frac{\pi}{4}]$ \cite{9291402}  are provided. Interestingly, the SN's reliability is enhanced by utilizing the RIS in the presence of the S-D link compared to other scenarios. This is due to the fact that the RIS can maximize the received SNR at D and thus improve the channel quality of the SN. It is also noteworthy that simulation and numerical
results match impeccably, verifying the correctness of our analysis. Furthermore, theoretical  results and simulation results agree perfectly, verifying the exactness of our analysis.

\section{Conclusion}
\label{s4}
In this work, the RIS technology is employed to simultaneously assist SN transmission and enhance the PN's secrecy performance in a CR environment. New analytical expressions are provided for the SN's outage probability and the PN's SOP, considering practical combining techniques. The accuracy of the provided expressions is confirmed via extensive Monte-Carlo simulations. Furthermore, the benefits of the proposed system model are verified through numerical and simulation results.
\balance