Performance Study of FSO/THz Dual-Hop System Based on Cognitive Radio and Energy Harvesting System

Abstract

In order to address the problems of low spectrum efficiency in current communication systems and extend the lifetime of energy-constrained relay devices, this paper proposes a novel dual-hop free-space optical (FSO) system that integrates cognitive radio (CR) and energy harvesting (EH). In this system, the source node communicates with two users at the terminal via FSO and terahertz (THz) hard-switching links, as well as a multi-antenna relay for non-orthogonal multiple access (NOMA). There is another link whose relay acts as both the power beacon (PB) in the EH system and the primary network (PN) in the CR system, achieving the double function of auxiliary transmission. In addition, based on the three possible practical working scenarios of the system, three different transmit powers of the relay are distinguished, thus enabling three different working modes of the system. Closed-form expressions are derived for the interruption outage probability per user for these three operating scenarios, considering the Gamma–Gamma distribution for the FSO link, the $\alpha -\mu$ distribution for the THz link, and the Rayleigh fading distribution for the radio frequency (RF) link. Finally, the numerical results show that this novel system can be adapted to various real-world scenarios and possesses unique advantages.

1. Introduction

With the maturation of fifth-generation (5G) and the ongoing development of sixth-generation (6G) systems, free-space optics (FSO) is emerging as a focal point of interest in academia and industry [1,2,3]. Compared with traditional radio frequency (RF), FSO has a larger bandwidth and higher transmission speed [4]. However, nowadays the demands for communication frequency, rate, and bandwidth are all increasing [5], and FSO systems are highly susceptible to atmospheric conditions like atmospheric turbulence and visibility. Therefore, a practical FSO/RF switching link was suggested, with these systems complementing each other, as early as 2001 by Isaac I. Kim [6]. RF serves as a backup link for the FSO system and is able to continue to secure the transmission of communications when the FSO system receives interference. In addition to RF, whose spectrum resources are scarce, terahertz (THz) has become an important complement to FSO communication systems [7,8]. THz is also important for the study of scattering environments [9,10]. The combination of the two can offer high transmission rates, expand spectrum resources, and enhance system reliability when FSO is impacted by atmospheric conditions [11].

Non-orthogonal multiple access (NOMA) has been proposed and applied as a solution to the low spectrum efficiency in communication systems [12]. This technique allows multiple signals from different users to be superimposed on a single transmission, enabling them to share the same portion of the spectrum simultaneously [13]. Marzieh Najaf [14] introduced a system of FSO combined with NOMA, and studied the performance of integrating NOMA into the FSO backhaul system, following the Gamma–Gamma distribution. The study also provided a closed expression for the outage probability (OP). Coincidentally, NOMA has also been introduced into FSO/RF systems to enhance spectrum utilization efficiency. Amina Girdher’s group [15] proposed and investigated the application of NOMA to intelligent reflecting surface (IRS)-assisted FSO/RF links. They derived exact expressions for OP, coverage probability (CP), and throughput (TP) for FSO obeying the Malaga distribution and RF obeying the Nakagami-m distribution. Liu et al. [16] applied NOMA to an FSO/THz dual-hop multiuser system and derived closed expressions for ergodic capacity (EC) and OP. The system achieves better performance with a smaller power allocation coefficient (PAC) mean deviation, confirming the superiority of the NOMA-assisted FSO/THz-RF dual-hop system.

Although the introduction of NOMA can alleviate the issue of spectrum resource scarcity, more efforts in other parts still need to be invested [17]. In this context, cognitive radio (CR) is applied to enhance the efficiency of spectrum utilization. By combining CR with FSO, the secondary network (SN) in the system can fully and efficiently utilize the underutilized spectrum in the primary network (PN) [18,19]. Eylem Erdogan [20] proposed a network that combines FSO/RF and CR, where the FSO link follows the Nakagami-m distribution and the RF link follows the Weibull fading channel. Erdogan provides the OP for the entire system as well as a closed expression for the traversal capacity. On the other hand, Soumen Mondal [19] investigated the system performance of FSO/RF links combined with CR under the Gamma–Gamma distribution with a Rayleigh fading channel model.

Energy harvesting (EH) is considered a solution to prolong the working life of energy storage components in communication networks [21]. Especially in today’s booming Internet era, relay and terminal devices in communication networks are becoming more diverse and miniaturized. Many of these devices are low-power devices with limited battery capacity or are unable to be connected to a power source. By utilizing an energy harvesting system, these devices can collect energy from power beacons (PBs) through relaying to ensure their proper operation [22]. In Ref. [23], Souad Labghough combined EH with FSO/RF links using amplify-and-forward (AF) and derived OPs for Nakagami-m fading channel and Malaga-M atmospheric turbulence channel conditions.

In addition, antenna selection (AS) has effectively improved the system transmission performance in traditional FSO/RF links [24]. Sanayei [25] proposed the use of AS to enhance system performance in communication systems using multiple-input multiple-output (MIMO) technology as early as 2004. This technique has also been implemented in RF, FSO, THz, and FSO/RF systems in recent years [26,27,28,29].

In this study, we investigated a dual-hop system combining FSO/THz hard switching, NOMA, CR, and AS to further improve the spectrum utilization efficiency of the system. Meanwhile, we offer an EH system for low-power terminals where a stabilized power supply is not available to ensure their normal operation within the entire system [30]. The proposed system can be applied in a variety of practical scenarios. For example, in an Internet of Things (IoT) network using an FSO system, the source needs to deliver a signal to two end users employing NOMA and MIMO [31]. The THz enhances the stability of message delivery in the initial hop of message transmission. Since the light source cannot directly transmit the message to the receiving device with an RF antenna, a relay device with an EH system capable of stable operation is added to facilitate the message forwarding. At the same time, this relay device acts as a secondary network of the CR and can utilize the idle state of the PN for message delivery, thereby enhancing the efficiency of spectrum utilization.

From the above papers, it is clear that most studies have focused on systems that combine FSO/RF dual-hop with either CR or EH, individually. To the best of our knowledge, there are no studies that combine both techniques simultaneously. The main contributions of this paper are summarized as follows:

• We propose and study a novel NOMA-based FSO/THz dual-hop system and utilize both CR and EH for the first time to enhance system performance. For the first time, we utilize the transmitting relay R as both the SN in CR and the energy harvesting relay in EH. Additionally, we employ the PU in the CR model as the PB in the EH model. In this case, R utilizes the EH to harvest energy from the PB for replenishment, and subsequently transmits the message to the users after AS in the N antennas.
• Under such a working model, we exemplify three different working modes to suit actual usage. The three operating modes are: PU idle and EH working normally, PU and EH working normally at the same time, and PU very busy or faulty and unable to act as the PB.
• On the basis of these three modes of operation, we assume that there are two users in the system. The FSO link follows the Gamma–Gamma distribution, while the RF link follows the Rayleigh fading channel distribution. Finally, the closed-form expressions for the exact OP of each of the two users in the three different operating modes is obtained.
• In addition, we study and test the impact of the first-hop signal-to-noise ratio, the number of antennas, and the unused power distribution coefficients on the user OP performance in the system proposed in this paper.

The rest of this article is organized as follows. Section 2 presents the system and channel model, deduces the probability density function (PDF) and cumulative distribution function (CDF) of FSO and THz links in the SR link, and the PDF and CDF of RF links in the EH system and RU link. In addition, we classify and discuss three different working scenarios of the system. In Section 3, we analyze the expressions of two users’ respective OPs in three work scenarios. In Section 4, we introduce some meaningful scenarios and carry out numerical analysis for some parameters of the system. In Section 5, we conclude this paper.

2. System and Channel Models

In this paper, we consider the following communication system scenario, as shown in Figure 1. The first-hop SR link involves S serving as the signal source through NOMA-assisted AF transmission, passing through relay R, which is integrated with the EH system. It transmits signals to both users simultaneously. Since a single FSO link is very susceptible to weather or atmospheric turbulence, the THz link was introduced as an alternative transmission option to support the first link. In this process, the relay R also needs to perform photoelectric signal conversion if it receives FSO signals, but not if it receives THz signals. After that, R will first collect the energy sent from the PB and combine it with its own stored energy to send the signal. The final transmission antenna is determined by AS among N antennas. We assume that the entire transmission process takes 2T, with the first and second links each consuming T time.

2.1. SR Link (FSO/THz Link)

In the first link, S needs to signal users $U_1$ and $U_2$ with messages $x_1$ and $x_2$ superimposed on each other. After receiving the signal, the relay R selects the better performing DF between AF and DF for the deployed protocol, and thus continues the transmission. The signal received at R is [32]:

$$y_R=h_{SR}\sqrt{A_1{P}_S}{x}_1+h_{SR}\sqrt{A_2{P}_S}{x}_2+n_{SR},$$

(1)

where $h_{SR}$ represents the channel transmission coefficient from S to R, and $n_{SR}$ is the additive white Gaussian noise (AWGN) with a mean of 0 variance generated by the intermediate information transmission from S to R.

R needs to perform a photoelectric conversion to convert the optical signal into an electrical signal when it receives the optical signal transmitted by S. This step can be skipped if an electrical signal is received. In the NOMA transmission protocol, the signals of different users are sent together in a superposition, so the interference caused by the signals of the remaining users needs to be eliminated. In this scenario, $x_1$ with higher transmit power does not require perfect successive interference cancellation (SIC) processing [33], while $x_2$ needs to perform SIC outside of R. The SIC process involves initially detecting the signal of $x_1$, then eliminating it from the overall signal, and finally decoding the signal of $x_2$. Thus, we can obtain the Signal-to-Interference-Noise Ratio (SINR) of $x_2$ and the Signal-to-Noise Ratio (SNR) of $x_1$, respectively [33]:

$$\begin{array}{cc}\hfill {\gamma }_{1\to 2}^{\left(1\right)}& ={\displaystyle \frac{P_S{h}_{SR}{\eta }_1\zeta A_1^2}{P_S{h}_{SR}{\eta }_1\zeta A_2^2+N_0}}\hfill \\ \hfill & ={\displaystyle \frac{A_1^2{\gamma }_{SR}}{A_2^2{\gamma }_{SR}+1}},\hfill \end{array}$$

(2)

and

$$\begin{array}{cc}\hfill {\gamma }_2^{\left(1\right)}& ={\displaystyle \frac{P_S{h}_{SR}{\eta }_1\zeta A_2^2}{N_0}}\hfill \\ \hfill & =A_2^2{\gamma }_{SR},\hfill \end{array}$$

(3)

where ${\gamma }_{SR}=P_S{h}_{SR}{\eta }_1\zeta /N_0={\overline{\gamma }}_{SR}{h}_{SR}$, ${\overline{\gamma }}_{SR}$ is the average signal-to-noise ratio of the SR link, $\zeta$ is the photoelectric conversion efficiency when the SR link is an FSO signal, and $\zeta$ is the receive conversion efficiency of the electrical signal when the SR link is a THz signal. Next, the PDF and CDF of the SR link need to be calculated. It should be noted that the system designed in this paper uses hard switches for THz communication and FSO communication. The $\alpha -\mu$ distribution is used in terahertz communication links because it can model and describe the complex fading characteristics of the channel flexibly and accurately. Because of its high frequency characteristics, transmission links in the terahertz band often face extremely complex propagation conditions and significant fading effects, which means traditional distributions such as Rayleigh distribution or Nakagami-M distribution may not be able to fully capture its characteristics. Therefore, the PDF and CDF of each of the two different modes of operation must first be obtained. The PDF of the THz link is as follows [34]:

