A PHY Layer Security Analysis of a Hybrid High Throughput Satellite with an Optical Feeder Link
Elmehdi Illi Faissal El Bouanani Fouad Ayoub Mohamed-Slim Alouini
cs.IT math.IT







Hybrid terrestrial-satellite (HTS) communication systems have gained a
tremendous amount of interest recently due to the high demand for global
high data rates. Conventional satellite communications operate in the
conventional Ku (12 GHz) and Ka (26.5-40 GHz) radio-frequency bands for
assessing the feeder link, between the ground gateway and the satellite.
Nevertheless, with the aim to provide hundreds of Mbps of throughput per
each user, free-space optical (FSO) feeder links have been proposed to
fulfill these high data rates requirements. In this paper, we investigate
the physical layer security performance for a hybrid very high throughput
satellite communication system with an FSO feeder link. In particular, the
satellite receives the incoming optical wave from an appropriate optical
ground station, carrying the data symbols of N users through various
optical apertures and combines them using the selection combining technique.
Henceforth, the decoded and regenerated information signals of the N users
are zero-forcing (ZF) precoded in order to cancel the interbeam
interference at the end-users. The communication is performed under the
presence of malicious eavesdroppers nodes at both hops. Statistical
properties of the signal-to-noise ratio of the legitimate and wiretap links
at each hop are derived, based on which the intercept probability metric is
evaluated. The derived results show that above a certain number of optical
apertures, the secrecy level is not improved further. Also, the system's
secrecy is improved using ZF precoding compared to the no-precoding scenario
for some specific nodes' positions. All the derived analytical expressions
are validated through Monte Carlo simulations.


A PHY Layer Security Analysis of a Hybrid High Throughput Satellite
with an Optical Feeder Link
    Elmehdi Illi, Member, IEEE, Faissal El Bouanani, 
Senior Member, IEEE, Fouad Ayoub,
Member,
IEEE, 


and Mohamed-Slim Alouini, Fellow, IEEE 
E. Illi and F. El Bouanani are with ENSIAS College of Engineering, Mohammed
V University, Rabat, Morocco (e-mails: {elmehdi.illi,
f.elbouanani}@um5s.net.ma). 
F. Ayoub is with CRMEF, Kenitra, Morocco (e-mail: ayoub@crmefk.ma). 
M.-S. Alouini is with Computer, Electrical, and Mathematical Sciences and
Engineering (CEMSE) Division, King Abdullah University of Science and
Technology (KAUST), Thuwal 23955-6900, Makkah Province, Saudi Arabia
(e-mail: slim.alouini@kaust.edu.sa).
    Final Version, December 30, 2023
==============================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================



§ INTRODUCTION


Throughout the last few years, satellite communication
(SatCom) has been a tremendously evolving segment of the wireless
communication industry, due to the increasing global demand on broadband
satellite communication links <cit.>. Interestingly, with the
arrival of the fifth-generation (5G) wireless cellular network, a variety of
satellite operators on the globe are developing broadband communications to
complement and compete with the terrestrial cellular networks <cit.>.
In this regard, multibeam SatCom has been widely advocated as an appropriate
way to assess very high-speed satellite-users link, in which a large number
of spot beams is used. Satellite links aim at providing high-speed
communications (order of hundreds of Mbps per user) to users in areas where
the traditional terrestrial networks offer very low quality of service <cit.>.

Among the critical challenges faced by the satellite communication industry
is the spectrum scarcity issue <cit.>, <cit.>. For
instance, with the high bandwidth requirement of the end-users, neither the
Ku-band (12 GHz) nor the Ka-band (26.5-40 GHz) seems to fulfill the hundreds
of Gbps aggregate user link throughput, due to the scarce frequency
resources in these bands <cit.>. In addition to this, multiple ground
stations (gateways) are needed to feed the satellite to assess the desired
data rates while operating on such bands, which results in significant
energy consumption <cit.>.

To this end, free-space optics (FSO) technology has been broadly endorsed as
an effective solution for providing very high data rate links on
terrestrial-satellite communications <cit.>. The overarching idea is to
carry on data from the optical ground station (OGS) in the form of conical
light beams using a powered laser device, operating either on the visible
(400-800 nm) or infrared (1500-1550 nm) spectrum, to a satellite, which
converts the optical signal to an electrical one, and serves the end-users
through radio-frequency (RF) spot beams <cit.>. From another
front, the use of optical bands does not require any regulation or license
fees as done in traditional Ku and Ka bands, where the International
Telecommunications Union (ITU) regulates their use. Also, besides its
immunity to interference and high security, FSO communication can provide a
data rate in the order of Tbps per optical beam, which renders it a viable
alternative solution to reach the desired high data rates <cit.>. On
the other hand, optical communication is highly affected by atmospheric and
weather losses along the propagation path, known as turbulence. Pointing
error due to transmitter/receiver misalignment, as well as free-space
path-loss are two other limiting factors of FSO in outdoor communications
<cit.>. Furthermore, another restricting phenomenon in
ground-to-satellite FSO links is cloud coverage. Indeed, optical
communication between the gateway and the satellite is blocked totally in
the presence of clouds <cit.>.

During the last decade, several works and researches have been conducted
onto the deployments of hybrid terrestrial-satellite (HTS) systems with an
optical feeder link. In <cit.>, the performance analysis of
optical feeder link with gateway diversity, in the presence of stochastic
cloud coverage model, is assessed. Furthermore, the authors in <cit.>
carried out a performance analysis of an HTS system with an optical feeder
link, with zero-forcing (ZF) precoding to remove the interbeam interference
(IBI), as well as amplification at the satellite are considered. In addition
to this, from an industrial perspective, several broadband optical
feeder-based HTS systems have been demonstrated in the last few years. For
instance, the first successful ground-to-satellite optical link was
performed between the earth and ETS-VI satellite in Konegi, Japan <cit.>. Besides this, other HTS systems experiments demonstrated the
achievability of great throughput records, such as NASA's 622 Mbps Laser
communication experiments in 2014 <cit.>, and DLR's Institute of
Communications and Navigation HTS experiments with a record data rate of
1.72 Tbit/s in 2016, and 13.16 Tbit/s in 2017 <cit.>. NICT is
planning in 2021 to launch a new test satellite with an aim to demonstrate a
10 Gbps speed on uplink and downlink of aggregate throughputs <cit.>
. Moreover, in <cit.>, an overview of implementation,
performance, technological aspects, and users' quality of service on HTS
systems with an optical feeder is assessed.

From another front, privacy and security are becoming a big concern in such
networks where considerable attention from the research community has been
paid. Importantly, the broadcast nature of the wireless RF link renders it
vulnerable to eavesdropping attacks <cit.>. While higher layers view
the security aspect as an implementation of cryptographic protocols, the
physical layer (PHY) security, introduced by Wyner, aims at establishing
secure transmissions, by ensuring that the data rate of the legitimate link
exceeds that of the wiretap one by a certain threshold.

From a multibeam satellite communication point of view, the legitimate
users-links are established through narrow multiple spot beams, where each
beam targets a single cell on the covered zone, rendering the communication
most unlikely to being intercepted from distant malicious nodes.
Nevertheless, potential eavesdroppers might be located in the same zone as
the legitimate users, and consequently, the secrecy level of the user link
is affected. To this end, several works in the literature have dealt with
the secrecy level of HTS systems as in <cit.>,
where the analysis carried out the performance of HTS relay-based networks,
where the feeder link operates on RF spectrum. On the other hand, several
works such as <cit.> assessed the PHY security of FSOlinks.



 §.§ Motivation and Contributions


Very few works in the literature have investigated the secrecy level of HTSsystems with optical feeder link. In particular, the authors in <cit.> dealt with the secrecy analysis of an HTS relay network, where
the optical link is used at the terrestrial side as a last-mile link. Also,
the wiretapper is considered only on the RF side. Distinctively,
contributions such as <cit.> dealt with the
secrecy level of mixed RF-FSO links where the RF link is considered in the
first hop. Furthermore, nodes parameters, a data precoding process, and an
eavesdropper on the optical link were not considered in such works. Besides,
differently from the works <cit.> where the authors proposed a
transmit Laser selection diversity techniques for FSO systems, the
statistics of the combined SNR using receive diversity selection combining
(SC) scheme were not carried out in closed-form expression. Capitalizing on
this, we aim at this work to investigate the PHY layer security of an HTSmulti-user relay-based system with an optical feeder link, where the
satellite, acting as a relay, converts the incoming optical wave carrying
the users' data, after combining the incoming beams to its photodetectors
through SC scheme, to an electrical signal, regenerate and conveys them to
the end-users. Two scenarios are analyzed, namely, i) the satellite
performs ZF precoding technique, after decoding information signals on the
second hop, before transmitting it to the end-users, ii) the satellite does
not perform ZF technique and delivers the processed received signal to the
end-users. Malicious eavesdroppers are considered on each hop of the
transmission. The main contributions of this paper can be summarized as
follows:



  * Statistical properties of the end-to-end secrecy capacity are
retrieved, in terms of legitimate and wiretap instantaneous SNRs of the S-R and R-D hops.


  * A novel expression for the IP of dual-hop DF relaying-based systems is
retrieved, where a decoding failure event at the satellite is considered.


  * Capitalizing on the above two results, a closed-form expression of the
intercept probability (IP) of the system is derived for the ZF and non-ZFscenarios.


  * An asymptotic analysis of the derived analytical result is performed
based on which the achievable coding gain and diversity order are quantified.




 §.§ Organization of the Paper


The remainder of this paper is organized as follows. Section II is
dedicated to present the system and channel model. Section III deals with
the statistical properties of the end-to-end SNRs, while Section IV depicts
the derived analytical expression of the IP. In Section V, illustrative
numerical results are shown to assess the effect of channel parameters on
the system's secrecy level. Finally, Section VI concludes the paper.



§ SYSTEM AND CHANNEL MODEL