$$\begin{array}{cc}\hfill & f_{|h_T{|}^2}\left({\lambda }_T\right)={\displaystyle \frac{g_T^2{\mu }_T^{\frac{g_T^2}{{\alpha }_T}\frac{g_T^2}{{\lambda }_T^2}-1}}{2{\left(A_{0,T}{h}_{l,T}{h}_{f,T}\right)}^{g_T^2}\Gamma \left({\mu }_T\right)}}\hfill \\ \hfill & \times {\mathbf{G}}_{1,2}^{2,0}\left[{\displaystyle \frac{{\mu }_T}{{\left(A_{0,T}{h}_{l,T}{h}_{T,f}\right)}^{{\alpha }_T}}}{\lambda }_T^{{\displaystyle \frac{{\alpha }_T}{2}}}\left|\begin{array}{c}-;1\\ 0,{\displaystyle \frac{{\alpha }_T{\mu }_T-g_T^2}{{\alpha }_T}};-\end{array}\right.\right].\hfill \end{array}$$

(4)

The CDF of the THz link is [34]:

$$\left.F_{|h_T{|}^2}\left({\lambda }_T\right)={\varphi }_T\times {\mathbf{G}}_{2,3}^{2,1}\left[{\psi }_T{\lambda }_T^{{\displaystyle \frac{{\alpha }_T}{2}}}\left|\begin{array}{c}1;1+{\displaystyle \frac{g_T^2}{{\alpha }_T}}\\ {\displaystyle \frac{g_T^2}{{\alpha }_T}},{\mu }_T;0\end{array}\right.\right.\right].$$

(5)

The PDF of the FSO link is [35]:

$$\begin{array}{cc}\hfill f_{{\left|h_F\right|}^2}\left({\lambda }_F\right)& ={\displaystyle \frac{g_F^2}{2{\lambda }_F\Gamma \left({\alpha }_F\right)\Gamma \left({\beta }_F\right)}}\hfill \\ \hfill & \times \left.\left.{\mathbf{G}}_{1,3}^{3,0}\left[{\displaystyle \frac{{\alpha }_F{\beta }_F}{A_{0,F}{h}_{l,F}}}\right.{\lambda }_F^{{\displaystyle \frac{1}{2}}}\right|\begin{array}{c}-;g_F^2+1\\ g_F^2,{\alpha }_F,{\beta }_F;-\end{array}\right].\hfill \end{array}$$

(6)

The Gamma–Gamma distribution is mainly used to describe signal fading in FSO links, especially for the effects of environmental factors such as atmospheric turbulence and rain and snow. The Gamma–Gamma distribution is used in FSO systems because it is able to more accurately simulate the attenuation of light signals in complex natural environments, which is especially important in long-distance, high-speed communications and dynamic climate conditions. By using GG distribution, FSO system design can be optimized to improve communication quality and reliability. The CDF of the FSO link is [35]:

$$\begin{array}{cc}\hfill & F_{{\left|h_F\right|}^2}\left({\lambda }_F\right)={\varphi }_F\hfill \\ \hfill & \times \left.{\mathbf{G}}_{3,7}^{6,1}\left[{\psi }_F\lambda \left|\begin{array}{c}1;{\displaystyle \frac{g_F^2+1}{2}},{\displaystyle \frac{g_F^2+2}{2}}\\ {\displaystyle \frac{g_F^2}{2}},{\displaystyle \frac{g_F^2+1}{2}},{\displaystyle \frac{{\alpha }_F}{2}},{\displaystyle \frac{{\alpha }_F+1}{2}},{\displaystyle \frac{{\beta }_F}{2}},{\displaystyle \frac{{\beta }_F+1}{2}};0\end{array}\right.\right.\right],\hfill \end{array}$$

(7)

where ${\mathbf{G}}_{p,q}^{m,n}\left[·\right]$ represents MeijerG function [36], ${\varphi }_F=2^{{\alpha }_F+{\beta }_F-3}{g}_F^2/\pi \Gamma \left({\alpha }_F\right)\Gamma \left({\beta }_F\right)$, ${\varphi }_T=g_T^2/{\alpha }_T\Gamma \left({\mu }_T\right)$, ${\psi }_F={\left({\alpha }_F{\beta }_F/4A_{0,F}{h}_{l,F}\right)}^2$, ${\psi }_T={\mu }_T/{\left(A_{0,T}{h}_{l,T}{\widehat{h}}_{f,T}\right)}^{{\alpha }_T}$, and ${\alpha }_F$ and ${\beta }_F$ represent the attenuation coefficients of atmospheric turbulence, while ${\alpha }_T$ and ${\mu }_T$ denote the attenuation parameters of the $\alpha -\mu$ distribution of the THz link. The ratios between the equivalent beam radius and the standard deviation of the pointing error displacement at the receiving antenna for the FSO and THz links are represented by $g_F$ and $g_T$, respectively. Variables ${\mathrm{A}}_{0,T}$ and ${\mathrm{A}}_{0,F}$ are the fractions of the collected power at the center of the detector for the two links. Path losses for the two links are denoted by $h_{l,T}$ and $h_{l,F}$, and $h_{T,f}$ represents the root-mean-square (RMS) value of the fading channel envelope.

In addition, since the SR link is hard switched and THz and FSO are statistically independent of each other in terms of probability, the PDF of the SR link is:

$$P_{{\left|h_{SR}\right|}^2}\left(\lambda \right)=\mathrm{Pr}\left[{\lambda }_F<\lambda ,{\lambda }_T<\lambda \right].$$

(8)

Then, by bringing (5) and (7) into (8), the PDF of the SR link can be obtained as:

$$\begin{array}{cc}\hfill & F_{|h_{SR}{|}^2}\left(\lambda \right)=F_{|h_T{|}^2}\left(\lambda \right)\times F_{|h_F{|}^2}\left(\lambda \right)\hfill \\ \hfill & \left.={\varphi }_T{G}_{2,3}^{2,1}\left[{\psi }_T{\lambda }^{{\displaystyle \frac{{\alpha }_T}{2}}}\left|\begin{array}{c}1;1+{\displaystyle \frac{g_T^2}{{\alpha }_T}}\\ {\displaystyle \frac{g_T^2}{{\alpha }_T}},{\mu }_T;0\end{array}\right.\right.\right]\hfill \\ \hfill & \times \left.{\varphi }_F{\mathbf{G}}_{3,7}^{6,1}\left[{\psi }_F\lambda \left|\begin{array}{c}1;{\displaystyle \frac{g_F^2+1}{2}},{\displaystyle \frac{g_F^2+2}{2}}\\ {\displaystyle \frac{g_F^2}{2}},{\displaystyle \frac{g_F^2+1}{2}},{\displaystyle \frac{{\alpha }_F}{2}},{\displaystyle \frac{{\alpha }_F+1}{2}},{\displaystyle \frac{{\beta }_F}{2}},{\displaystyle \frac{{\beta }_F+1}{2}};0\end{array}\right.\right.\right].\hfill \end{array}$$

(9)

At the same time that the SR link is operating, R also obtains energy from the PB. The expression for the energy gained is:

$$E_r=T{P}_{B\eta }|h_{BR}{|}^2,$$

(10)

where $\eta$ represents the energy conversion efficiency, $P_B$ denotes the transmit power of the PB, and $|h_{BR}{|}^2$ represents the channel coefficients from PB to R. In this article, all RF links are assumed to follow and obey the Rayleigh distribution. Therefore, the channel coefficients $h_{BR}$ of the BR link and the squares of the channel coefficients $h_{RU1}$ and $h_{RU2}$ of the RU link follow an exponential distribution, and their PDFs and CDFs can be expressed as [32]:

$$f_{|h_b{|}^2}\left(x\right)={\displaystyle \frac{1}{g_b}}\exp \left(-{\displaystyle \frac{x}{g_b}}\right),$$

(11)

and

$$F_{\left(\right|h_b{|}^2)}\left(x\right)=1-\exp (-{\displaystyle \frac{x}{g_b}}),$$

(12)

where $b\in \left\{BR,R{U}_1,R{U}_2\right\}$ and $g_b$ is the desired power gain of the link. Let ${\lambda }_b=1/g_b$ and since this system involves AS, the PDF of ${\left|h_{BR}\right|}^2$ for that link for the nth ($n\in \left\{1,2,\dots ,N\right\}$) antenna is finally selected for transmission and can be expressed as [37]:

$$f_{|h_{BR}{|}^2}\left(x\right)={\displaystyle \frac{x^{N-1}{\lambda }_{B{R}_n}^N}{\Gamma \left(N\right)}}\exp \left(-{\lambda }_{B{R}_n}x\right),$$

(13)

and CDF can be expressed as:

$$F_{|h_{BR}{|}^2}\left(x\right)=1-\exp \left(-{\lambda }_{B{R}_n}x\right)\sum _{d=0}^{N-1}{\displaystyle \frac{x^d{\lambda }_{B{R}_n}^d}{d!}}.$$

(14)

2.2. RU Link (RF Link)

In the second link, R needs to choose one of the N antennas as the transmitter antenna and then utilize the gathered energy along with its own energy for transmission. Although AS increases the complexity of the system, it can significantly improve the performance of the communication system and reduce the impact of channel fading, thereby enhancing the signal quality and system capacity. By selecting the optimal antenna, the system can also achieve a good balance between power consumption and hardware complexity. In this scenario, we assume that the statistics of the link performance between each antenna and the user are independent of each other. The AS scheme consists of two stages. The first stage involves calculating the channel coefficients between each of the two users and the various antennas. We assume that $h_n^{min}$ represents the smaller value of the channel coefficients between the nth antenna and the two users:

$$\left|h_n^{min}\right|=\min \left(\left|h_{R_n{U}_1}\right|,\left|h_{R_n{U}_2}\right|\right).$$

(15)

The second step selects the antenna with the largest channel coefficients among all the antennas, and the selected antenna can be denoted as:

$$n^{*}=\arg \underset{N}{\max }\left(\left|h_n^{min}\right|\right).$$

(16)

Thus, we can obtain the PDF of the channel coefficient of R transmitted through this antenna as [38]:

$$f_{|h_{R{U}_1}{|}^2}\left(x\right)={\lambda }_{R{U}_A}N{e}^{-{\lambda }_{R{U}_A}x}{(1-e^{-{\lambda }_{R{U}_A}x})}^{N-1},$$

(17)

and CDF can be expressed as:

$$\begin{array}{cc}\hfill F_{|h_{R{U}_1}{|}^2}\left(x\right)& ={\left[F_{|h_n^{\min }{|}^2}\left(x\right)\right]}^N\hfill \\ \hfill & ={\left(1-e^{-{\lambda }_{R{U}_A}x}\right)}^N.\hfill \end{array}$$

(18)

The PDF of the channel coefficients $h_{R{U}_2}$ transmitted by R to $U_2$ through this antenna is:

$$\begin{array}{cc}\hfill f_{|h_{R{U}_2}{|}^2}\left(x\right)& ={\lambda }_{R_n{U}_1}{\lambda }_{R_n{U}_2}N\sum _{i=0}^{N-1}\left(\begin{array}{c}N-1\\ i\end{array}\right){(-1)}^{N-i-1}\hfill \\ \hfill & \times \left\{{\displaystyle \frac{1}{V}}\left[e^{-{\lambda }_{R_n{U}_1}x}-e^{-\left(N-i\right){\lambda }_{R{U}_A}x}\right]\right.\hfill \\ \hfill & \left.+{\displaystyle \frac{1}{W}}\left[e^{-{\lambda }_{R_n{U}_{2}x}}-e^{-\left(N-i\right){\lambda }_{R{U}_{A}x}}\right]\right\},\hfill \end{array}$$