We consider in this analysis a single OGS communicating with a GEOsatellite through an optical feeder link. The transmit OGS, acting a source (
S) node, is assumed to have a clear line of sight (LOS) link with the
satellite unblocked by clouds. The selected source node transmits data
symbols of the N users through a turbulent optical channel to the
satellite. This latter, acting as a relay (R), combines through K
optical photodetectors the incoming optical wave, converts it to the
electrical domain, performs SC technique, decodes, and delivers it in the
form of N beams to N user earth stations ( T_i) _1≤
i≤ N. If the received SNR at the satellite is greater than a
predefined decoding threshold, it can successfully decode the information
signal. Otherwise, the decoding cannot be ensured correctly. For the former
case, the satellite can either forward the regenerated signals directly or
precode them using ZF technique. We assume that an eavesdropper attempting
to overhear the divergent optical beam coming from the OGS to the satellite
(FSOS-R hop). In addition to this, another potential wiretapper
attempts to overhear the signal carried to the earth-station T_i. 





 §.§ S-R Link




  §.§.§ Legitimate Link


An OGS with a LOS unblocked by clouds is considered for data transmission.
Without loss of generality, all users' signals are assumed to have the same
power.

The received beam at the satellite contains N multiplexed signals
transmitted from the transmit gateway. Hence, the received electrical signal
vector at the satellite's k-th aperture, after demultiplexing, is
expressed as <cit.>

    𝐲_1^(k)=√(ω _lP_S)( η I_k) ^
    r/2𝐮+𝐧_k,k=1,… ,K,

with I_k=I_k^(a)I_k^(p)I_k^(ℓ ), being the product of the
irradiance fluctuation due to atmospheric turbulence, the pointing error due
to the beam misalignment, and the free-space path loss, respectively, with I_k^(ℓ )=I_te^-ϕ d, and I_t, d, and ϕ denote the
laser emittance, the station-satellite distance, and the path loss exponent,
respectively. r∈{1,2} being a detection-technique dependent
parameter, with r=1 referring to coherent detection, and r=2 stands for
direct detection, and ω _l denotes the portion of power received by
the satellite's photodetectors, among the total radiated power from the
selected OGS. Also:



  * 𝐮=[u_1,u_2,… ,u_N]^T denotes the transmitted
signal vector of N signals, with the superscript ^T refering to the
transpose of a vector. Each unit-power signal is modulated onto one of the
N optical sub-carriers, where 𝔼[ 𝐮^H𝐮
] =N, with 𝔼[ .] denoting the
expectation operator.


  * P_S is the OGS transmit power.


  * 𝐧_k=[n_k^(1),n_k^(2),… ,n_k^(N)]^T
stands for the additive white Gaussian noise (AWGN) process at the satellite
with zero mean and the same variance σ _1^2.


The satellite converts the received optical waves at the K apertures into
electrical signals, and then uses SC to choose the branch with highest
instantaneous SNR as: s^∗=1≤ k≤ Kmaxγ
_1^(k), with γ _1^(k)=P_Sω _l( η
I_k) ^r/σ _1^2 denoting the instantaneous
signal-to-noise ratio (SNR) received at the k-th satellite's aperture.
Consequently, the combined SNR at the satellite is: γ _1=
1≤ k≤ Kmaxγ _1^(k). Afterwards, the satellite performs
a decoding process on the combined electrical signal to regenerate the
information signal again.



  §.§.§ Wiretap Link


The hybrid ground-satellite communication is performed under the malicious
attempt of eavesdroppers per each one of the two hops to intercept the
legitimate message. For the FSO link, we consider the presence of one
wiretapper e^(1) located at some altitude within the divergence region
of the OGS beam, being able to capture a portion ω _e=1-ω _l
of the optical power. That is, the received SNR at e^(1) is γ
_1^(e)=ω _eP_S( η I_e^(1)) ^r/σ
_e^(1)^2, with I_e^(1) and σ _e^(1)^2 denote the
respective optical channel gain, and the variance of the AWGN at e^(1),
respectively.

In terrestrial FSO communication, the altitude-dependent refractive index
structure parameter C_n,ϖ^2(h) ( ϖ∈{1,e}) as well as the Rytov variance σ _R,ϖ^2 are two crucial
parameters that represent the atmospheric turbulence and pointing error loss
impairments. Such parameters are expressed in key system and environment
quantities. In the context of vertical optical links in HTS systems, the
altitude-dependent refractive structure index parameter in m^-2/3
and the Rytov variance can be expressed for the uplink using the
Hufnagel-Valley Boundary model as given in <cit.>, in
terms of the altitude h in meters, the wind speed V_w in m/s, the
satellite and OGSaltitudes d_ϖ and h_0, the satellite's
zenith angle with respect to the OGS, and the operating wavelength λ
.

Interestingly, ξ _ϖ^2=W_eq/2σ _s^2 denotes
the pointing error strength, which is the ratio between the equivalent beam
waist at the satellite altitude and the beam wander displacement variance,
where the last-mentioned quantities are defined in <cit.>, <cit.>, and <cit.>.



 §.§ R-D Link


The satellite generates N adjacent beams. The received signal vector at
the N legitimate end-users and wiretap nodes can be formulated as

    𝐲_2^(ϰ )=𝐇^(ϰ )𝐱+𝐧
    _ϰ,

with ϰ equals either l for the legitimate earth stations or e^(2) for the second hop's wiretappers, P_SAT is the satellite
transmit power, 𝐇^(ϰ ) is the channel matrix between
the N satellite antennas and the nodes [ ϰ _i]
_i=1,..,N[
The subscripts/superscripts "l" and "e^(2)" are used to denote the
legitimate and wiretap links at the second hop, respectively.], and 𝐧_ϰ=[ n_ϰ^(i)] _i=1,..,N is the AWGNvector whose elements are zero mean and with the same variance σ
_ϰ^2.

It is known that the channel matrix 𝐇^(ϰ ) can be
decomposed as <cit.>

    𝐇^(ϰ )=𝐃^(ϰ )𝐕^(ϰ ),

with 𝐃^(ϰ ) being a diagonal matrix containing
real-valued random fading coefficients, and 𝐕^(ϰ )
entries are path-loss and radiation pattern coefficients, defined by V_i,j^(ϰ )=c√(G_T^(j)G_R^(ϰ_i)
a_j^(ϰ _i))/4π fr^(ϰ _i)√(κ TB_W)
<cit.>, where c being the light celerity in the free space, G_T^(j) and G_R^(ϰ _i) are the respective gains of the
satellite's j-th transmit antenna and the receive antenna of the node ϰ _i, κ is the Boltzmann constant, T is the receiver
noise temperature, f is the operating frequency, r^(ϰ _i) is
the distance between the satellite and the node ϰ _i, and a_j^(ϰ _i) is the normalized beam radiation gain, which can
be approximated as a_j^(ϰ _i)≈( J_1(
ϑ _j^(ϰ _i)) /2ϑ _j^(ϰ _i)
+36J_3( ϑ _j^(ϰ _i)) /(
ϑ _j^(ϰ _i)) ^3) ^2 <cit.>,
with ϑ _j^(x_i)=2.07123sin( ϕ
_j^(ϰ _i)) /ϕ _j,3dB, and J_n(.)
denotes the n-th order Bessel function of the first kind <cit.>, and ϕ _j^(ϰ _i) denote the angle between the
node ϰ _i and the j-th beam boresight as indicated in (<ref>), given as

    ϕ _j^(ϰ _i)={[ ϕ( i-j) +ϕ _j^(ϰ _j),j≥ i;  ϕ( i-j) -ϕ _j^(ϰ _j),j<i ]
    . ,

with ϕ = 𝒟/r referring to the beam width angle,
with 𝒟 being the cells' diameter assumed to be equal for all
cells, and r is the distance between each node ϰ _i and the
satellite, assumed to be equal for all nodes (i.e., r^(ϰ
_i)≈ r,∀ i≤ N). Also, ϕ _j,3dB denotes the
angle corresponding to 3dB power loss of the j-th beam.



  §.§.§ With ZF Precoding




  
Legitimate Link

We consider the case of a ZF receiver employed in <cit.>, where the
transmit signals are precoded, after decoding, at the satellite before
transmitting them to the end-users. Thus, the transmit vector 𝐱
is expressed as: 𝐱=𝐌𝐮, where 𝐌=√(φ)
𝐀 is a N× N precoding matrix with 𝐀=
( 𝐕^(l)) ^-1, and

    φ =P_SAT/N[ ( 𝐕^(l)( 𝐕
    ^(l)) ^H) ^-1] ,

where the superscript ^H refers to the Hermitian operator defined as the
conjugate of the transpose matrix, and [.] stands for the matrix's
trace.

After performing ZF precoding at the satellite, the received signal and
SNR at the i-th earth-station, without considering the first hop, can be
formulated as follows

    y_ZF^(l_i)   =   √(φ)D_i^(l)u_i+n^(l_i), 
    γ _ZF^(l_i)   =   γ_ZF^(l_i)|
    D_i^(l)| ^2,

with D_i^(ϰ ) denoting the i-th element of the diagonal of 𝐃^(ϰ ), and γ_ZF^(l_i)=
φ/σ _l^2. Since the satellite performs DF protocol, the
equivalent received SNR at the i-th earth-station is given as

    γ _eq,ZF^(l_i)=min( γ _1,γ
    _ZF^(l_i)) .




  
Wiretap Link

As far as the second hop is concerned, the satellite-stations links are
ensured under the potential presence of oneeavesdropper (
e_i^(2)) _1≤ i≤ N per each cell i, aiming to intercept
the RF beam transmitted to it. The received signal and instantaneous SNRat e_i^(2) can be expressed, respectively, as follows

    y_e_i^(2),ZF    =√(φ)D_i^( e^(2))
    ∑_l=1^Nu_l∑_j=1^NV_i,j^( e^(2))
    A_j,l+n^(e_i^(2)), 
    γ _e_i^(2),ZF    =ψ _iγ _ZF^(
    e_i^(2)) /θ _iγ _ZF^( e_i^(2)) +1
    ,

where n^(e_i^(2)) denotes an AWGN process at e_i^(2) with zero
mean and the same variance for all wiretap nodes σ _e^(2)^2, γ _ZF^( e_i^(2)) =γ
_ZF^(e_i)| D_i^( e^(2)) | ^2, γ_ZF^( e_i^(2)) = φ/
σ _e^(2)^2, ψ _i=(
∑_j=1^NV_i,j^(e^(2))A_j,i) ^2, and θ
_i=( ∑_m=1,m≠
i^N∑_j=1^NV_i,j^(e^(2))A_j,m) ^2.



  §.§.§ Without ZF Precoding


When no ZF precoding is performed at the transmit OGS, the received signal
and SNR at the legitimate earth station l_i as well as the eavesdropper e_i^(2) are given as [
The subscript "NZF" is used to denote the "non-ZF case", while the "ZF"
subscript denotes the ZF adoption scenario.]

    y_NZF^(ϰ _i)   =   √(P_SAT/N)D_i^(ϰ
    )∑_j=1^NV_i,j^(ϰ )u_j+n^(ϰ _i), 
    γ _ϰ _i,NZF   =   Ψ ^(ϰ _i)γ
    _NZF^(ϰ _i)/Θ ^(ϰ _i)γ
    _NZF^(ϰ _i)+1,

with γ _NZF^(ϰ _i)=γ_NZF^(ϰ
_i)| D_i^(ϰ )| ^2,γ
_NZF^(ϰ _i)=P_SAT/Nσ _ϰ^2, Ψ
^(ϰ _i)=( V_i,i^(ϰ )) ^2, and Θ
^(ϰ _i)=( ∑_j=1,j≠ i^NV_i,j^(ϰ
)) ^2.



§ STATISTICAL PROPERTIES


In this section, the cumulative distribution function (CDF) of the SNR of
the legitimate link both hops, as well as the wiretap link, is expressed in
terms of the system and channel parameters.



 §.§ S-R Link


The Gamma-Gamma distribution is considered for modeling the atmospheric
turbulence induced-fading with pointing error for each received optical
beam. The respective probability density function (PDF) and CDF of the
SNRs γ _1^(k) and γ _1^(e), received respectively at
the satellite's k-th aperture and the received one at E_1, are given
as <cit.>

    f_γ _1^(ϖ )(z)    =𝒫_ϖ/rz
    G_1,3^3,0( Υ _ϖ( z/μ _r^(ϖ )
    ) ^1/r|[         -;ξ _ϖ^2+1; ξ _ϖ^2,α _ϖ,β _ϖ;- ]
    . ) ,  
    
    F_γ _1^(ϖ )(z)    =( 2π) ^1-r𝒫
    _ϖ/r^2-α _ϖ-β _ϖG_r+1,3r+1^3r,1
    ( Υ _ϖ^rz/r^2rμ _r^(ϖ )|[ 1;κ _1^(ϖ ); κ _2^(ϖ );0 ]
    . ) ,

respectively, with Υ _ϖ=ξ _ϖ^2α
_ϖβ _ϖ/ξ _ϖ^2+1,ϖ∈{k,e}, 𝒫_ϖ=ξ _ϖ^2/Γ( α
_ϖ) Γ( β _ϖ) , μ
_r^(ϖ )=𝔼[ γ _1^(ϖ )], κ
_1^(ϖ )=( ξ _ϖ^2+i/r) _i=1,..,r, κ _2^(ϖ )=( ξ _ϖ^2+i/r,α
_ϖ+i/r,β _ϖ+i/r) _i=0,..,r-1, and G_p,q^m,n( .| .. ) refers to the Meijer's G
-function <cit.>. Moreover, the
turbulence-induced fading parameters α _ϖ and β _
 ϖ are expressed terms of the Rytov variance σ
_R,ϖ^2 given in <cit.> using <cit.>. Without loss of generality, we consider that the fading
amplitudes are i.i.d on all branches, that is α _k=α _1, β _k=β _1, μ _1^(k)=μ _1, ξ _k=ξ _1 for k=1,..,K.


The CDF of the combined SNR at the satellite can be expressed as follows

    F_γ _1(z)=𝒫_1^K∑_h_1+h_2+h_3=K∑
    _l=0^∞ℱ_h_1,h_2,h_3,l/(
    Υ z) ^-ϱ _l,h_1,h_2,h_3,

with ℱ_h_1,h_2,h_3,l=K!/h_1!h_2!h_3!(
a_0^(1))
^h_1∑_q_2+q_3=lc_q_2^(2)c_q_3^(3), a_0^(1)=Γ( α _1-ξ _1^2) Γ( β _1-ξ _1^2) /ξ _1^2, Υ =ξ
_1^2α _1β _1/μ _1( ξ _1^2+1) , ϱ _l,h_1,h_2,h_3=l+h_1ξ _1^2+h_2α
_1+h_3β _1,

    c_m^(i)={[                            ( a_0^(i)) ^h_i,m=0; 1/ma_0∑_j=1^m( jh_i-m+j)
    a_j^(i)c_m-j^(i);m≥ 1 ]
    . ,

and a_l^(2)=(-1)^lΓ( β _1-α _1-l)
/l!( ξ _1^2-α _1-l) ( l+α _1) 
,a_l^(3)=(-1)^lΓ( α _1-β _1-l) /
l!( ξ _1^2-β _1-l) ( l+β _1) .

The proof is provided in Appendix A.





 §.§ R-D Link




  §.§.§ With ZF Precoding


On the other hand, shadowed-Rician fading channel is considered for the
satellite RF links, where the PDF of the SNR γ _ZF^(ϰ
_i) can be calculated by applying Jacobi transform on the fading
envelope PDF given in <cit.> as follows

    f_γ _ZF^(ϰ _i)( z) =λ
    _ϰ _i/γ_ZF^(ϰ _i)e^-ρ
    _ϰ _iz/γ_ZF^(ϰ _i)
    _1F_1( m_s^(ϰ _i);1;δ _ϰ _iz/
    γ_ZF^(ϰ _i)) ,

with_1F_1(.;.;.) denoting the confluent hypergeometric function <cit.>Eq. (9.210)]integrals, and λ _ϰ _i=1/2b_ϰ _i( 
2b_ϰ _im_s^(ϰ _i)/2b_ϰ
_im_s^(ϰ _i)+Ω _s^(ϰ _i))
^m_s^(ϰ _i), ρ _ϰ _i=1/2b_ϰ
_i, δ _ϰ _i=Ω _s^(ϰ _i)/
2b_ϰ _i( 2bm_s^(ϰ _i)+Ω _s^(ϰ
_i)) . Also, Ω _s^(ϰ _i), 2b_ϰ
_i, and m_s^(ϰ _i) stand for the average power of LOS
and multipath components, and the fading severity parameter, respectively.

Consequently, by representing the hypergeometric function _1F_1(
.;.;.) through the finite series <cit.> for integer
parameter values, and using <cit.>, the
respective CDF is expressed as

    F_γ _ZF^(ϰ _i)( z) =λ _ϰ
    _i∑_n=0^m_s^(ϰ _i)-1m_s^(ϰ
    _i)-1nδ _ϰ _i^n/v_ϰ _i^n+1n!
    γ _inc( n+1,v_ϰ _iz/γ
    _ZF^(ϰ _i)) ,

where γ _inc( .,.) stands for the lower incomplete
Gamma function <cit.>, with v_ϰ
_i=ρ _ϰ _i-δ _ϰ _i. Furthermore, the
CDF and PDF of the end-to-end SINR in (<ref>) at e_i^(2) are
expressed as

    F_γ _e_i^(2),ZF( z)     ={[ F_γ _ZF^( e_i^(2)) ( z/ψ
    _i-θ _iz) ;z<ℒ_i;                                   1;z≥ℒ_i ]
    . ,  
    
    f_γ _e_i^(2),ZF( z)     ={[ ψ _i/( ψ _i-θ _iz) ^2∂
    F_γ _ZF^( e_i^(2)) ( z/ψ
    _i-θ _iz) /∂ z;z<ℒ_i;                                                              0;z≥ℒ_i ]
    . .

with ℒ_i=ψ _i/θ _i.



  §.§.§ Without ZF Precoding


Importantly, when no ZF precoding is performed, the CDF of the received
SNR at node ϰ _i, given in (<ref>), can be expressed as

    F_γ _ϰ _i,NZF( z) =[ γ
    _NZF^(ϰ _i)( Ψ ^(ϰ _i)-Θ ^(ϰ
    _i)z) <z] ,

where two cases are distinguished, namely z<ℒ^(ϰ _i)
and z≥ℒ^(ϰ _i), with ℒ^(ϰ
_i)=Ψ ^(ϰ _i)/Θ ^(ϰ _i).
Similarly to the ZF precoding case, the CDF and PDF of the received SNRat node ϰ _i are given as 
    F_γ _ϰ _i,NZF( z)     ={[ F_γ _NZF^(ϰ _i)( z/Ψ ^(ϰ
    _i)-Θ ^(ϰ _i)z) ;z<ℒ^(ϰ _i);                                         1;z≥ℒ^(ϰ _i) ]
    . ,  
    f_γ _ϰ _i,NZF( z)     ={[ Ψ ^(ϰ _i)/( Ψ ^(ϰ _i)-Θ
    ^(ϰ _i)z) ^2∂ F_γ _NZF^(ϰ
    _i)( z/Ψ ^(ϰ _i)-Θ ^(ϰ _i)z
    ) /∂ z;z<ℒ^(ϰ _i);                                                                                    0;z≥ℒ^(ϰ _i) ]
    . .

 with F_γ _NZF^(ϰ _i)( .) is obtained
from (<ref>) by replacing γ_ZF^(ϰ
_i) with γ_NZF^(ϰ _i).



§ SECRECY PERFORMANCE ANALYSIS


The intercept probability metric is defined as the probability that the
secrecy capacity, which is the difference between the capacity of the
legitimate links and that of the wiretap channels (i.e., {S-R, S-e^(1)}, {R-l_i, R-e_i^(2)}), equals to zero.
Additionally, from <cit.> and <cit.>, the
respective channel capacities of the FSO legitimate and wiretap links are
given for coherent detection techniques as

    C_ϱ=log _2( 1+Y_ϱ) ,

with ( ϱ ,Y_ϱ) equals either ( 1,γ
_1) or ( e^(1),γ _1^(e)) .

On the other hand, the system's intercept probability is defined as

    P_int,Ξ^(i)    =( C_s,Ξ^(i)=0)  
        =( . C_s,Ξ^(i)=0|γ _1>γ
    _th) ( γ _1>γ _th)  
        +( . C_s,Ξ^(i)=0|γ _1<γ
    _th) ( γ _1<γ _th) ,

where γ _th is a decoding threshold SNR, below which the decoding
process fails, Ξ ∈{ ZF,NZF}, and

    C_s,Ξ^(i)=min( C_s^(1),C_s,eq,Ξ^(2,i)) ,

with

    C_s^(1)=log _2( 1+γ _1/1+γ _1^(e))
    ,

and

    C_s,eq,Ξ^(2,i)=min( C_s,Ξ^(1,2,i),C_s,Ξ^(2,i)) ,

stand for the end-to-end, first hop, and second hops' secrecy capacities,
respectively, considering DF relaying protocol, with

    C_s,Ξ^(1,2,i)   =   log _2( 1+γ _1/1+γ
    _e_i^(2),Ξ) ,  
    
    C_s,Ξ^(2,i)   =   log _2( 1+γ _l_i,Ξ/1+γ
    _e_i^(2),Ξ) .


Therefore, the overall secrecy capacity in (<ref>) can be expressed as

    C_s,Ξ^(i)=min( C_s^(1),C_s,Ξ^(2,i),C_s,Ξ^(1,2,i)) .



The system's intercept probability for DF relaying scheme is given as

    P_int,Ξ^(i)=1-∫_y=γ _th^∞f_γ _1(
    y) F_γ _1^(e)( y) 𝒥_Ξ(y)dy,

where 𝒥_Ξ(y)=∫_z=0^yf_γ _e_i^(2),Ξ( z) F_γ _l_i,Ξ^c( z) dz, with F_.^c( .) accounts for the complementary CDF [
Both notations γ _l_i,ZF and γ _ZF^( l_i)
stand for the same random variable.].

The proof is provided in Appendix B.





 §.§ Exact Analysis



The integral 𝒥_Ξ(y) in the expression (<ref>) is
given for the ZF and NZF cases as 

    𝒥_ZF(y)   
    =∑_n_1=0^m_s^(l_i)-1∑_n_2=0^m_s^
    ( e_i^(2)) -1𝒰_i( n_1,n_2)
    /exp( r_1^( e_i^(2)) ) ∑
    _k_1=0^n_1∑_j=0^∞(-1)^j/j!k_1!
       ×( ψ _ir_1^(l_i)r_1^( e_i^(2))
    ) ^j+k_1∑_p=0^k_1+n_2+j
    k_1+n_2+jp/( -r_1^( e_i^(2)) )
    ^p-n_2
       ×[
    [ θ _ir_1^( e_i^(2)) Γ(
    p+1-k_1-j,r_1^( e_i^(2)) );   -Γ( p+1-k_1-j,r_1^( e_i^(2)) ψ _i
    /ψ _i-θ _iy) ]
    ] ,


and in (<ref>) at the top of the next page  for y<s_i, with
s_i equals ℒ_i or ℒ^( e_i^(2)) for the ZF and NZF cases, respectively, r_1^(ϰ _i)=
v_ϰ _i/γ_2,ZF^(ϰ _i)θ _i, 𝒢^(a_i,b_i)=v_a_iΨ _i^(b)/𝒯
_iγ_2,NZF^(a_i),( a,b∈{l,e^(2)})
, 𝒯_i=Ψ ^( e_i^(2)) Θ ^(l_i)-Ψ
^(l_i)Θ ^( e_i^(2)) , Γ( .,.)
is the upper-incomplete Gamma function <cit.>, and

𝒰_i( n_1,n_2) =λ
_l_iλ _e_i^(2)/n_2!(
[ m_s^(l_i)-1;         n_1 ]
) (
[ m_s^( e_i^(2)) -1;               n_2 ]
) δ _l_i^n_1/v_l_i^n_1+1δ
_e_i^(2)^n_2/v_e_i^(2)^n_2+1.



It is noteworthy that the system's IP will be computed for the ZF precoding
scenario for two cases, namely ℒ_i<γ _th and ℒ_i≥γ _th, while for the non-ZF scenario, the IPexpression is limited only to the case when ℒ^(
e_i^(2)) <γ _th, due to the toughness of the encountered
mathematical computation.

The proof is provided in Appendix C.




The integral K( φ) =∫_y=φ^∞f_γ
_1( y) F_γ _1^(e)dy can be expressed as follows

    K( φ) =1-F_γ _1( φ) F_γ
    _1^(e)( φ) -ℋ_1+ℋ_2(
    φ) ,

with

    ℋ_1=ξ _e^2/𝒫_1^-K
    ∑_h_1+h_2+h_3=K∑_l=0^∞ℱ
    _h_1,h_2,h_3,l( α _e) _ϱ
    _l,h_1,h_2,h_3( β _e) _ϱ
    _l,h_1,h_2,h_3/( ξ _e^2+ϱ
    _l,h_1,h_2,h_3) Υ _e^ϱ _l,h_1,h_2,h_3
    ,


    ℋ_2( φ)     =𝒫_1^K𝒫
    _e∑_h_1+h_2+h_3=K∑_l=0^∞ϱ
    _l,h_1,h_2,h_3ℱ_h_1,h_2,h_3,l/φ ^-ϱ
    _l,h_1,h_2,h_3
       × G_3,5^3,2( Υ _eφ|[        1-ϱ _l,h_1,h_2,h_3,1;ξ _e^2+1; ξ _e^2,α _e,β _e;0,-ϱ _l,h_1,h_2,h_3 ]
    . ) ,

where ( .) _. stands for the Pochhammer symbol <cit.>, and Υ _e=ξ _e^2α
_eβ _e/μ _1^(e)( ξ _e^2+1) .

The proof is provided in Appendix D.




The system's IP closed-form expression is given in (<ref>) at
the top of the next page, as well as


    P_int,NZF^(i)=1-𝒥_NZF( min[ ℒ
    ^(l_i),ℒ^( e_i^(2)) ] ) 𝒦
    ( γ _th) ,

for ℒ^( e_i^(2)) <γ _th, with

    𝒪_1^(i)=𝒥_ZF( ℒ_i) [
    𝒦( γ _th) -𝒦( ℒ
    _i) ] ,


    𝒪_m+1^(i)    =𝒫_1^K𝒫_e/exp( -r_1^( e_i^(2)) ) 
    ∑_h_1+h_2+h_3=K∑_l=0^∞ℱ
    _h_1,h_2,h_3,lϱ _l,h_1,h_2,h_3
       ×∑_n_1=0^m_s^(l_i)-1∑_n_2=0^m_s^
    ( e_i^(2)) -1𝒰_i( n_1,n_2)
    ∑_k_1=0^n_1∑_j=0^∞(
    r_1^(l_i)ψ _i) ^j+k_1/j!k_1!
       ×∑_p=0^k_1+n_2+jk_1+n_2+jp
    [ ℬ_1^(m,i)+ℬ_α _e^(m,i)+ℬ
    _β _e^(m,i)] /( -r_1^( e_i^(2))
    ) ^p-n_2-k_1-j;m=1,2,

and

    ℬ_ξ _e^2^(m,i)   =   {[ 𝒮_i( ξ _e^2,0) ;m=1; ℛ_i( ξ _e^2,0) ;m=2 ]
    . ,  
    ℬ_x^(m,i)   =   {[ ∑_v=0^∞𝒮_i( x,v) ;m=1; ∑_v=0^∞ℛ_i( x,v) ;m=2 ]
    . ,x∈{α _e,β _e} ,

with b_ξ _e^2^(0)=Γ( α _e-ξ
_e^2) Γ( β _e-ξ _e^2) /ξ _e^2, b_α _e^(v)=(-1)^vΓ( β _e-α
_e-v) /v!( ξ _e^2-α _e-v) ( v+α
_e) ,b_β _e^(v)=(-1)^vΓ( α
_e-β _e-v) /v!( ξ _e^2-β _e-v) (
l+β _e),

𝒞(x,q,v)=(-1)^q( g(x,q,v)-1) q!, g(x,q,v)=1+q+v+x+ϱ _l,h_1,h_2,h_3,

    𝒮_i( y,v)     =Υ _e^y+v/ℒ
    _i^-𝒞(y,0,v)b_x^(v)∑_q=0^∞
    ( 1-ℒ_i/γ _th) ^q+1/q!
       × G_0,3^4,0( . r_1^( e_i^(2))
    γ _th/ℒ_i-γ _th|Λ _1(
    y,q,v) ) ,


    ℛ_i( y,v)     =Υ _e^y+vb_y^(v)/
    ℒ_i^-𝒞(y,0,v)∑_q=0^∞
    ( ℒ_i/γ _th-1) ^-𝒞(y,0,v)-q
    /𝒞(y,q,v)
       × G_0,2^3,0( . r_1^( e_i^(2))
    |Λ _2( y,q,v) ) -Υ
    _e^y+vb_y^(v)Γ( 𝒞(y,0,v)) /ℒ
    _i^-𝒞(α _e,0,v)
       × G_0,1^2,0( r_1^( e_i^(2)) |[     -;g(y,0,v); 0,-k_1-j+p+1;- ]
    . ) ,

Λ _1( x,q,v) =(
[            -;g(x,0,v),1,2+q; g(x,q,v),1+q,0,-k_1-j+p+1;- ]
), and

Λ _2( x,q,v) =(
[            -;g(x,0,v),1; 0,-k_1-j+p+1,g(x,q,v);- ]
)

The proof is provided in Appendix E.





 §.§ Asymptotic Analysis



The system's IP can be asymptotically expanded at high SNR regime (i.e., γ_Ξ^(l_i),μ _1→∞ ) as follows

    P_int,Ξ^( i,∞) =G_c,Ξ( γ
    _Ξ^(l_i)) ^-G_d,

with

    G_c,Ξ    ={[ ( a_0^(d)𝒫_1( ϵ _Ξ^(i)Υ
    ^'a^(Ξ )) ^x_d) ^KF_γ _1^(e)(
    a^(Ξ ));                                       +𝒳( a^(Ξ )) , if Kx_d<1;                                               𝒬^(i,Ξ ),Kx_d>1 ]
    . ,  
    
    a^(Ξ )    ={[ γ _th,Ξ =NZF or ( Ξ =ZF and γ
    _th≥ℒ_i);                ℒ_i,Ξ =ZF and γ _th<ℒ_i ]
    . ,  
    
    G_d    =min( 1,Kx_d) ,

where Υ ^'=ξ _1^2α _1β _1/ξ
_1^2+1, ϵ _Ξ^(i)=γ_Ξ^(l_i)/μ _1, x_d=i=1,2,3min( x_i)
, and

    𝒳( φ)     =𝒫_e( a_0^(d)
    𝒫_1( ϵ _Ξ^(i)Υ ^')
    ^x_d) ^K
       ×[
    [                            Γ( α _e+Kx_d) Γ( β
    _e+Kx_d) /Υ _e^Kx_d( ξ _e^2+Kx_d); -φ ^Kx_dG_2,4^3,1( Υ _eφ|[        1-Kx_d;ξ _e^2+1; ξ _e^2,α _e,β _e;-Kx_d ]
    . ) ]
    ] .

and 𝒬^(i,ZF), 𝒬^(i,NZF) are given in (<ref>)
and (<ref>) at the top of the next page, with 𝒵^(i)=
v_l_iΨ ^( e_i^(2)) /𝒯_iexp(
𝒢^( e_i^(2),l_i) ) , E_i(.) is the
exponential integral function <cit.>.


The proof is provided in Appendix F.





§ NUMERICAL RESULTS


The derived analytic results, which are validated with their respective
Monte Carlo simulations are depicted in this section to analyze the
performance of the considered set up. Some illustrative numerical examples
are depicted to inspect the effects of the key system parameters on the
overall secrecy performance of the considered hybrid terrestrial-satellite
link. To this end, the system and channel parameters are set as detailed in
Table <ref>. The FSO hop turbulence parameters were computed based
on OGS-satellite distance, wavelength, and aperture radius using <cit.>, <cit.>, <cit.>, and <cit.>. Furthermore, the positions of the legitimate and wiretap
nodes were set using Table I and based on (<ref>) with  ϕ
_j^(l_j)=3× 10^-3 rad, ϕ _j^( e_j^(2))
=6.66× 10^-4 rad , respectively, except Figs. 7-11, where ϕ
_j^( e_j^(2)) =6.66× 10^-4 rad, ϕ
_j^(l_j)=5× 10^-4 rad. Also, the simulation is performed by
generating 2× 10^6 random samples.






In Fig. <ref>, the IP is shown as a function of the average SNR P_SAT/σ _ϰ^2 (ϰ ={l,e^(2)}), by
considering P_SAT/σ _l^2=P_SAT/σ
_e^(2)^2, for various optical apertures at the satellite. One can
ascertain evidently that the secrecy performance decreases as a function of P_SAT/σ _ϰ^2, particularly below a certain
threshold SNR and for high values of K. However, the secrecy gets steady
for higher average SNR values. Also, the number of optical apertures at the
satellite impacts the secrecy level of the system. The greater the number of
apertures, the higher is the legitimate link capacity, and consequently, the
more secure is the whole link. Nevertheless, one can note also that above K=7 apertures, the IP is not improved significantly.

Fig. <ref> shows the IP evolution for the ZF case versus P_SAT/σ _ϰ^2 for various values of legitimate
receiver antenna gains. One can remark that for increasing the antennas
gains of the legitimate users can ensure a significant improvement of the
received legitimate SNR, and consequently results in better secrecy
performance of the system.

Fig. <ref> depicts the IP for the non-ZF case as a function of Ω _s^(l_i) for various values b_l_i, with Ω
_s^( e_i^(2)) =3, m_s^(x_i)=2, b_e_i^(2)=1, and P_SAT/σ _l^2=P_SAT/
σ _e^(2)^2=40 dB. One can ascertain that the higher the power
of the legitimate link LOS and multipath components Ω _s^(l_i)
and b_l_i, the higher the second hop's legitimate link SNR, which
results in a greater overall secrecy capacity. Thus, the system's secrecy
gets better.

Fig. <ref> shows the IP evolution versus the power portion received
at the satellite ω _l, for a fixed transmit power-to-noise ratio of
the second hop P_SAT/σ _l^2=P_SAT/σ
_e^(2)^2=40 dB, and P_S/σ _1^2=50, 60, and 70 dB with P_S/σ _1^2=P_S/σ
_e^(1)^2+20 dB. One can notice that the system's secrecy improves by
increasing the portion of power received by the legitimate node (i.e., the
satellite). In fact, as ω _l is increasing, the SNR of the first
hop's legitimate link as well as the equivalent end-to-end SNR increase. As
a result, the overall secrecy capacity in (<ref>) gets greater too,
which leads to a decrease in the system's IP. Nevertheless, at higher ω _l values (higher γ _1), the IP gets steady at a certain
level. In fact, from (<ref>), the overall secrecy capacity will be
restricted to the minimum of the two hops' secrecy capacities (i.e., C_s,eq,Ξ^(2,i)), regardless the increase in ω _l (i.e., C_s^(1) increase).

Figs. <ref> and <ref> depict the IP evolution versus P_SAT/σ _l^2, for the ZF and the non-ZF precoding
scenarios, respectively, for three different scenarios of the first hop
SNRs, namely P_S/σ _1^2=80 dB, and P_S/σ
_e^(1)^2=60, 70,and 80 dB. One can remarkably note that for the
both abovementioned scenarios, the IP is improved with increasing the
difference between the FSO legitimate and wiretapper average SNRs.Nevertheless, the greater the gap, the lesser the secrecy improvement.

In Fig. <ref>, the IP is shown versus P_SAT/σ _l^2
for two scenarios: (i) ZF precoding case, and (ii) non-ZF precoding. We
fixed again P_S/σ _1^2=P_S/σ _e^(1)^2
=30 dB. Importantly, we depict the IP for the two abovementioned
scenarios, as computed in the non-ZF scenario in (<ref>), for
case 1 ( ℒ^(l_i)<ℒ^( e_i^(2))
) and case 2 ( ℒ^(l_i)≥ℒ^(
e_i^(2)) ), based on Table I and (<ref>) with ϕ
_j^( e_j^(2)) =6.66× 10^-4 rad, ϕ
_j^(l_j)=5× 10^-4 rad ( i.e., case 1) and
 ϕ _j^(l_j)=3× 10^-3 rad, ϕ _j^(
e_j^(2)) =6.66× 10^-4 rad ( i.e., 
case 2), respectively. On the other hand, the two
aforementioned cases were implemented for the ZF precoding expression by
adopting the same two abovementioned scenarios of ϕ _j^(ϰ
_i). One can ascertain evidently that the ZF precoding case outperforms
its non-ZF counterpart for case 1. Thus, the closer the legitimate node
to the beams boresight (cells' centers) compared to the wiretap node, the
worse is the secrecy. Additionally, the non-ZF scenario admits a steady IPbehavior at high SNR regime regardless of the increase in P_SAT/
σ _l^2. In fact, from (<ref>), one can notice clearly that at
high γ_NZF^(ϰ _i) values, the received SNRsγ _ϰ _i,NZF reduce to ℒ^(l_i)and ℒ^( e_i^(2)) , for ϰ =l,e^(2),
respectively. Given that ℒ^(l_i)<ℒ^(
e_i^(2)) , it follows that the wiretapper received SNR is
always greater than the legitimate one for high average SNR values, leading
to steady IP performance. On the other hand, the IP shows a distinguished
behavior for the second case ( i.e., ϕ _j^(l_i)≥ϕ _j^( e_i^(2)) ), where the non-ZF case
outperforms the ZF precoding one. In fact, and similarly to the previous
case, as ℒ^(l_i)≥ℒ^( e_i^(2))
, the legitimate SNR is most likely greater than the wiretapper one. Thus,
the farther the legitimate node from the cells' center compared to e_i^(2), the better the secrecy performance. Additionally, due to ZFSNR normalization in (<ref>) and (<ref>), the ZF precoding
average SNR most likely drops below ℒ^(l_i), resulting in
secrecy performance degradation.

Fig. <ref> represents the IP evolution versus μ _1, when ZF
precoding is performed, for various SNRproportionality values (
i.e., μ _1=ϵ _iγ_ZF^(
l_i) ,ϵ _i=1/2,1,2), and for
eavesdropper average SNR value P_S/σ _e^(1)^2=
P_S/σ _e^(2)^2=30 dB. We can ascertain that the higher the
first hop average SNR (i.e., higher OGS transmit power or lower satellite
reception noise), the better the overall system's secrecy. Furthermore, the
asymptotic curves match tightly the exact ones, which proves the accuracy of
the retrieved expressions in high SNR regime.



§ CONCLUSION


The secrecy performance of a hybrid terrestrial-satellite communication
system, operating with an optical feeder in the presence of potential
wiretappers in both hops is assessed. Statistical properties of the per-hop
SNR of the legitimate, as well as the wiretap link, are derived.
Furthermore, a novel expression for the IP for a dual-hop DF-based
communication system is retrieved, based on which the analyzed system's IP
metric is investigated in closed-form and asymptotic expressions in terms of
key system and channel parameters for several parameters' values cases. The
analysis was performed for two scenarios, namely by considering ZFprecoding at the satellite, and by assuming direct satellite delivery to the
earth-stations (i.e., non-ZF). The obtained analytical results show that the
system's secrecy can be significantly improved by increasing the number of
optical photodetectors at the satellite, gateway and satellite transmit
powers, as well as antennas gains. Furthermore, ZF precoding technique
improves the system's secrecy level compared to the non-ZF case for some
specific nodes' positions scenarios.

A potential extension of this work might be the consideration of
amplify-and-forward (AF) relaying scheme as well as investigating the impact
of cloud coverage on the overall secrecy performance.



§ APPENDIX A: PROOF OF PROPOSITION 1


The SNR γ _1 is the maximum among individual SNRs γ
_1^(k) at each of the K photodetectors at the satellite.
Consequently, its respective CDF can be expressed in the case of i.i.d
fading amplitudes, using (<ref>) for coherent detection (i.e., r=1) as

    F_γ _1(z)=[ 𝒫_1G_2,4^3,1( Υ
    z|[         1;ξ _1^2+1; ξ _1^2,α _1,β _1;0 ]
    . ) ] ^K,

with Υ is defined in Proposition 1. By making use of the residues
theorem <cit.> on the Meijer's G-function above, one
obtains

    F_γ _1(z)=𝒫_1^K[ Δ _1+Δ _2+Δ
    _3] ^K,

with

    Δ _i=( Υ z) ^x_iℰ_i,

and ℰ_i=∑_l=0^∞a_l^(i)( Υ z)
^l, and x_i∈{ξ _1^2,α _1,β _1}.

Using the multinomial theorem, the CDF can be formulated as

    F_γ _1(z)=𝒫_1^K∑_h_1+h_2+h_3=KK!/
    h_1!h_2!h_3!( Υ z) ^h_1ξ
    _1^2+h_2α _1+h_3β _1/ℰ_1^-h_1
    ℰ_2^-h_2ℰ_3^-h_3.


Involving the identity <cit.>, we have the following:
ℰ_i^h_i=∑_l=0^∞c_l^(i)z^lΥ
^l, where the coefficients c_l^(i) are defined in Proposition 1.
Consequently, we have the following

    ℰ_1^h_1ℰ_2^h_2ℰ
    _3^h_3=c_1^(1)∑_l=0^∞(
    ∑_q_2+q_3=lc_q_2^(2)c_q_3^(3))
    z^lΥ ^l.


By involving (<ref>) into (<ref>), Proposition 1 is attained.



§ APPENDIX B: PROOF OF LEMMA 1


 Relying on the probability theory, we have the following

    ( C_s,Ξ^(i)=0,γ _1>γ _th) =(
    γ _1>γ _th) -( C_s,Ξ^(i)>0,γ
    _1>γ _th) .


When γ _1<γ _th, the satellite fails at decoding the
information message. Therefore, no signal will be transmitted to the
legitimate earth-stations l_i as well as the eavesdroppers e_i^(2).
 Hence, we have γ _l_i,Ξ=γ _e_i^(2),Ξ=0, which
yields from (<ref>) and (<ref>) that C_s,Ξ^(i)=0 and ( . C_s,Ξ^(i)=0|γ _1<γ
_th) =1. As a result, by making use of (<ref>) into (<ref>), the overall IP expression reduces to

    P_int,Ξ^(i)=1-( C_s,Ξ^(i)>0,γ _1>γ
    _th) .


By involving (<ref>) into the above equation and relying on
probability theory and some algebraic manipulations, one obtains

    P_int,Ξ^(i)=1-ℐ( γ
    _1>γ _1^(e),γ _1>γ _e_i^(2),Ξ,γ
    _l_i,Ξ>γ _e_i^(2),Ξ,γ _1>γ _th) .


The above probability ℐ can be expressed as ℐ
=∑_i=1^6ℐ_i, with I_i are the probabilities
associated with the events E_i shown in Table II.


Relying on Table II, one can see that

    ℐ_1   =∫_y=γ _th^∞F_γ
    _1^c( y) f_γ _1^(e)( y)
    dy∫_z=γ _th^yf_γ _e_i^(2),Ξ( z)
    F_γ _l_i,Ξ^c( z) dz,  
    ℐ_2    =∫_z=γ _th^∞F_γ _1^c(
    z) f_γ _e_i^(2),Ξ( z) F_γ _l_i,Ξ^c( z) ∫_y=γ _th^zf_γ _1^(e)(
    y) dydz,  
    ℐ_3    =∫_x=γ _th^∞∫_y=0^γ
    _th∫_z=y^γ _th∫_t=z^∞f_γ _1(
    x) f_γ _1^(e)( y) f_γ _e_i^(2),Ξ( z)  
       × f_γ _l_i,Ξ( t) dxdydzdt,  
    ℐ_4    =∫_x=γ _th^∞∫_z=0^γ
    _th∫_y=z^γ _th∫_t=z^∞f_γ _1(
    x) f_γ _1^(e)( y) f_γ _e_i^(2),Ξ( z)  
       × f_γ _l_i,Ξ( t) dxdydzdt,  
    ℐ_5    =∫_y=0^γ _th∫_x=γ _th^∞∫_z=γ _th^x∫_t=z^∞f_γ _1( x)
    f_γ _1^(e)( y) f_γ _e_i^(2),Ξ(
    z)  
       × f_γ _2,Ξ^(l_i)( t) dxdydzdt,


    ℐ_6   =∫_x=γ _th^∞∫_y=γ
    _th^x∫_z=0^γ _th∫_t=z^∞f_γ _1(
    x) f_γ _1^(e)( y) f_γ _e_i^(2),Ξ( z)  
       × f_γ _l_i,Ξ( t) dxdydzdt.


By applying an integration by parts on (<𝓇ℯ𝒻>) and (<𝓇ℯ𝒻>) with u^'(y)=f_γ _1^(e)( y) , and
on (<𝓇ℯ𝒻>) with u^'(x)=f_γ _1(
x) alongside with some algebraic manipulations, one obtains

    ℐ_1    =∫_y=γ _th^∞f_γ _1(
    y) F_γ _1^(e)( y) dy∫_z=γ
    _th^yf_γ _e_i^(2),Ξ( z) F_γ
    _l_i,Ξ^c( z) dz  
        -∫_y=γ _th^∞F_γ _1^c( y)
    f_γ _e_i^(2),Ξ( y) F_γ _l_i,Ξ^c( y) F_γ _1^(e)( y) dy.


    ℐ_3=F_γ _1^c( γ _th)
    ∫_y=0^γ _thF_γ _1^(e)( y) f_γ
    _e_i^(2),Ξ( y) F_γ _l_i,Ξ^c(
    y) dy.


    ℐ_5=F_γ _1^(e)( γ _th) {[             ∫_z=γ _th^∞f_γ _e_i^(2),Ξ( z)
    F_γ _2^c( z) dz; -∫_x=γ _th^∞F_γ _1( x) f_γ
    _e_i^(2),Ξ( x) F_γ _2^c( x) dx ]
    }


Additionally, by using the basic definition of the CDF in terms of
the respective PDF, it yields

    ℐ_2    =∫_z=γ _th^∞F_γ _1^c(
    z) f_γ _e_i^(2),Ξ( z) F_γ _l_i,Ξ^c( z) F_γ _1^(e)( z) dz  
        -F_γ _1^(e)( γ _th) ∫_z=γ
    _th^∞F_γ _1^c( z) f_γ
    _e_i^(2),Ξ( z) F_γ _l_i,Ξ^c(
    z) dz.



    ℐ_4    =F_γ _1^c( γ _th) F_γ
    _1^(e)( γ _th) ∫_z=0^γ _thf_γ
    _e_i^(2),Ξ( z) F_γ _l_i,Ξ^c(
    z) dz  
        -F_γ _1^c( γ _th) ∫_z=0^γ
    _thf_γ _e_i^(2),Ξ( z) F_γ _l_i,Ξ^c( z) F_γ _1^(e)( z) dz.



    ℐ_6    =[ ∫_x=γ _th^∞f_γ
    _1( x) F_γ _1^(e)( x) dx]
    ∫_z=0^γ _thf_γ _e_i^(2),Ξ( z)
    F_γ _l_i,Ξ^c( z) dz  
        -F_γ _1^c( γ _th) F_γ _1^(e)(
    γ _th) ∫_z=0^γ _thf_γ _e_i^(2),Ξ( z) F_γ _l_i,Ξ^c( z) dz.


By summing the terms ℐ_3, ℐ_4, and ℐ
_6, one obtains

    ℐ_3+ℐ_4+ℐ_6    =[ ∫_x=γ
    _th^∞f_γ _1( x) F_γ _1^(e)(
    x) dx]   
       ×[ ∫_z=0^γ _thf_γ _e_i^(2),Ξ( z) F_γ _l_i,Ξ^c( z) dz] .


In a similar manner, summing the terms ℐ_1, ℐ_2,
and ℐ_5 gives the following

    ℐ_1+ℐ_2+ℐ_5    =∫_y=γ
    _th^∞f_γ _1( y) F_γ _1^(e)(
    y) dy∫_z=γ _th^yf_γ _e_i^(2),Ξ
       × F_γ _l_i,Ξ^c( z) dz.


Thus, by summing the abovementioned two formulas, (<ref>) is achieved.



§ APPENDIX C: PROOF OF LEMMA 2




 §.§ ZF Case


By involving (<ref>) and (<ref>) alongside with <cit.> and <cit.> into 𝒥_ZF(y) in Lemma 1, one obtains

    𝒥_ZF(y)    =F_γ _e_i^(2),ZF( y)
    -∑_n_1=0^m_s^(l_i)-1∑_n_2=0^m_s^
    ( e_i^(2)) -1𝒰_i( n_1,n_2)
    /n_1!
       ×( 𝒴_e_i^(2)^(ZF))
    ^n_2+1∫_0^yz^n_2γ _inc( n_1+1,𝒴
    _l_i^(ZF)z) /exp( 𝒴_e_i^(2)^(ZF)z
    /ψ _i-θ _iz) ( ψ _i-θ _iz)
    ^n_2+2dz.

with 𝒴_ϰ _i^( Ξ) =v_ϰ
_i/γ_Ξ^( ϰ _i). By using
the lower-incomplete Gamma sum representation in <cit.> alongside with some algebraic manipulations, it yields

    𝒥_ZF(y)   
    =∑_n_1=0^m_s^(l_i)-1∑_n_2=0^m_s^
    ( e_i^(2)) -1𝒰_i( n_1,n_2)
    /( 𝒴_e_i^(2)^(ZF)) ^-n_2-1
    ∑_k_1=0^n_1( 𝒴_l_i^(ZF)) ^k_1
    /k_1!
       ×∫_0^yz^n_2+k_1exp[ -( 1/
    𝒴_e_i^(2)^(ZF)( ψ _i-θ _iz) +
    𝒴_l_i^(ZF)) z] /( ψ _i-θ
    _iz) ^n_2+2dz.


Based on the change of variable w=z/ψ _i-θ _iz+1/
θ _i, we have the following

    𝒥_ZF(y)    =exp( 𝒴_e_i^(2)^(ZF)
    θ _i^-1) /ψ _i∑
    _n_1=0^m_s^(l_i)-1∑_n_2=0^m_s^(
    e_i^(2)) -1𝒰_i( n_1,n_2) /
    ( 𝒴_e_i^(2)^(ZF)) ^-n_2-1
       ×∑_k_1=0^n_1(-1)^k_1/k_1!
    ∑_j=0^∞( ψ _i𝒴
    _l_i^(ZF)) ^j+k_1/j!∑_p=0^k_1+n_2+j
    k_1+n_2+jp
       ×( -θ _i) ^p-n_2∫_1/θ _i^
    y/ψ _i-θ _iy+1/θ _iw^-k_1-j+pexp( -𝒴_e_i^(2)^(ZF)w) dw.


Finally, using the upper-incomplete Gamma function definition in <cit.>, (<ref>) is attained.



 §.§ Non-ZF Case


By incorporating the CDF and PDF expressions from (<ref>) jointly
with <cit.> and <cit.> into 𝒥_NZF(y) in Lemma 1, one obtains

    𝒥_NZF(y)    =F_γ _e_i^(2),NZF( y) -Ψ
    ^( e_i^(2))
    ∑_n_1=0^m_s^(l_i)-1∑_n_2=0^m_s^
    ( e_i^(2)) -1𝒰_i( n_1,n_2)
    /( 𝒴_e_i^(2)^(NZF)) ^-n_2-1
       ×∫_z=0^yz^n_2exp( -𝒴
    _e_i^(2)^(ZF)z/( Ψ ^( e_i^(2)) -Θ
    ^( e_i^(2)) z) ) /( Ψ ^(
    e_i^(2)) -Θ ^( e_i^(2)) z) ^n_2+2
       ×γ _inc( n_1+1,𝒴_l_i^(NZF)z/
    Ψ ^(l_i)-Θ ^(l_i)z) dz.


By transforming the lower-incomplete Gamma function in (<ref>) to
an upper incomplete one alongside with the finite sum expression of the
upper-incomplete Gamma function <cit.>, one can
see that

    𝒥_NZF(y)    =F_γ _e_i^(2),NZF( y)
    -F_γ _e_i^(2),NZF( y) +Ψ ^(
    e_i^(2)) 
       ×∑_n_1=0^m_s^(l_i)-1∑_n_2=0^m_s^
    ( e_i^(2)) -1𝒰_i( n_1,n_2)
    /( 𝒴_e_i^(2)^(NZF)) ^-n_2-1
    ∑_k_1=0^n_1( 𝒴_l_i^(NZF))
    ^k_1/k_1!
       ×∫_z=0^yexp[ -z( 𝒴
    _l_i^(NZF)/Ψ ^(l_i)-Θ ^(l_i)z+𝒴
    _e_i^(2)^(NZF)/Ψ ^( e_i^(2)) -Θ ^(
    e_i^(2)) z) ] /( Ψ ^(
    e_i^(2)) -Θ ^( e_i^(2)) z) ^n_2+2
       ×z^n_2+k_1/( Ψ ^(l_i)-Θ
    ^(l_i)z) ^k_1dz.


At this level, two subcases for the result are distinguished, namely ℒ^( e_i^(2)) >ℒ^(l_i), and ℒ^( e_i^(2)) <ℒ^(l_i). By making
use of the change of variable w={[                  Ψ ^(l_i)𝒯_i/Ψ ^(l_i)-Θ ^(l_i)z
+Ψ ^(l_i)Θ ^( e_i^(2)) , if ℒ
^( e_i^(2)) >ℒ^(l_i); -Ψ ^( e_i^(2)) 𝒯_i/Ψ
^( e_i^(2)) -Θ ^( e_i^(2)) z+Θ
^(l_i)Ψ ^( e_i^(2)) , if ℒ^(
e_i^(2)) <ℒ^(l_i) ]
. alongside with <cit.> and the
binomial theorem, yields (<ref>) at the top of the next page.
Finally, by using <cit.>, 𝒥_NZF(y)
in (<ref>) is attained. 



For the first study subcase (i.e., ℒ^( e_i^(2)) >
ℒ^(l_i)), and from the integral definition 𝒥
_Ξ(y) in Lemma 1, we have the following

    𝒥_NZF( ℒ^( e_i^(2)) )    
    =∫_0^ℒ^(l_i)f_γ _e_i^(2),NZF(
    z) ( 1-F_γ _e_i^(2),NZF( z) ) dz
    
        +∫_ℒ^(l_i)^ℒ^( e_i^(2))
    f_γ _e_i^(2),NZF( z) 0
    F_γ _l_i,NZF^c( y) dz  
        =𝒥_NZF( ℒ^(l_i)) .


Thus, when evaluating 𝒥_NZF( .) for y=ℒ
^( e_i^(2)) when ℒ^( e_i^(2))
>ℒ^(l_i), it should be evaluated at ℒ^(l_i)
within the final result in (<ref>) for the non-ZF case.




§ APPENDIX D: PROOF OF LEMMA 3


Relying on integration by parts in (<ref>), the integral K(
φ) can be expressed as

    K( φ) =1-[ F_γ _1( y) F_γ
    _1^(e)( y) ] _φ^∞-ℋ_1+
    ℋ_2( φ) .

with ℋ_1=∫_y=0^∞F_γ _1( y)
f_γ _1^(e)( y) dy and ℋ_2( φ) =∫_y=0^φF_γ _1( y) f_γ
_1^(e)( y) dy. Involving the CDF and PDF expressions in (
<ref>) with parameters r=1, α _e, β _e, ξ
_e^2, and μ _e, and (<ref>) into the abovementioned
equations, and making use of the Mellin transform <cit.>,
the integral ℋ_1 is computed as

    ℋ_1    =𝒫_e𝒫_1^K
    ∑_h_1+h_2+h_3=K∑_l=0^∞ℱ
    _h_1,h_2,h_3,l
       ×Γ( α _e+ϱ
    _l,h_1,h_2,h_3) Γ( β _e+ϱ
    _l,h_1,h_2,h_3) /Υ _e^ϱ
    _l,h_1,h_2,h_3( ξ _e^2+ϱ
    _l,h_1,h_2,h_3) .


Furthermore, by using <cit.>, the above
equation reduces to (<ref>). Again, by using <cit.> in ℋ_2( φ), (
<ref>) is attained.



§ APPENDIX E: PROOF OF PROPOSITION 2




 §.§ ZF Case


By using Lemma 1 result as well as the PDF and CDF expressions of γ
_e_i^(2),ZF and γ _l_i,ZF in (<ref>) and (<ref>), we distinguish two cases, namely ℒ _i<γ
_th<y and ℒ _i>γ _th.



  * First case ( ℒ_i<γ _th): as
the PDF f_γ _e_i^(2),ZF( z) =0 for z>ℒ
_i, the integral in (<ref>) becomes

    ℐ=𝒦( γ _th) 
    [ ∫_y=γ _th^∞f_γ _1( y)
    F_γ _1^(e)( y) dy] 𝒥
    _ZF( ℒ_i) [ ∫_z=0^
    ℒ_if_γ _e_i^(2),ZF( z) F_γ
    _l_i,ZF^c( z) dz] ,

for ℒ_i<γ _th. Thus, using Lemma 2 result for y=
ℒ_i as well as Lemma 3 result with φ =γ _th and
involving it into (<ref>), one obtains

    P_int,ZF^(i)=1-𝒦( γ _th) 𝒥
    _ZF( ℒ_i) ;ℒ_i<γ _th,



  * Second case ( ℒ_i≥γ _th): in such
an instance, two sub-cases are distinguished, namely ℒ_i≥
y≥γ _th and y>ℒ_i≥γ _th. Hence, the
integral (<ref>) becomes

    ℐ    =𝒪^(i)∫_y=γ
    _th^ℒ_if_γ _1( y) F_γ
    _1^(e)( y) ∫_z=0^yf_γ _e_i^(2),ZF(
    z) F_γ _l_i,ZF^c( z) dzdy
        +𝒩^(i)[ ∫_y=ℒ
    _i^∞f_γ _1( y) F_γ _1^(e)(
    y) dy] [ ∫_z=0^ℒ_if_γ
    _e_i^(2),ZF( z) F_γ _l_i,ZF^c( z)
    dz] .



Interestingly, the integrals of the second term 𝒩^(i)
can be computed readily from Lemma 2 and Lemma 3, with y=ℒ_i
and φ =ℒ_i, respectively. On the other hand, by
involving Lemma 2 result in (<ref>) as well as the derivative of (
<ref>) and (<ref>) with ( ϖ =e, r=1)
into 𝒪^(i), it produces the following

    𝒪^(i)   =   ∫_y=γ _th^ℒ
    _if_γ _1( y) F_γ _1^(e)( y)
    𝒥(y)dy  
       =   𝒪_1^(i)-𝒪( ℒ_i) .

with 𝒪_1^(i) and 𝒪( ℒ_i)
given in (<ref>) and (<ref>) at the top of the next page,
respectively.



Note that 𝒪( ℒ_i)equals 𝒪_2^(i) or 𝒪_3^(i) given in (<ref>) for γ _th≤ℒ_i<2γ _th and ℒ_i>2γ
_th, respectively.

One can notice evidently that (<ref>) yields from (<ref>), (<ref>), (<ref>), as well as (<ref>) and (<ref>). On
the other hand, the Meijer's G-function in (<ref>) can be expressed
using residues theorem as <cit.> 

    G_2,4^3,1( Υ _ey|[         1;ξ _e^2+1; ξ _e^2,α _e,β _e;0 ]
    . )     =b_ξ _e^2^(0)( Υ _ey) ^ξ
    _e^2+∑_v=0^∞b_α _e^(v)( Υ
    _ey) ^α _e+v
        +∑_v=0^∞b_β _e^(v)( Υ _ey)
    ^β _e+v.

under the assumption: x_i-x_j∉

ℤ
,i≠ j, with x_i being defined in Proposition 1.



  * Thus, incorporating the above residues expansion into (<ref>),
making use of the Meijer's G representation of Γ( .,.)
given in <cit.> as well as performing a
change of variable y/ψ _i-θ _iy+1/θ _i=t
produces: ℬ^(1,i)=ℬ_ξ _e^2^(1,i)+ℬ
_α _e^(1,i)+ℬ_β _e^(1,i), where ℬ
_x^(1,i) is defined in (<ref>) for γ _th<ℒ
_i≤ 2γ _th, with

    𝒮_i( y,v)     =ℒ_i/θ _i
    b_y^(v)Υ _e^y∫_t=γ _th/ψ _i-θ
    _iγ _th+1/θ _i^∞t^-g(y,0,v)
       ×G_1,2^2,0( 𝒴_e_i^(2)^(ZF)t
    |[            -;1; 0,-k_1-j+p+1;- ]
    . ) /( ℒ_i( t-1/θ _i
    ) ) ^2-g(y,0,v)dt.


Again, using the Meijer's G-function definition in <cit.> as well as <cit.> and
through some algebraic manipulations, one obtains:

    𝒮_i( y,v)     =ℒ_i^g( y,0,v)
    -1Υ _e^yb_y^(v)/θ _i1/2π i
    ∫_C( γ _th/ψ _i-θ _iγ _th)
    ^-s-1
       ×Γ( s) Γ( -k_1-j+p+1+s) /
    Γ( 2+s) ( 𝒴_e_i^(2)^(ZF))
    ^-s
       ×_2F_1( s+g( y,0,v) ,1+s;2+s,1-
    ℒ_i/γ _th) ds.


with _2F_1( .,.,.;.) denotes the Gauss hypergeometric
function <cit.>. These
last mentioned identities define this function when the absolute value of
the argument 1-ℒ_i/γ _th is either less or
greater than 1. For the former case ( i.e., γ _th<
ℒ_i≤ 2γ _th) , using the function definition
using eq. (07.23.02.0001.01) of <cit.>, and based on the Pochhammer
symbol simplification <cit.> as well as <cit.>Eq. (07.34.02.0001.01)]wolfram , 𝒮_i( y,v)
given for the first case of (<ref>) (i.e., m=1) is attained.

Importantly, when ( i.e., ℒ_i/γ _th
>2), and using the second definition of _2F_1(
.,.,.;z) , z>1 <cit.> alongside
with eqs. (06.10.02.0001.01, 07.34.02.0001.01) of <cit.> and some
algebraic manipulation, 𝒮_i( y,v) in (<ref>)
is substituted by the notation ℛ_i( y,v), where
the resulting expression of ℛ_i( y,v) are obtained
for the second case of (<ref>) (i.e., m=2).




 §.§ Non ZF Case


By using Lemma 1 result alongside with the PDF/CDF of γ
_e_i^(2),NZF and γ _l_i,NZF in (<ref>), we
consider only the case when ℒ^( e_i^(2))
<γ _th<y . As the PDF f_γ _E_2,NZF^(i)(
z) =0 for z≥ℒ^( e_i^(2)), the
integral definition of 𝒥_NZF( y) in Lemma 1 is
positive only for z<ℒ^( e_i^(2)). Thus, in
such an instance, the integral in (<ref>) becomes

    ℐ   =𝒦( γ _th) 
    [ ∫_y=γ _th^∞f_γ _1( y)
    F_γ _1^(e)( y) dy] 
       ×𝒥_NZF( ℒ^(
    e_i^(2)) ) [ ∫_z=0^ℒ
    ^( e_i^(2)) f_γ _e_i^(2),NZF( z)
    F_γ _l_i,NZF^c( z) dz] ;ℒ^(
    e_i^(2))  <γ _th.

Hence, using Lemma 2 result for y=ℒ^( e_i^(2))
and Lemma 3 result with φ =γ _th and involving it into (<ref>), one obtains (<ref>).



§ APPENDIX F: PROOF OF PROPOSITION 3


From (<ref>), the systems IP can be expressed at high SNR (i.e., γ_Ξ^(l_i),μ _1→∞ ,
γ_Ξ^(l_i)=ϵ _Ξ^(i)μ _1,ϵ _Ξ^(i)>0) as 

    P_int,Ξ^( i,∞)    ≃   1-∫_y=γ
    _th^∞f_γ _1^∞( y) F_γ
    _1^(e)( y) dy∫_z=0^yf_γ _e_i^(2),Ξ( z)  
          ×( 1-F_γ _l_i,Ξ^∞( z) )
    dz.


As μ _1→∞ , Υ→ 0 as can be seen
after (<ref>), which yields that Δ _i given in (<ref>
) will be asymptotically represented by considering only least powers of Υ in ℰ_i's expression given right after (<ref>
), corresponding to l=0. Therefore


    Δ _i^(∞ )∼ a_0^(i)( Υ z) ^x_i,

with x_i denotes ξ _1^2, α _1, and β _1, for i=1,2, and 3, respectively. Consequently, the CDF F_γ _1(z) (
<ref>) can be asymptotically approximated using (<ref>) as

    F_γ _1^∞(z)∼( a_0^(d)( Υ z)
    ^x_d𝒫_1) ^K,

with x_d=i=1,2,3min( x_i).

On the other hand, by using the lower-incomplete Gamma expansion <cit.> in (<ref>) and (<ref>), taking
only the least powers of 1/γ_Ξ^(l_i) (
i.e., n_1=0), one obtains

    F_γ _l_i,ZF^∞( z)    ≃   λ
    _l_iv_l_iz/γ_ZF^( l_i) ,
    
    
    F_γ _l_i,NZF^∞( z)    ≃   λ
    _l_iv_l_iz/γ_NZF^(l_i)( Ψ
    ^(l_i)-Θ ^(l_i)z) .




 §.§ ZF case




  §.§.§ First Case: γ _th>ℒ _i


In this case, by making use of integration by parts, (<ref>) is
expressed as


    P_int,ZF^( i,∞)     =1-( ∫_y=γ
    _th^∞f_γ _1^∞( y) F_γ
    _1^(e)( y) dy)  
       ×( 1-∫_z=0^ℒ_if_γ
    _e_i^(2),ZF( z) F_γ _l_i,ZF^∞(
    z) dz) 
        =F_γ _1^∞( γ _th) F_γ
    _1^(e)( γ _th) +∫_y=γ _th^∞F_γ _1^∞( y) f_γ _1^(e)(
    y) dy  
        +∫_z=0^ℒ_if_γ _e_i^(2),ZF( z)
    F_γ _l_i,ZF^∞( z) dz  
        -( F_γ _1^∞( γ _th) F_γ
    _1^(e)( γ _th) +∫_y=γ _th^∞F_γ _1^∞( y) f_γ _1^(e)(
    y) dy)  
       ×∫_z=0^ℒ_if_γ _e_i^(2),ZF(
    z) F_γ _l_i,ZF^∞( z) dz.


By involving (<ref>) and (<ref>) into (<ref>), it
can be seen that the diversity order of the first two terms above is G_d=Kx_d, G_d=1 for the third term, and G_d=Kx_d+1 for the
fourth one. Therefore, the IP will be expanded by either the first two
terms if Kx_d<1, or the third term when Kx_d>1. Hence, the IP is
expressed as

    P_int,ZF^( i,∞) ≃ G_c,ZF( γ
    _ZF^( l_i) ) ^-G_d,γ _th>ℒ_i,

with G_c,ZF and G_d are given in Proposition 3, where 𝒬
^( i,Ξ) =γ_Ξ^( l_i)
∫_z=0^ρf_γ _e_i^(2),Ξ( z) F_γ
_l_i,Ξ^∞( z) dz, with ρ equals ℒ
_i or ℒ^( e_i^(2)) for ZF and NZF
scenarios, respectively.


    𝒳( φ)     =( γ_Ξ^( l_i) /ϵ _Ξ^(i))
    ^Kx_d∫_y=φ^∞F_γ _1^∞( y)
    f_γ _1^(e)( y) dy  
        =ξ _e^2/Γ( α _e) Γ( β
    _e) ( a_0^(d)ξ _1^2( ϵ _Ξ^(i)Υ ^') ^x_d/Γ( α _1)
    Γ( β _1) ) ^K
       ×[
    [  ∫_y=0^∞z^Kx_d-1G_1,3^3,0( Υ _ez|[         -;ξ _e^2+1; ξ _e^2,α _e,β _e;- ]
    . ); -∫_y=0^φz^Kx_d-1G_1,3^3,0( Υ
    _ez|[         -;ξ _e^2+1; ξ _e^2,α _e,β _e;- ]
    . ) ]
    ] .


By using the Mellin transform <cit.> and the identity <cit.>Eq. (07.34.21.0003.01)]wolfram for 𝒳( φ),
respectively, alongside with some manipulations, one obtains (<ref>). On
the other hand, by plugging (<ref>) and (<ref>) into 𝒬^( i,ZF) and by using a change of variable t=ψ
_i/ψ _i-θ _iz with <cit.>, (<ref>) is reached.



  §.§.§ Second Case: γ _th<ℒ _i


Likewise, using integration by parts, the IP in (<ref>) is expressed
at high SNR when γ _th<ℒ_i as given in (<ref>
) at the top of the next page,


where in a similar way, it can be seen that the IP is expanded as given in (
<ref>), (<ref>) and (<ref>) for Ξ =ZF.



 §.§ Non-ZF case


In a similar manner to the ZF case, the IP can be expanded as given in (<ref>) by replacing ℒ_i by ℒ^(
e_i^(2)) as P_int,NZF^( i,∞) ≃
G_c,NZF( γ_NZF^( l_i) )
^-G_d for γ _th>ℒ^( e_i^(2)) , with
G_c,NZF and G_d are defined in Proposition 3. By involving (<ref>) and (<ref>) into 𝒬^( i,NZF)
defined in the previous subsection, it yields (<ref>) at the top of
the next page. Finally, by using the change of variable x={[                     ( 𝒯_iΨ ^(l_i)/Ψ ^(l_i)-Θ
^(l_i)z+Θ ^( e_i^(2)) Ψ ^(l_i))
^-1,ℒ^( e_i^(2)) >ℒ^(l_i); ( -𝒯_iΨ ^( e_i^(2)) /Ψ
^( e_i^(2)) -Θ ^( e_i^(2)) z+Θ
^(l_i)Ψ ^( e_i^(2)) ) ^-1,ℒ
^( e_i^(2)) <ℒ^(l_i) ]
. with the binomial theorem alongside with <cit.> for ℒ^( e_i^(2)) >
ℒ^(l_i), and <cit.>, <cit.> for ℒ^(
e_i^(2)) <ℒ^(l_i), and performing some algebraic
manipulations, one obtains (<ref>) given in the proposition.


IEEEtran