(19)

and CDF can be expressed as:

$$\begin{array}{cc}\hfill F_{|h_{R{U}_2}{|}^2}\left(x\right)& ={\lambda }_{R_n{U}_1}{\lambda }_{R_n{U}_2}N\sum _{i=0}^{N-1}\left(\begin{array}{c}N-1\\ i\end{array}\right){(-1)}^{N-i-1}\hfill \\ \hfill & \times \left\{{\displaystyle \frac{1}{V}}\left[{\displaystyle \frac{e^{-\left(N-i\right){\lambda }_{R{U}_A}x}}{\left(N-i\right){\lambda }_{R{U}_A}}}-{\displaystyle \frac{e^{-{\lambda }_{R_n{U}_1}x}}{{\lambda }_{R_n{U}_1}}}\right]\right.\hfill \\ \hfill & \left.+{\displaystyle \frac{1}{W}}\left[{\displaystyle \frac{e^{-\left(N-i\right){\lambda }_{R{U}_{A}x}}}{\left(N-i\right){\lambda }_{R{U}_A}}}-{\displaystyle \frac{e^{-\lambda R_n{U}_{2}x}}{{\lambda }_{R_n{U}_2}}}\right]\right\},\hfill \end{array}$$

(20)

where $V=\left(N-i-1\right){\lambda }_{R{U}_A}+{\lambda }_{R_n{U}_2}$ and $W=\left(N-i-1\right){\lambda }_{R{U}_A}+{\lambda }_{R_n{U}_1}$.

When relay R sends information to the user, the second-hop link, as a sub-network of CR, will be constrained by the main network. Therefore, based on the different operating states of the main network, the whole system can be distinguished into three different operating scenarios. In order to make the model clearer and better integrate with the practical application, we assume a remote environment inspection system as the real application of our model.

In the first scenario, we assume that the PN is idle or the CR of the PN is not available. At this point, the PN can provide energy to R as a PB, and at the same time, it can allocate its spectrum resources to RU links. Since the PN is not working at this time, it will not be affected. Therefore, the PN also does not limit the maximum transmit power of R. For example, when a hypothetical remote environmental detection system is deployed in a remote mountainous area where the network is unstable, the weather monitoring sensor node relies entirely on EH. Sensor nodes work by collecting energy at the relay and storing the data locally at the right time. Once the energy storage is sufficient, the equipment can ensure the continuity of monitoring work by regularly transmitting meteorological data. Once the network is unstable or until the base station is repaired and maintained, the CR cannot work normally. At this time, all data transmission of the sensor relies on its own limited transmission frequency band. In this scenario, the energy transmitted by R is the energy $E_R$ received from the PB and the energy $E_{self}$ stored by itself or supplied by its own power supply, which can be expressed as:

$$E_R=T{P}_B\eta |h_{BR}{|}^2+E_{self},$$

(21)

which gives the transmit power $P_R=E_R/T$ of R as [39]:

$$P_R=\left\{\begin{array}{cc}\eta P_B{\left|h_{BR}\right|}^2+P_{self},\hfill & \mathrm{if}{P}_B{\left|h_{BR}\right|}^2\le \Gamma \\ \eta \Gamma +P_{self},\hfill & \mathrm{if}{P}_B{\left|h_{BR}\right|}^2>\Gamma \end{array}\right.$$

(22)

where $\Gamma$ is the saturation threshold of the transmission energy that can be received from the PB and $P_{self}=E_{self}/\Gamma$ is the transmission power that it can provide itself.

In the second scenario, we assume that the main network PN is in normal operation, and R is experiencing low energy or an unstable power supply. In order for R to function properly, the PN can serve as an energy harvesting PB to transmit energy to R. The following section provides an analysis of the working conditions under such circumstances. First, R, acting as the PU of a CR, will utilize the free spectrum resources of the PN to transmit messages to $U_1$ and $U_2$ after selecting the transmission antenna. In the remote meteorological monitoring station, the sensor node EH collects energy from it. At the same time, the node dynamically selects the idle frequency band for data transmission through CR. For example, the sensor is charged by solar energy during the day and relies on radio frequency signals from nearby communication base stations to continue charging at night. The base station at night is also in a relatively idle state, and CR enables the monitoring station to intelligently access the idle frequency band, thereby transmitting real-time meteorological data, which not only reduces the dependence on the battery but also improves the efficiency of spectrum use. Because the transmission power of R cannot affect the normal operation of the PN itself, it is first necessary to obtain the CDF of the R-PN link channel coefficient $h_{RP}$, which can be seen from (12); the CDF of the channel coefficient obeying the Rayleigh fading distribution is [32]:

$$F_{|h_{RP}{|}^2}\left(x\right)=1-\exp \left(-{\displaystyle \frac{x}{g_{RP}}}\right).$$

(23)

The constraints on the conduct can then be obtained as [40]:

$$\begin{array}{c}Pr\left({\left|h_{RP}\right|}^2{P}_R>I\right)\le \delta \hfill \\ \to 1-F_{|h_{RP}{|}^2}\left({\displaystyle \frac{I}{P_R}}\right)\le \delta \hfill \\ \to P_R\le {\displaystyle \frac{I}{F_{|h_{RP}{|}^2}^{-1}\left(1-\delta \right)}},\hfill \end{array}$$

(24)

where $F_{|h_{RP}{|}^2}^{-1}\left(x\right)=-1/\ln \left(1-x\right)$, $g_{RP}$ is the desired power gain of the RP link, and $\delta$ is the saturation threshold of the interference that the PN can withstand when operating normally. Since the energy that R itself can provide to the transmit power is small or unstable in this scenario, we consider $P_{self}<{\displaystyle \frac{-I}{g_{RP}\ln \delta }}$. This means that the transmit power provided by R’s own energy will not affect the normal operation of the PN. The final transmission power of R can be expressed as:

$$P_R=\left\{\begin{array}{cc}\eta P_B|h_{BR}{|}^2+P_{self},\hfill & \mathrm{if}{P}_B{\left|h_{BR}\right|}^2\le \Phi \\ {\displaystyle \frac{I}{|h_{BR}{|}^2}},\hfill & \mathrm{if}{P}_B{\left|h_{BR}\right|}^2>\Phi ,\end{array}\right.$$

(25)

where $\Phi =-I/g_{RP}\ln \delta -P_{self}$.

In the third scenario, we assume that R itself has a sufficiently stable power supply, and the PN network is operational at this time or the energy transmitting system of the PN network fails. For example, the weather monitoring station is equipped with temporary external power supply, such as a connection to the grid or backup battery system, and does not rely on EH at all. The researchers hope to obtain more data faster by providing a temporary stable power source. Sensor nodes use CR to detect and utilize idle frequency bands for data transmission. With CR technology, monitoring stations are able to intelligently access unused frequency bands to transmit meteorological data without interfering with the main network. This not only effectively uses the spectrum resources but also ensures the real-time data transmission. In this case, R will not obtain additional energy from the PN for signaling but can still utilize the free spectrum resources in the PN with CR. Since the transmission power of R in this scenario still needs to be limited without affecting the normal operation of the PN, the final transmission power of R can be obtained as:

$$P_R=\min \left\{P_{self},{\displaystyle \frac{-I}{g_{RP}\ln \delta }}\right\}.$$

(26)

In addition, the signal received by $U_k$ ($k=1,2$) can be expressed as:

$$y_k=\sqrt{P_R}{h}_{R{U}_k}\left(a_1{x}_1+a_2{x}_2\right)+n_k.$$

(27)

3. Performance Analysis

In this section, we will provide closed expressions for the OPs of each of the two users in three different work scenarios. OP is an important metric to measure the reliability of a communication system, which indicates the probability that the system cannot meet the minimum performance requirements at a specific SNR. The study of OP is helpful to optimize the system design, improve its anti-interference ability, and ensure the stability of data transmission and quality of service. It has practical application value especially in dynamic or harsh channel environments.

Based on Section 2, we divide the communication process of the entire system into two links. Consequently, we can determine the OPs of $U_1$ and $U_2$ as:

$$P_1=Pr(P_{1,1}^{\left(1\right)}\cup P_{1,2}^{\left(1\right)}),$$

(28)

and

$$P_2=Pr(P_{2,1}^{\left(1\right)}\cup P_{2,1}^{\left(2\right)}\cup P_{2,2}^{\left(1\right)}\cup P_{2,2}^{\left(2\right)}),$$

(29)

where $P_{k,l}^{\left(j\right)}$ indicates that user k was unable to decode the information from user $j(j\le k\le 2)$ on the lth link ($l=\{1,2\}$, where $l=1$ for SR links and $l=2$ for RU links). Specifically, $P_{1,1}^{\left(1\right)}\stackrel{\mathsf{\Delta }}{=}\left\{{\gamma }_1^{\left(1\right)}<{\gamma }_{th}^{1,1}\right\}$, $P_{2,1}^{\left(1\right)}\stackrel{\mathsf{\Delta }}{=}\left\{{\gamma }_{1\to 2}^{\left(1\right)}<{\gamma }_{th}^{2,1}\right\}$, $P_{2,1}^{\left(2\right)}\stackrel{\mathsf{\Delta }}{=}\left\{{\gamma }_2^{\left(1\right)}<{\gamma }_{th}^{2,1}\right\}$, $P_{1,2}^{\left(1\right)}\stackrel{\mathsf{\Delta }}{=}\left\{{\gamma }_1^{\left(2\right)}<{\gamma }_{th}^{1,2}\right\}$, $P_{2,2}^{\left(1\right)}\stackrel{\mathsf{\Delta }}{=}\left\{{\gamma }_{1\to 2}^{\left(2\right)}<{\gamma }_{th}^{2,2}\right\}$, and $P_{2,2}^{\left(2\right)}\stackrel{\mathsf{\Delta }}{=}\left\{{\gamma }_2^{\left(2\right)}<{\gamma }_{th}^{2,2}\right\}$, where ${\gamma }_{th}^{k,l}$ denotes the threshold of the kth user on the lth link, and furthermore: [16],

$$\begin{array}{c}{\gamma }_1^{\left(1\right)}={\gamma }_{1\to 2}^{\left(1\right)}={\displaystyle \frac{A_1^2{\gamma }_{SR}}{{\gamma }_{SR}{A}_2^2+1}},{\gamma }_2^{\left(1\right)}=A_2^2{\gamma }_{SR},\hfill \\ {\gamma }_1^{\left(2\right)}={\displaystyle \frac{a_1^2{P}_R{\left|h_{R{U}_1}\right|}^2}{P_R{\left|h_{R{U}_1}\right|}^2{a}_2^2+N_0}},{\gamma }_{1\to 2}^{\left(2\right)}={\displaystyle \frac{a_1^2{P}_R{\left|h_{R{U}_2}\right|}^2}{P_R{\left|h_{R{U}_2}\right|}^2{a}_2^2+N_0}}.\hfill \end{array}$$

(30)

For the convenience of subsequent calculations, we also need to rewrite $P_{k,l}^{\left(j\right)}$. Specifically, the first-hop SR link is rewritten as:

$$P_{k,1}^{\left(j\right)}\stackrel{\mathsf{\Delta }}{=}\left\{{\gamma }_{SR}<{\theta }_{k,1}^j\right\}.$$

(31)

The second-hop RU link is then rewritten as:

$$P_{k,2}^{\left(j\right)}\stackrel{\mathsf{\Delta }}{=}\left\{{\gamma }_{R{U}_k}<{\theta }_{k,2}^j\right\},$$

(32)

where

$$\begin{array}{cc}\hfill {\theta }_{1,1}^{\left(1\right)}=& {\displaystyle \frac{{\gamma }_{th}^{1,1}}{A_1^2-{\gamma }_{th}^{1,1}{A}_2^2}},{\theta }_{2,1}^{\left(1\right)}={\displaystyle \frac{{\gamma }_{th}^{2,1}}{A_1^2-{\gamma }_{th}^{2,1}{A}_2^2}},\hfill \\ \hfill {\theta }_{2,1}^{\left(2\right)}=& {\displaystyle \frac{{\gamma }_{th}^{2,1}}{A_2^2}},{\theta }_{1,2}^{\left(1\right)}={\displaystyle \frac{{\gamma }_{th}^{1,2}}{a_1^2-{\gamma }_{th}^{1,2}{a}_2^2}},\hfill \\ \hfill {\theta }_{2,2}^{\left(1\right)}=& {\displaystyle \frac{{\gamma }_{th}^{2,2}}{a_1^2-{\gamma }_{th}^{2,2}{a}_2^2}},{\theta }_{2,2}^{\left(2\right)}={\displaystyle \frac{{\gamma }_{th}^{2,2}}{a_2^2}}.=I_k^1-I_k^2.\hfill \end{array}$$

(33)

Bringing (31) and (32) into (28) and (29) obtains:

$$P_k=Pr({\gamma }_{SR}<{\theta }_{k,1}^{*}\cup {\gamma }_{R{U}_k}<{\theta }_{k,2}^{*}),$$

(34)

where

$$\begin{array}{c}{\theta }_{1,1}^{*}={\displaystyle \frac{{\gamma }_{th}^{1,1}}{A_1^2-{\gamma }_{th}^{1,1}{A}_2^2}},{\theta }_{2,1}^{*}=max\left\{{\displaystyle \frac{{\gamma }_{th}^{2,1}}{A_1^2-{\gamma }_{th}^{2,1}{A}_2^2}},{\displaystyle \frac{{\gamma }_{th}^{2,1}}{A_2^2}}\right\},\\ {\theta }_{1,2}^{*}={\displaystyle \frac{{\gamma }_{th}^{1,2}}{a_1^2-{\gamma }_{th}^{1,2}{a}_2^2}},{\theta }_{2,2}^{*}=max\left\{{\displaystyle \frac{{\gamma }_{th}^{2,2}}{a_1^2-{\gamma }_{th}^{2,2}{a}_2^2}},{\displaystyle \frac{{\gamma }_{th}^{2,2}}{a_2^2}}\right\},\end{array}$$

(35)

and $k\in \{1,2\}$.

Next, we compute the OP for the first work scenario:

$$\begin{array}{cc}\hfill P_k& =1-Pr\left({\gamma }_{SR}\ge {\theta }_{k,1}^{*}\right)Pr\left({\gamma }_{R{U}_k}\ge {\theta }_{k,2}^{*}\right)\hfill \\ \hfill & =1-\left[1-F_{{\gamma }_{SR}}\left({\theta }_{k,1}^{*}\right)\right]\left(M_k^1+M_k^2\right),\hfill \end{array}$$

(36)

where

$$\begin{array}{cc}\hfill M_k^1& =Pr\left({\gamma }_{R{U}_k}\ge {\theta }_{k,2}^{*},|h_{BR}{|}^2\le E\right)\hfill \\ \hfill & =Pr({\displaystyle \frac{D_k}{z}}-{\displaystyle \frac{P_{peak}}{\eta P_B}}\le |h_{BR}{|}^2\le E,{\displaystyle \frac{D_k}{E+\frac{P_{self}}{\eta P_B}}}\le z)\hfill \\ \hfill & ={\int }_{{\displaystyle \frac{D_k}{E_1}}}^{\infty }\left[F_{|h_{BR}{|}^2}\left(E\right)-F_{|h_{BR}{|}^2}\left({\displaystyle \frac{D_k}{z}}-E_2\right)\right]\hfill \\ \hfill & \times f_{|h_{R{U}_k}{|}^2}\left(z\right)dz\hfill \\ \hfill & =I_k^1-I_k^2,\hfill \end{array}$$

(37)

and

$$\begin{array}{cc}\hfill M_k^2& =Pr\left({\gamma }_{R{U}_k}\ge {\theta }_{k,2}^{*},{\left|h_{BR}\right|}^2>E\right)\hfill \\ \hfill & =Pr\left(z\ge F_k,{\left|h_{BR}\right|}^2>E\right)\hfill \\ \hfill & =\left[1-F_{|h_{R{U}_k}{|}^2}\left(F_k\right)\right]\left[1-F_{|h_{BR}{|}^2}\left(E\right)\right],\hfill \end{array}$$

(38)

in which

$$\begin{array}{cc}\hfill E=& {\displaystyle \frac{\Gamma }{P_B}},E_1={\displaystyle \frac{\Gamma }{P_B}}+P_{self},E_2={\displaystyle \frac{P_{peak}}{\eta P_B}},D_1={\displaystyle \frac{N_0{\theta }_{1,2}^{*}}{\eta P_B}},\hfill \\ \hfill F_1=& {\displaystyle \frac{N_0{\theta }_{1,2}^{*}}{\eta \Gamma -P_{self}}},D_2={\displaystyle \frac{N_0{\theta }_{2,2}^{*}}{\eta P_B}},F_2={\displaystyle \frac{N_0{\theta }_{2,2}^{*}}{\eta \Gamma -P_{self}}}.\hfill \end{array}$$

(39)

Then, bringing (17) into (37) yields:

$$\begin{array}{cc}\hfill I_1^1=& {\int }_{{\displaystyle \frac{D_1}{E_1}}}^{\infty }{e}^{{\displaystyle \frac{-{\lambda }_{B{R}_n{D}_1}}{z}}}\sum _{d=0}^{N-1}{\displaystyle \frac{{\left({\lambda }_{B{R}_n}\right)}^d}{d!}}{({\displaystyle \frac{D_1}{z}}-E_2)}^d{f}_{|h_{R{U}_1}{|}^2}\left(z\right)dz\hfill \\ \hfill =& {\lambda }_{R{U}_A}N\sum _{d=0}^{N-1}\sum _{i=0}^{N-1}{\displaystyle \frac{{\left({\lambda }_{B{R}_n}\right)}^d}{d!}}\left(\begin{array}{c}N-1\\ i\end{array}\right){(-1)}^{N-i-1}\hfill \\ \hfill =& {\int }_{{\displaystyle \frac{D_1}{E_1}}}^{\infty }{({\displaystyle \frac{D_1}{z}}-E_2)}^d{e}^{-\left(N-i\right){\lambda }_{R{U}_A}z-{\displaystyle \frac{{\lambda }_{B{R}_n{D}_1}-E_{2}z}{z}}dz}\hfill \\ \hfill =& {\lambda }_{R{U}_A}N\sum _{d=0}^{N-1}\sum _{i=0}^{N-1}{\displaystyle \frac{{\left({\lambda }_{B{R}_n}\right)}^d}{d!}}\left(\begin{array}{c}N-1\\ i\end{array}\right){(-1)}^{N-i-1}\hfill \\ \hfill & \times \sum _{r=0}^d\left(\begin{array}{c}r\\ d\end{array}\right)D_1^{d-r}{e}^{E_2}\hfill \\ \hfill & \times {\int }_{{\displaystyle \frac{D_1}{E_1}}}^{\infty }{z}^{r-d}{e}^{-\left(N-i\right){\lambda }_{R{U}_A}z-{\displaystyle \frac{{\lambda }_{B{R}_n{D}_1}}{z}}dz},\hfill \end{array}$$

(40)

and

$$\begin{array}{cc}\hfill I_1^2& =e^{-{\lambda }_{B{R}_n}E}\sum _{d=0}^{N-1}{\displaystyle \frac{{\left({\lambda }_{B{R}_n}E\right)}^d}{d!}}{\int }_{{\displaystyle \frac{D_1}{E_1}}}^{\infty }{f}_{|h_{R{U}_1}{|}^2}\left(z\right)dz\hfill \\ \hfill & ={\lambda }_{R{U}_A}N\sum _{d=0}^{N-1}\sum _{i=0}^{N-1}{\displaystyle \frac{{\left({\lambda }_{B{R}_n}E\right)}^d}{d!}}\left(\begin{array}{c}N-1\\ i\end{array}\right){(-1)}^{N-i-1}\hfill \\ \hfill & \times {\displaystyle \frac{1}{\left(N-i\right){\lambda }_{R{U}_A}}}{e}^{-{\lambda }_{B{R}_n}E-{\displaystyle \frac{D_1\left(N-i\right){\lambda }_{R{U}_A}}{E_1}}}.\hfill \end{array}$$

(41)

To further solve (40), we let:

$$J_{1,1}={\int }_{{\displaystyle \frac{D_1}{E_1}}}^{\infty }{z}^{r-d}{e}^{-\left(N-i\right){\lambda }_{R{U}_A}z-{\displaystyle \frac{{\lambda }_{B{R}_n{D}_1}}{z}}dz}.$$

(42)

Then, based on [41] and using Gaussian Chebyshev interpolation integrals (where $N_{1,1}$ is the number of terms of the Gaussian Chebyshev polynomial), we obtain:

$$\begin{array}{cc}\hfill J_{1,1}=& {\int }_0^{\infty }{z}^{-d+r}{e}^{-\left(N-i\right){\lambda }_{R{U}_A}z-{\displaystyle \frac{{\lambda }_{B{R}_B{D}_1}}{Z}}dz}\hfill \\ \hfill & -{\int }_0^{{\displaystyle \frac{D_1}{E_1}}}{z}^{-d+r}{e}^{-\left(N-i\right){\lambda }_{R{U}_A}z-{\displaystyle \frac{{\lambda }_{B{R}_n{D}_1}}{z}}dz}\hfill \\ \hfill =& 2{\left({\displaystyle \frac{{\lambda }_{B{R}_n{D}_1}}{Q_{1,1}}}\right)}^{{\displaystyle \frac{1-d+r}{2}}}{K}_{1-d+r}\left(2\sqrt{B_{1,1}}\right)\hfill \\ \hfill & -{\displaystyle \frac{D_1}{2E_1}}{\displaystyle \frac{\pi }{N_{1,1}}}\sum _{n_1=1}^{N_{1,1}}{\left({\tau }_{1,1}\right)}^{-d+r}{e}^{-Q_{1,1}{\tau }_{1,1}-{\displaystyle \frac{{\lambda }_{B{R}_{nD}1}}{{\tau }_{1,1}}}}\hfill \\ \hfill & \times \sqrt{1-{\left(y_{1,1}\right)}^2},\hfill \end{array}$$

(43)

where

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill Q_{1,1}=& \left(N-i\right){\lambda }_{R{U}_A},B_{1,1}={\lambda }_{B{R}_n}{D}_1{Q}_{1,1},\hfill \\ \hfill {\tau }_{1,1}=& {\displaystyle \frac{D_1}{2E}}\left(y_{1,1}+1\right),y_{1,1}=cos\left({\displaystyle \frac{2n_1-1}{2N_{1,1}}}\pi \right).\hfill \end{array}\end{array}$$

(44)

Next, bring (18) into (38) to obtain:

$$\begin{array}{cc}\hfill M_1^2& ={\lambda }_{R{U}_A}N\sum _{d=0}^{N-1}\sum _{i=0}^{N-1}{\displaystyle \frac{{\left({\lambda }_{B{R}_n}E\right)}^d}{d!}}\left(\begin{array}{c}N-1\\ i\end{array}\right){(-1)}^{N-i-1}\hfill \\ \hfill & \times {\displaystyle \frac{1}{\left(N-i\right){\lambda }_{R{U}_A}}}{e}^{-{\lambda }_{B{R}_n}E-\left(N-i\right){\lambda }_{R{U}_A}{F}_1}.\hfill \end{array}$$

(45)

Afterwards, the closed expression for the OP of $U_1$ for the first working scenario can be obtained by substituting (9), (40), (41), and (45) into (36):

$$\begin{array}{cc}\hfill & P_1=1-\left(1-{\varphi }_T{\mathrm{G}}_{2,3}^{2,1}\left[{\psi }_T{\lambda }^{{\displaystyle \frac{{\alpha }_T}{2}}}\left|\begin{array}{c}1;1+{\displaystyle \frac{g_T^2}{{\alpha }_T}}\\ {\displaystyle \frac{g_T^2}{{\alpha }_T}},{\mu }_T;0\end{array}\right.\right.\right]\hfill \\ \hfill & \left.\left.\times {\varphi }_F{\mathbf{G}}_{3,7}^{6,1}\left[{\psi }_F\lambda \left|\begin{array}{c}1;{\displaystyle \frac{g_F^2+1}{2}},{\displaystyle \frac{g_F^2+2}{2}}\\ {\displaystyle \frac{g_F^2}{2}},{\displaystyle \frac{g_F^2+1}{2}},{\displaystyle \frac{{\alpha }_F}{2}},{\displaystyle \frac{{\alpha }_F+1}{2}},{\displaystyle \frac{{\beta }_F}{2}},{\displaystyle \frac{{\beta }_F+1}{2}};0\end{array}\right.\right.\right]\right)\hfill & \\ \hfill & \times \left\{H_1\left[\sum _{r=0}^d\left(\begin{array}{c}r\\ d\end{array}\right)D_1^{d-r}{e}^{E_2}{J}_{1,1}-E^d{\displaystyle \frac{1}{Q_{1,1}}}{e}^{-{\lambda }_{B{R}_n}E}\right.\right.\hfill \\ \hfill & \times \left.\left.\left(e^{-Q_{1,1}{F}_1}-e^{{\displaystyle \frac{-Q_{1,1}{D}_1}{E}}}\right)\right]\right\},\hfill \end{array}$$

(46)

where

$$H_1={\lambda }_{R{U}_A}N\sum _{d=0}^{N-1}\sum _{i=0}^{N-1}{\displaystyle \frac{{\left({\lambda }_{B{R}_n}\right)}^d}{d!}}\left(\begin{array}{c}N-1\\ i\end{array}\right){(-1)}^{N-i-1}.$$

(47)

After that, we compute the OP for $U_2$, which can be obtained as:

$$\begin{array}{cc}\hfill I_2^1& ={\int }_{{\displaystyle \frac{D_2}{E_1}}}^{\infty }{e}^{{\displaystyle \frac{-{\lambda }_{B{R}_n{D}_2}+E_{2}z}{z}}}\sum _{d=0}^{N-1}{\displaystyle \frac{{\left({\lambda }_{B{R}_n}\right)}^d}{d!}}\hfill \\ \hfill & \times {({\displaystyle \frac{D_2}{z}}-E_2)}^d{f}_{|h_{R{U}_1}{|}^2}\left(z\right)dz\hfill \\ \hfill & ={\lambda }_{R_n{U}_1}{\lambda }_{R_n{U}_2}N\sum _{r=0}^d\left(\begin{array}{c}d\\ r\end{array}\right)D_2^{d-r}{e}^{\lambda B_{R_n}{E}_2}\hfill \\ \hfill & \times {(-E_2)}^r\sum _{d=0}^{N-1}\sum _{i=0}^{N-1}{\displaystyle \frac{{\left({\lambda }_{B{R}_n}\right)}^d}{d!}}{(-1)}^{N-i-1}\hfill \\ \hfill & \times \left(\begin{array}{c}N-1\\ i\end{array}\right)\left[{\displaystyle \frac{1}{V}}\left(J_{2,1}-J_{2,2}\right)+{\displaystyle \frac{1}{W}}\left(J_{2,3}-J_{2,2}\right)\right],\hfill \end{array}$$

(48)

and

$$\begin{array}{cc}\hfill I_2^2& ={\int }_{{\displaystyle \frac{D_2}{E_1}}}^{\infty }{e}^{{\displaystyle \frac{-{\lambda }_{B{R}_n{D}_2}+E_{2}z}{z}}}\sum _{d=0}^{N-1}{\displaystyle \frac{{\left({\lambda }_{B{R}_n}\right)}^d}{d!}}\hfill \\ \hfill & \times {({\displaystyle \frac{D_2}{z}}-E_2)}^d{f}_{|h_{R{U}_1}{|}^2}\left(z\right)dz\hfill \\ \hfill & ={\lambda }_{R_n{U}_1}{\lambda }_{R_n{U}_2}N\sum _{r=0}^d\left(\begin{array}{c}d\\ r\end{array}\right)D_2^{d-r}{e}^{\lambda B_{R_n}{E}_2}{(-E_2)}^r\hfill \\ \hfill & \times \sum _{d=0}^{N-1}\sum _{i=0}^{N-1}{\displaystyle \frac{{\left({\lambda }_{B{R}_n}\right)}^d}{d!}}{(-1)}^{N-i-1}\left(\begin{array}{c}N-1\\ i\end{array}\right)\hfill \\ \hfill & \times \left[{\displaystyle \frac{1}{V}}\left(J_{2,1}-J_{2,2}\right)+{\displaystyle \frac{1}{W}}\left(J_{2,3}-J_{2,2}\right)\right],\hfill \end{array}$$

(49)

where

$$\begin{array}{cc}\hfill J_{2,m}& =2{\left({\displaystyle \frac{{\lambda }_{B{R}_n}{D}_2}{Q_{2,m}}}\right)}^{{\displaystyle \frac{1-d+r}{2}}}{K}_{1-d+r}\left(2\sqrt{B_{2,m}}\right)\hfill \\ \hfill & -{\displaystyle \frac{D_2}{2E}}{\displaystyle \frac{\pi }{N_{2,m}}}\sum _{n_m=1}^{N_{2,m}}{\left({\tau }_{2,m}\right)}^{-d+r}\hfill \\ \hfill & \times e^{-Q_{2,m}{\tau }_{2,m}-{\displaystyle \frac{{\lambda }_{B{R}_n{D}_2}}{{\tau }_{2,m}}}}\sqrt{1-{\left(y_{2,m}\right)}^2},\hfill \end{array}$$

(50)

in which $m=\left\{1,2,3\right\}$ and $N_{2,m}$ is the number of terms of the Gauss Chebyshev polynomials, while:

$$\begin{array}{cc}\hfill Q_{2,1}=& {\lambda }_{R_n{U}_1},B_{2,1}={\lambda }_{B{R}_n}{D}_2{Q}_{2,1},\hfill \\ \hfill {\tau }_{2,1}=& {\displaystyle \frac{D_2}{2E_1}}\left(y_{2,1}+1\right),y_{2,1}=\cos \left({\displaystyle \frac{2n_1-1}{2N_{2,1}}}\pi \right),\hfill \\ \hfill Q_{2,2}=& \left(N-i\right){\lambda }_{R{U}_A},B_{2,2}={\lambda }_{B{R}_n}{D}_2{Q}_{2,2},\hfill \\ \hfill {\tau }_{2,2}=& {\displaystyle \frac{D_2}{2E_1}}\left(y_{2,2}+1\right),y_{2,2}=\cos \left({\displaystyle \frac{2n_2-1}{2N_{2,2}}}\pi \right),\hfill \\ \hfill Q_{2,3}=& {\lambda }_{R_n{U}_2},B_{2,3}={\lambda }_{B{R}_n}{D}_2{Q}_{2,3},\hfill \\ \hfill {\tau }_{2,3}=& {\displaystyle \frac{D_2}{2E_1}}\left(y_{2,3}+1\right),y_{2,3}=\cos \left({\displaystyle \frac{2n_3-1}{2N_{2,3}}}\pi \right).\hfill \end{array}$$

(51)

After this, bringing (20) into (38) gives:

$$\begin{array}{cc}\hfill M_{2,2}& ={\lambda }_{R_n{U}_1}{\lambda }_{R_n{U}_2}N\sum _{d=0}^{N-1}\sum _{i=0}^{N-1}{\displaystyle \frac{{\left({\lambda }_{B{R}_n}E\right)}^d}{d!}}\hfill \\ \hfill & \times {(-1)}^{N-i-1}\left(\begin{array}{c}N-1\\ i\end{array}\right)e^{-{\lambda }_{B{R}_n}E}\hfill \\ \hfill & \times \left[{\displaystyle \frac{1}{V}}\left({\displaystyle \frac{1}{Q_{2,1}}}{e}^{-Q_{2,1}{F}_2}-{\displaystyle \frac{1}{Q_{2,2}}}{e}^{-Q_{2,2}{F}_2}\right)\right.\hfill \\ \hfill & \left.+{\displaystyle \frac{1}{W}}\left({\displaystyle \frac{1}{Q_{2,3}}}{e}^{-Q_{2,3}{F}_2}-{\displaystyle \frac{1}{Q_{2,2}}}{e}^{-Q_{2,2}{F}_2}\right)\right].\hfill \end{array}$$

(52)

Eventually, bringing (9), (48), (49), and (52) into (36) gives the OP of $U_2$:

$$\begin{array}{cc}\hfill P_2=& 1-\left(1-{\varphi }_T{\mathrm{G}}_{2,3}^{2,1}\left[{\psi }_T{\lambda }^{{\displaystyle \frac{{\alpha }_T}{2}}}\left|\begin{array}{c}1;1+{\displaystyle \frac{g_T^2}{{\alpha }_T}}\\ {\displaystyle \frac{g_T^2}{{\alpha }_T}},{\mu }_T;0\end{array}\right.\right.\right]\hfill \\ \hfill & \left.\left.\times {\varphi }_F{\mathbf{G}}_{3,7}^{6,1}\left[{\psi }_F\lambda \left|\begin{array}{c}1;{\displaystyle \frac{g_F^2+1}{2}},{\displaystyle \frac{g_F^2+2}{2}}\\ {\displaystyle \frac{g_F^2}{2}},{\displaystyle \frac{g_F^2+1}{2}},{\displaystyle \frac{{\alpha }_F}{2}},{\displaystyle \frac{{\alpha }_F+1}{2}},{\displaystyle \frac{{\beta }_F}{2}},{\displaystyle \frac{{\beta }_F+1}{2}};0\end{array}\right.\right.\right]\right)\hfill \\ \hfill & \times H_2\left\{\left[{\displaystyle \frac{D_2^d}{V}}\left(J_{2,1}-J_{2,2}\right)+{\displaystyle \frac{D_2^d}{W}}\left(J_{2,3}-J_{2,2}\right)\right]\right.\hfill \\ \hfill & -e^{-{\lambda }_{B{R}_n}E}\left\{{\displaystyle \frac{E^d}{V}}\left[{\displaystyle \frac{1}{Q_{2I,1}}}\left(e^{-{\displaystyle \frac{Q_{2,1}{D}_2}{E}}}+e^{-Q_{2,1}{F}_2}\right)\right.\right.\hfill \\ \hfill & \left.-{\displaystyle \frac{1}{Q_{2,2}}}\left(e^{-{\displaystyle \frac{Q_{2,2}{D}_2}{E}}}+e^{-Q_{2,2}{F}_2}\right)\right]\hfill \\ \hfill & +{\displaystyle \frac{E^d}{W}}\left[{\displaystyle \frac{1}{Q_{2,3}}}\left(e^{-{\displaystyle \frac{Q_{2,3}{D}_2}{E}}}+e^{-Q_{2,3}{F}_2}\right)\right.\hfill \\ \hfill & \left.\left.\left.-{\displaystyle \frac{1}{Q_{2,2}}}\left(e^{-{\displaystyle \frac{Q_{2,2}{D}_2}{E}}}+e^{-Q_{2,2}{F}_2}\right)\right]\right\}\right\},\hfill \end{array}$$

(53)

where $H_2={\lambda }_{RU}N{\displaystyle \sum _{d=0}^{N-1}}{\displaystyle \sum _{i=0}^{N-1}}{\displaystyle \frac{{\lambda }_{B{R}_n}^d}{d!}}\left(\begin{array}{c}N-1\\ i\end{array}\right){(-1)}^{N-i-1}$, and ${\lambda }_{RU}={\lambda }_{R_n{U}_1}{\lambda }_{R_n{U}_2}$.

Next is the derivation of the OPs for the two users in the second work scenario, which is given by (17), (19), (22), and (37):

$$\begin{array}{cc}\hfill M_k^{S,1}& =\mathrm{Pr}({\gamma }_{R{U}_k}\ge {\theta }_{k,2}^{*},|h_{BR}{|}^2\le E_3)\hfill \\ \hfill & =\mathrm{Pr}({\displaystyle \frac{D_k}{z}}-{\displaystyle \frac{P_{peak}}{\eta P_B}}\le |h_{BR}{|}^2\le E_3,{\displaystyle \frac{D_k}{E_3+\frac{P_{peak}}{\eta P_B}}}\le z)\hfill \\ \hfill & ={\int }_{{\displaystyle \frac{D_k}{E_4}}}^{\infty }\left[F_{|h_{BR}{|}^2}\left(E_3\right)-F_{|h_{BR}{|}^2}\left({\displaystyle \frac{D_k}{z}}-E_2\right)\right]\hfill \\ \hfill & \times f_{|h_{R{U}_k}{|}^2}\left(z\right)dz\hfill \\ \hfill & =I_k^{S,1}-I_k^{S,2},\hfill \end{array}$$

(54)

where S denotes the second scene; in addition:

$$\begin{array}{cc}\hfill I_1^{S,1}=& {\int }_{{\displaystyle \frac{D_1}{E_4}}}^{\infty }{e}^{{\displaystyle \frac{-{\lambda }_{B{R}_n{D}_1}}{z}}}\sum _{d=0}^{N-1}{\displaystyle \frac{{\left({\lambda }_{B{R}_n}\right)}^d}{d!}}\hfill \\ \hfill & \times {({\displaystyle \frac{D_1}{z}}-E_2)}^d{f}_{|h_{R{U}_1}{|}^2}\left(z\right)dz\hfill \\ \hfill =& {\lambda }_{R{U}_A}N\sum _{d=0}^{N-1}\sum _{i=0}^{N-1}{\displaystyle \frac{{\left({\lambda }_{BR}\right)}^d}{d!}}\left(\begin{array}{c}N-1\\ i\end{array}\right){(-1)}^{N-i-1}\hfill \\ \hfill & \times \sum _{r=0}^d\left(\begin{array}{c}d\\ r\end{array}\right)D_1^{d-r}{e}^{E_2}\hfill \\ \hfill & \times {\int }_{{\displaystyle \frac{D_1}{E_4}}}^{\infty }{z}^{r-d}{e}^{-\left(N-i\right){\lambda }_{R{U}_A}z-{\displaystyle \frac{{\lambda }_{B{R}_n{D}_1}}{z}}dz},\hfill \end{array}$$

(55)

and

$$\begin{array}{cc}\hfill I_1^{S,2}=& e^{-{\lambda }_{B{R}_n}{E}_3}\sum _{d=0}^{N-1}{\displaystyle \frac{{\left({\lambda }_{B{R}_n}{E}_3\right)}^d}{d!}}{\int }_{{\displaystyle \frac{D_1}{E_4}}}^{\infty }{f}_{|h_{R{U}_1}{|}^2}\left(z\right)dz\hfill \\ \hfill =& {\lambda }_{R{U}_A}N\sum _{d=0}^{N-1}\sum _{i=0}^{N-1}{\displaystyle \frac{{\left({\lambda }_{B{R}_n}{E}_3\right)}^d}{d!}}\left(\begin{array}{c}N-1\\ i\end{array}\right)\hfill \\ \hfill & \times {(-1)}^{N-i-1}{\displaystyle \frac{1}{\left(N-i\right){\lambda }_{R{U}_A}}}{e}^{-{\lambda }_{B{R}_n}{E}_3-{\displaystyle \frac{D_1\left(N-i\right){\lambda }_{R{U}_A}}{E_4}}},\hfill \end{array}$$

(56)

where $E_3={\displaystyle \frac{-I}{g_{RP}ln\theta }}$, $E_4={\displaystyle \frac{D_1}{E_2+E_3}}$, $F_3={\displaystyle \frac{N_0{\theta }_{1,2}^{*}}{\frac{-I}{g_{RP}ln\theta }}}$ and $F_4={\displaystyle \frac{N_0{\theta }_{2,2}^{*}}{\frac{-I}{g_{RP}ln\theta }}}$.

Similarly, let:

$$J_{1,1}^S={\int }_{{\displaystyle \frac{D_1}{E_4}}}^{\infty }{z}^{r-d}{e}^{-\left(N-i\right){\lambda }_{R{U}_A}z-{\displaystyle \frac{{\lambda }_{B{R}_n{D}_1}}{z}}dz},$$

(57)

and we can obtain:

$$\begin{array}{cc}\hfill J_{1,1}^S& =2{\left({\displaystyle \frac{{\lambda }_{B{R}_n{D}_1}}{Q_{1,1}}}\right)}^{{\displaystyle \frac{1-d+r}{2}}}{K}_{1-d+r}\left(2\sqrt{B_{1,1}}\right)\hfill \\ \hfill & -{\displaystyle \frac{D_1}{2E_4}}{\displaystyle \frac{\pi }{N_{1,1}^S}}\sum _{n_1=1}^{N_{1,1}^S}{\left({\tau }_{1,1}\right)}^{-d+r}{e}^{-Q_{1,1}{\tau }_{1,1}-{\displaystyle \frac{{\lambda }_{B{R}_n{D}_1}}{{\tau }_{1,1}}}}\hfill \\ \hfill & \times \sqrt{1-{\left(y_{1,1}^S\right)}^2},\hfill \end{array}$$

(58)

where ${\tau }_{1,1}^S=D_1/\left(2E\right)\left(y_{1,1}+1\right)$, $y_{1,1}^S=\cos \left(\left[2n_1-1\left]/\right[2N_{1,1}^S\right]\pi \right)$, and $N_{1,1}^S$ is the number of terms of the Gauss Chebyshev polynomial.

Then, we can bring (18), (20), and (22) into (38) to obtain:

$$\begin{array}{cc}\hfill M_1^{S,2}& ={\lambda }_{R{U}_A}N\sum _{d=0}^{N-1}\sum _{i=0}^{N-1}{\displaystyle \frac{{\left({\lambda }_{B{R}_n}{E}_3\right)}^d}{d!}}\left(\begin{array}{c}N-1\\ i\end{array}\right){(-1)}^{N-i-1}\hfill \\ \hfill & \times {\displaystyle \frac{1}{\left(N-i\right){\lambda }_{R{U}_A}}}{e}^{-{\lambda }_{B{R}_n}E-\left(N-i\right){\lambda }_{R{U}_A}{F}_3}.\hfill \end{array}$$

(59)

By taking (9), (55), (56), and (59) into (36), we obtain the OP of $U_1$ as:

$$\begin{array}{cc}\hfill P_1& =1-\left(1-{\varphi }_T{\mathrm{G}}_{2,3}^{2,1}\left[{\psi }_T{\lambda }^{{\displaystyle \frac{{\alpha }_T}{2}}}\left|\begin{array}{c}1;1+{\displaystyle \frac{g_T^2}{{\alpha }_T}}\\ {\displaystyle \frac{g_T^2}{{\alpha }_T}},{\mu }_T;0\end{array}\right.\right.\right]\hfill \\ \hfill & \left.\left.\times {\varphi }_F{\mathbf{G}}_{3,7}^{6,1}\left[{\psi }_F\lambda \left|\begin{array}{c}1;{\displaystyle \frac{g_F^2+1}{2}},{\displaystyle \frac{g_F^2+2}{2}}\\ {\displaystyle \frac{g_F^2}{2}},{\displaystyle \frac{g_F^2+1}{2}},{\displaystyle \frac{{\alpha }_F}{2}},{\displaystyle \frac{{\alpha }_F+1}{2}},{\displaystyle \frac{{\beta }_F}{2}},{\displaystyle \frac{{\beta }_F+1}{2}};0\end{array}\right.\right.\right]\right)\hfill \\ \hfill & \times \left\{H_1\left[\sum _{r=0}^d\left(\begin{array}{c}r\\ d\end{array}\right)D_1^{d-r}{e}^{E_2}{J}_{1,1}-{E_3}^d{\displaystyle \frac{1}{Q_{1,1}}}\right.\right.\hfill \\ \hfill & \times \left.\left.\left(e^{{\lambda }_{B{R}_n}{E}_3-Q_{1,1}{F}_3}-e^{{\lambda }_{B{R}_n}E-{\displaystyle \frac{Q_{1,1}{D}_1}{E_4}}}\right)\right]\right\}.\hfill \end{array}$$

(60)

Next, calculate the OP of $U_2$, which is obtained by bringing (20) and (25) into (37):

$$\begin{array}{cc}\hfill I_2^{S,1}& ={\int }_{{\displaystyle \frac{D_2}{E_4}}}^{\infty }{e}^{{\displaystyle \frac{-{\lambda }_{B{R}_n{D}_2}+E_{2}z}{z}}}\sum _{d=0}^{N-1}{\displaystyle \frac{{\left({\lambda }_{B{R}_n}\right)}^d}{d!}}\hfill \\ \hfill & \times {({\displaystyle \frac{D_2}{z}}-E_2)}^d{f}_{|h_{R{U}_1}{|}^2}\left(z\right)dz\hfill \\ \hfill & ={\lambda }_{R_n{U}_1}{\lambda }_{R_n{U}_2}N\sum _{r=0}^d\left(\begin{array}{c}d\\ r\end{array}\right)D_2^{d-r}{e}^{\lambda B_{R_n}{E}_2}\hfill \\ \hfill & \times {(-E_2)}^r\sum _{d=0}^{N-1}\sum _{i=0}^{N-1}{\displaystyle \frac{{\left({\lambda }_{B{R}_n}\right)}^d}{d!}}{(-1)}^{N-i-1}\left(\begin{array}{c}N-1\\ i\end{array}\right)\hfill \\ \hfill & \times \left[{\displaystyle \frac{1}{V}}\left(J_{2,1}^S-J_{2,2}^S\right)+{\displaystyle \frac{1}{W}}\left(J_{2,3}^S-J_{2,2}^S\right)\right],\hfill \end{array}$$

(61)

and

$$\begin{array}{cc}\hfill I_2^{S,2}& ={\int }_{{\displaystyle \frac{D_2}{E_4}}}^{\infty }{e}^{-{\lambda }_{B{R}_n}{E}_3}\sum _{d=0}^{N-1}{\displaystyle \frac{{\left({\lambda }_{B{R}_n}E\right)}^d}{d!}}{f}_{|h_{R{U}_2}{|}^2}\left(z\right)dz\hfill \\ \hfill & ={\lambda }_{R_n{U}_1}{\lambda }_{R_n{U}_2}N\sum _{d=0}^{N-1}\sum _{i=0}^{N-1}{\displaystyle \frac{{\left({\lambda }_{B{R}_n}E\right)}^d}{d!}}\hfill \\ \hfill & \times {(-1)}^{N-i-1}\left(\begin{array}{c}N-1\\ i\end{array}\right)e^{-{\lambda }_{B{R}_n}{E}_3}\hfill \\ \hfill & \times \left[{\displaystyle \frac{1}{V}}\left({\displaystyle \frac{1}{Q_{2,1}}}{e}^{{\displaystyle \frac{-Q_{2,1}{D}_2}{E_4}}}-{\displaystyle \frac{1}{Q_{2,2}}}{e}^{{\displaystyle \frac{-Q_{2,2}{D}_2}{E_4}}}\right)\right.\hfill \\ \hfill & \left.+{\displaystyle \frac{1}{W}}\left({\displaystyle \frac{1}{Q_{2,3}}}{e}^{{\displaystyle \frac{-Q_{2,3}{D}_2}{E_4}}}-{\displaystyle \frac{1}{Q_{2,2}}}{e}^{{\displaystyle \frac{-Q_{2,2}{D}_2}{E_4}}}\right)\right],\hfill \end{array}$$

(62)

where

$$\begin{array}{cc}\hfill J_{2,m}^S& =2{\left({\displaystyle \frac{{\lambda }_{B{R}_n}{D}_2}{Q_{2,m}}}\right)}^{{\displaystyle \frac{1-d+r}{2}}}{K}_{1-d+r}\left(2\sqrt{B_{2,m}}\right)\hfill \\ \hfill & -{\displaystyle \frac{D_2}{2E}}{\displaystyle \frac{\pi }{N_{2,m}^S}}\sum _{n_m=1}^{N_{2,m}}{\left({\tau }_{2,m}^S\right)}^{-d+r}\hfill \\ \hfill & \times e^{-Q_{2,m}{\tau }_{2,m}-{\displaystyle \frac{{\lambda }_{B{R}_n{D}_2}}{{\tau }_{2,m}}}}\sqrt{1-{\left(y_{2,m}^S\right)}^2}.\hfill \end{array}$$

(63)

In this equation, $m=\left\{1,2,3\right\}$, $N_{2,m}^S$ represents the number of terms of the Gauss–Chebyshev polynomial, and:

$$\begin{array}{c}{\tau }_{2,1}^S={\displaystyle \frac{D_2}{2E_4}}\left(y_{2,1}^S+1\right),{\tau }_{2,2}^S={\displaystyle \frac{D_2}{2E_4}}\left(y_{2,2}^S+1\right),\\ {\tau }_{2,3}^S={\displaystyle \frac{D_2}{2E_4}}\left(y_{2,3}^S+1\right),y_{2,1}^S=cos\left({\displaystyle \frac{2n_1-1}{2N_{2,1}^S}}\pi \right),\\ y_{2,2}^S=cos\left({\displaystyle \frac{2n_2-1}{2N_{2,2}^S}}\pi \right),y_{2,3}^S=cos\left({\displaystyle \frac{2n_3-1}{2N_{2,3}^S}}\pi \right).\end{array}$$

(64)

Bring (20) and (22) into (38) to obtain:

$$\begin{array}{cc}\hfill M_{2,2}& ={\lambda }_{R_n{U}_1}{\lambda }_{R_n{U}_2}N\sum _{d=0}^{N-1}\sum _{i=0}^{N-1}{\displaystyle \frac{{\left({\lambda }_{B{R}_n}{E}_3\right)}^d}{d!}}\hfill \\ \hfill & \times {(-1)}^{N-i-1}\left(\begin{array}{c}N-1\\ i\end{array}\right)e^{-{\lambda }_{B{R}_n}{E}_3}\hfill \\ \hfill & \times \left[{\displaystyle \frac{1}{V}}\left({\displaystyle \frac{1}{Q_{2,1}}}{e}^{-Q_{2,1}{F}_4}-{\displaystyle \frac{1}{Q_{2,2}}}{e}^{-Q_{2,2}{F}_4}\right)\right.\hfill \\ \hfill & \left.+{\displaystyle \frac{1}{W}}\left({\displaystyle \frac{1}{Q_{2,3}}}{e}^{-Q_{2,3}{F}_4}-{\displaystyle \frac{1}{Q_{2,2}}}{e}^{-Q_{2,2}{F}_4}\right)\right].\hfill \end{array}$$

(65)

Then, we can bring (9), (61), (62), and (65) into (36) to obtain the OP of $U_2$:

$$\begin{array}{cc}\hfill P_2& =1-\left(1-{\varphi }_T{\mathrm{G}}_{2,3}^{2,1}\left[{\psi }_T{\lambda }^{{\displaystyle \frac{{\alpha }_T}{2}}}\left|\begin{array}{c}1;1+{\displaystyle \frac{g_T^2}{{\alpha }_T}}\\ {\displaystyle \frac{g_T^2}{{\alpha }_T}},{\mu }_T;0\end{array}\right.\right.\right]\hfill \\ \hfill & \left.\left.\times {\varphi }_F{\mathbf{G}}_{3,7}^{6,1}\left[{\psi }_F\lambda \left|\begin{array}{c}1;{\displaystyle \frac{g_F^2+1}{2}},{\displaystyle \frac{g_F^2+2}{2}}\\ {\displaystyle \frac{g_F^2}{2}},{\displaystyle \frac{g_F^2+1}{2}},{\displaystyle \frac{{\alpha }_F}{2}},{\displaystyle \frac{{\alpha }_F+1}{2}},{\displaystyle \frac{{\beta }_F}{2}},{\displaystyle \frac{{\beta }_F+1}{2}};0\end{array}\right.\right.\right]\right)\hfill \\ \hfill & \times H_2\left\{\left[{\displaystyle \frac{D_2^d}{V}}\left(J_{2,1}^S-J_{2,2}^S\right)+{\displaystyle \frac{D_2^d}{W}}\left(J_{2,3}^S-J_{2,2}^S\right)\right]\right.\hfill \\ \hfill & -e^{-{\lambda }_{B{R}_n}{E}_3}\left\{{\displaystyle \frac{E^d}{V}}\left[{\displaystyle \frac{1}{Q_{2,1}}}\left(e^{-{\displaystyle \frac{Q_{2,1}{D}_2}{E_4}}}+e^{-Q_{2,1}{F}_4}\right)\right.\right.\hfill \\ \hfill & \left.-{\displaystyle \frac{1}{Q_{2,2}}}\left(e^{-{\displaystyle \frac{Q_{2,2}{D}_2}{E_4}}}+e^{-Q_{2,2}{F}_4}\right)\right]\hfill \\ \hfill & +{\displaystyle \frac{E^d}{W}}\left[{\displaystyle \frac{1}{Q_{2,3}}}\left(e^{-{\displaystyle \frac{Q_{2,3}{D}_2}{E_4}}}+e^{-Q_{2,3}{F}_4}\right)\right.\hfill \\ \hfill & \left.\left.\left.-{\displaystyle \frac{1}{Q_{2,2}}}\left(e^{-{\displaystyle \frac{Q_{2,2}{D}_2}{E_4}}}+e^{-Q_{2,2}{F}_4}\right)\right]\right\}\right\}.\hfill \end{array}$$

(66)

Next, we derive the OP for the two users in the third working scenario, which is given by (36):

$$\begin{array}{cc}\hfill P_k^l& =1-\mathrm{Pr}\left({\gamma }_{SR}\ge {\theta }_{k,1}^{*}\right)\mathrm{Pr}\left(\right|h_{R{U}_k}{|}^2\ge {\theta }_{k,2})\hfill \\ \hfill & =1-\left[1-F_{{\gamma }_{SR}}\left({\theta }_{k,1}^{*}\right)\right]\left[1-F_{|h_{R{U}_k}{|}^2}\left({\theta }_{k,2}\right)\right],\hfill \end{array}$$

(67)

where ${\theta }_{1,2}=N_0{\gamma }_{th}^{1,2}/P_R\left(a_1^2-{\gamma }_{th}^{1,2}{a}_2^2\right)$, ${\theta }_{2,2}=\max \left\{N_0{\gamma }_{th}^{2,2}/P_R\left(a_1^2-{\gamma }_{th}^{2,2}{a}_2^2\right),N_0{\gamma }_{th}^{2,2}/P_R\left(a_2^2\right)\right\}.$

Bringing (18) and (20) into (36) gives:

$$\begin{array}{cc}\hfill P_1& =1-\left(1-{\varphi }_T{\mathrm{G}}_{2,3}^{2,1}\left[{\psi }_T{\lambda }^{{\displaystyle \frac{{\alpha }_T}{2}}}\left|\begin{array}{c}1;1+{\displaystyle \frac{g_T^2}{{\alpha }_T}}\\ {\displaystyle \frac{g_T^2}{{\alpha }_T}},{\mu }_T;0\end{array}\right.\right.\right]\hfill \\ \hfill & \left.\left.\times {\varphi }_F{\mathbf{G}}_{3,7}^{6,1}\left[{\psi }_F\lambda \left|\begin{array}{c}1;{\displaystyle \frac{g_F^2+1}{2}},{\displaystyle \frac{g_F^2+2}{2}}\\ {\displaystyle \frac{g_F^2}{2}},{\displaystyle \frac{g_F^2+1}{2}},{\displaystyle \frac{{\alpha }_F}{2}},{\displaystyle \frac{{\alpha }_F+1}{2}},{\displaystyle \frac{{\beta }_F}{2}},{\displaystyle \frac{{\beta }_F+1}{2}};0\end{array}\right.\right.\right]\right)\hfill \\ \hfill & \times {\left(1-e^{-{\lambda }_{R{U}_A}{\theta }_{1,2}}\right)}^N,\hfill \end{array}$$

(68)

and

$$\begin{array}{cc}\hfill P_2& =1-\left(1-{\varphi }_T{\mathrm{G}}_{2,3}^{2,1}\left[{\psi }_T{\lambda }^{{\displaystyle \frac{{\alpha }_T}{2}}}\left|\begin{array}{c}1;1+{\displaystyle \frac{g_T^2}{{\alpha }_T}}\\ {\displaystyle \frac{g_T^2}{{\alpha }_T}},{\mu }_T;0\end{array}\right.\right.\right]\hfill \\ \hfill & \left.\left.\times {\varphi }_F{\mathbf{G}}_{3,7}^{6,1}\left[{\psi }_F\lambda \left|\begin{array}{c}1;{\displaystyle \frac{g_F^2+1}{2}},{\displaystyle \frac{g_F^2+2}{2}}\\ {\displaystyle \frac{g_F^2}{2}},{\displaystyle \frac{g_F^2+1}{2}},{\displaystyle \frac{{\alpha }_F}{2}},{\displaystyle \frac{{\alpha }_F+1}{2}},{\displaystyle \frac{{\beta }_F}{2}},{\displaystyle \frac{{\beta }_F+1}{2}};0\end{array}\right.\right.\right]\right)\hfill \\ \hfill & \times {\lambda }_{R_n{U}_1}{\lambda }_{R_n{U}_2}N\sum _{i=0}^{N-1}\left(\begin{array}{c}N-1\\ i\end{array}\right){(-1)}^{N-i-1}\hfill \\ \hfill & \times \left\{{\displaystyle \frac{1}{V}}\left[{\displaystyle \frac{e^{-\left(N-i\right){\lambda }_{R{U}_A}{\theta }_{2,2}}}{\left(N-i\right){\lambda }_{R{U}_A}}}-{\displaystyle \frac{e^{-{\lambda }_{R_n{U}_1}{\theta }_{2,2}}}{{\lambda }_{R_n{U}_1}}}\right]\right.\hfill \\ \hfill & \left.\left.+{\displaystyle \frac{1}{W}}[{\displaystyle \frac{e^{-\left(N-i\right){\lambda }_{R{U}_A{\theta }_{2,2}}}}{\left(N-i\right){\lambda }_{R{U}_A}}}-{\displaystyle \frac{e^{-\lambda R_n{U}_2{\theta }_{2,2}}}{{\lambda }_{R_n{U}_2}}}\right]\right\}.\hfill \end{array}$$

(69)

4. Numerical Results

In this section, we present the numerical results and discuss the system. We present

Table 1. Parameter settings for FSO, THz, and RF links without special instructions.

	Parameter	Value
FSO Link	Wavelength	${\lambda }_F=1550\mathrm{nm}$
	Link length	$L_F=200\mathrm{m}$
	Visibility	$V=20\mathrm{km}$
	Weak turbulence	$C_n^2=1\times {10}^{-13}{\mathrm{m}}^{-2/3}$
	Moderate turbulence	$C_n^2=1\times {10}^{-12}{\mathrm{m}}^{-2/3}$
	Small pointing errors	${\mathcal{A}}_F=0.2\mathrm{m}$, ${\omega }_{z,F}=0.4\mathrm{m}$,
		${\sigma }_{s,F}=0.06\mathrm{m}$, $g_F=3.8122$
	Large pointing errors	${\mathcal{A}}_F=0.2\mathrm{m}$, ${\omega }_{z,F}=0.4\mathrm{m}$,
		${\sigma }_{s,F}=0.14\mathrm{m}$, $g_F=1.6338$
THz Link	Frequency	$f_T=100\mathrm{GHz}$
	Link length	$d_T=200\mathrm{m}$
	Antenna gain	$G_t^T=G_r^T=55\mathrm{dB}$
	Weak turbulence	${\alpha }_T=1.7$, ${\mu }_T=1.7$
	Small turbulence	${\alpha }_T=1.7$, ${\mu }_T=1.5$
RF Link	Power Allocation Coefficient (PAC)	$a_1=0.35$, $a_2=0.33$, $a_3=0.32$
	Saturation threshold	${\gamma }_{th}^{1.1}={\gamma }_{th}^{2.1}={\gamma }_{th}^{1.2}={\gamma }_{th}^{2.2}=-2\mathrm{dB}$

, which includes typical data for the system used, including performance parameters in the FSO, THz, and RF links. Specifically, regarding the data presented: ${\mathcal{A}}_F$ is the receiver radius, ${\omega }_{z,F}$ is the equivalent beam width, and ${\sigma }_{s,F}$ is the pointing error displacement standard deviation at R. If new conditions need to be introduced in the qualitative analysis, we will explain them.

The beamwidth and jitter standard deviation play a crucial role in antenna transmission. To clearly observe the impact of the beamwidth and jitter standard deviation on the OPs of two users under three different scenarios, Figure 2 illustrates the comparison between different beamwidth and jitter standard deviations versus OPs. The first row shows the OP of $U_1$ under three operating scenarios, and the second row shows the OP of $U_2$ under the same scenarios. From Figure 2, we can observe that the OP performance of both users decreases with an increase in beamwidth or an increase in jitter standard deviation. However, since the final equivalent beamwidth is determined by the beamwidth along with the receiving radius, the OP also decreases slightly in the range where the jitter standard deviation equals the receiving radius. Furthermore, we observe that in the first two operating scenarios, the OP is significantly better than in the third operating scenario due to the energy received from the PB. In the same operating scenario, the OP performance of $U_1$ is better than that of $U_2$.

Ensuring high performance while transmitting messages over longer distances is a fundamental requirement for every communication system. To investigate the performance of the proposed system at different distances, Figure 3 illustrates the impact of THz frequency and link transmission distance on the OPs of two users in three scenarios. The first row displays the OP of $U_1$ in the three operating scenarios, while the second row shows the OP of $U_2$ in the same scenarios. The findings indicate that a shorter THz link transmission distance and a lower THz frequency result in improved performance for both users in the system. For the same reason as in Figure 2, both the first and second operating scenarios are significantly better than the third operating scenario.

In practical communication transmission, weather conditions often have a significant impact. Therefore, Figure 4 illustrates the impact of two different visibility levels on the OPs for two users across three scenarios. The significance of its practical application is that the small V scene simulates the communication condition of low visibility on rainy days, and the large V scene simulates the communication condition of high visibility on sunny days. From Figure 4, it can be seen that the scenario with better visibility exhibits superior OP performance for a specific SR link SNR condition. However, as the SR link SNR increases, the entire system becomes more limited by the second-hop RU link, which ultimately stabilizes at a fixed value, namely the outage probability of the $RU$ link.

Figure 5 illustrates the impact of atmospheric turbulence and pointing errors on the OPs of two users in three scenarios. We split the various atmospheric turbulence conditions and pointing errors into two cases for comparison. In the first case, the THz link is parameterized by ${\alpha }_T=1.7$, ${\mu }_T=1.5$, and $g_T=2.0859$, while the FSO link is parameterized by $C_n^2=1\ast {10}^{-12}$ and $g_T=1.6338$. In contrast, the THz link in the second case is parameterized by ${\alpha }_T=1.7,{\mu }_T=1.7$, and $g_T=2.0859$, while the FSO link is parameterized by $C_n^2=1\ast {10}^{-13}$ and $g_T=1.6338$. The first case represents the actual scenario of strong turbulence with large pointing errors, while the second case represents the scenario of low turbulence with small pointing errors. As can be seen from Figure 5, the link performance is better in the scenario with low turbulence. However, once the $SR$ link SNR is good enough, the overall system performance for both cases is limited by the $RU$ link performance and stabilizes at the same value.

The number of transmitting antennas also affects the performance of the communication system. To explore the impact of OP versus N among three working scenarios, we obtain Figure 6. From Figure 6, it can be seen that the OP performance of the system improves as the number of antennas increases. However, this trend of improvement is not infinite. When the number of antennas increases to a certain level ($N=7$), the overall performance starts to be affected by the $SR$ link and the $RU$ link. Without changing the other parameters, the enhancement in OP performance gradually decreases, eventually leading to the OP for both users converging to a stable value.

Different power gain results in varying OP values; in order to explore the impact of the relationship between OP performance and the channel coefficient $g_{u1}$ for two users in three scenarios, we obtain Figure 7. Based on the trends in Figure 7, we categorize the three working scenarios into two groups. The first group includes the first and second scenarios, both utilizing the EH system and therefore exhibiting similar curve shapes, whereas the second group pertains to the third working scenario. In both categories, when $g_{u1}$ is very small, the OP of $U_1$ is close to 1, i.e., the link will almost certainly fail to communicate. As $g_{u1}$ increases, the OP performance of $U_1$ improves earlier and more significantly because $g_{u1}$ is the channel coefficient of the $R{U}_1$ link. Due to the NOMA system, the $R{U}_1$ link also transmits a portion of $U_2$ information. As $g_{u1}$ increases, $U_2$’s OP performance also improves. However, in general, and particularly in the first type of scenarios, it converges to a stable value earlier. All six curves in Figure 7 eventually converge to a stabilized value. This convergence occurs because the overall system’s OP performance is constrained by the various conditions of the first-hop SR link when $g_{u1}$ is sufficiently boosted.

In practical application scenarios, two users need to allocate a portion of the transmission power for sending messages. This allocation ratio will impact their respective OP. Figure 8 shows the relationship between ${\alpha }_1$ and the OP performance of the three scheme systems. As the value of ${\alpha }_1$ gradually increases, more power will be allocated to $U_1$, which will improve the OP performance of $U_1$, while the OP performance of $U_2$ gradually deteriorates. However, it can be seen that the OP performance of $U_1$ improves relatively slowly. This is because simply increasing the power allocation factor without changing the total transmit power has a very limited impact on the system as a whole. Therefore, the OP of $U_1$ is more constrained by the rest of the system. Conversely, the OP performance of $U_2$ deteriorates drastically after ${\alpha }_1$ exceeds 0.9, and the OP eventually approaches 1. This is because the allocated power is so small that it is almost impossible to achieve normal communication transmission.

In order to explore the relationship between the collected energy and the OP of the system more clearly, we decide to plot the first working scenario, which is not constrained by CR, and set three different $P_{self}$ values for distinction. Figure 9 illustrates the relationship between $\eta /P_B$ and the OP performance of the two users in the first working scenario. As can be seen from Figure 9, initially, increasing the saturation threshold improves the system’s performance. This is because the energy used for transmission also increases. After the saturation threshold is increased to a certain level (approximately 10 dB), the OP performance for both users no longer changes. This is because the transmit power of the PB does not change, and increasing the saturation threshold does not increase the energy received by the relay after the saturation threshold exceeds the transmit power. In addition, the three different $P_{self}$ values are 1.5 dB for case 1, 2 dB for case 2, and 2.5 dB for case 3. It can be seen that the different $P_{self}$ values affect the initial system performance when the saturation threshold is very low (−20 dB), but these differences become progressively smaller with the increase in received energy. Moreover, the entire system is influenced by the first hop, and eventually, all three different cases stabilize at a similar value.

Figure 10 illustrates the relationship between I and the OP performance of the two users in the third working scenario. The third working scenario is chosen because there is no influence of the EH system in this scenario. It is a very classical CR model that allows for a clearer exploration of the impact of I on the transmission network. Since the third working scenario is energy-rich or has a stabilized power supply, we distinguish here three cases with higher $P_{self}$ values, namely $P_{self}$ = 5 dB, $P_{self}$ = 10 dB, and $P_{self}$ = 15 dB. As can be seen in Figure 10, as I is increased, i.e., the PN is able to allow the system to transmit higher power, the overall performance of the system is improved. However, this improvement is not infinite. When I exceeds $P_{self}$, the actual power sent by R will no longer change, and the OP performance of the two users will remain constant.

5. Conclusions

In this paper, we investigate a hybrid-assisted, dual-hop link system for DF relaying based on the deployment of the NOMA protocol, incorporating EH and CR. The first hop of the system consists of a hybrid FSO/THz link switched by a hard switching. In this setup, the FSO link follows Gamma–Gamma distribution and the THz link follows $\alpha -\mu$ distribution. The second-hop link utilizes an RF link following a Rayleigh fading distribution. The channel for EH energy transmission also exhibits the characteristics of Rayleigh fading distribution. These statistical distributions are used to derive exact expressions for the PDF and CDF of the channel power coefficients for each link. By incorporating three different real-world operating scenarios, the exact closed expressions for the OPs of two users under the three operating scenarios are further derived. In addition, the study investigates the impact of various system parameters on the OP performance of the system. These parameters include atmospheric turbulence, pointing error, visual range, number of antennas, power allocation coefficient, EH saturation threshold, and channel power gain coefficient. The numerical results confirm the superiority of the hybrid system and its suitability for various operating scenarios. In the process of research, we mainly face the following challenges: First, due to resource limitations, the model proposed in this article is difficult to validate in practical environments and cannot achieve comprehensive experimental verification. Secondly, the adoption of some assumptions simplifies the complex problem, but there may be some differences from the real scenario. In addition, computational complexity is an issue that needs further optimization. To sum up, a new scheme is proposed in this paper to address the issues of spectrum resource scarcity and communication equipment energy life in communication systems simultaneously. This model offers a valuable approach for future research.