# Modulating Retroreflector Based Free Space Optical Link for UAV-to-Ground
Communications

M. T. Dabiri, M. Rezaee, L. Mohammadi, F. Javaherian, V. Yazdanian,
M. O. Hasna, Senior Member, IEEE, and M. Uysal, Fellow, IEEE

###### Abstract

Weight reduction and low power consumption are key requirements in the next
generation of unmanned aerial vehicle (UAV) networks. Employing modulating
retro-reflector (MRR)-based free space optical (FSO) technology is an
innovative technique for UAV-to-ground communication in order to reduce the
payload weight and power consumption of UAVs which leads to increased
maneuverability and flight time of UAV. In this paper, we consider an MRR-
based FSO system for UAV-to-ground communication. We will show that the
performance of the considered system is very sensitive to tracking errors.
Therefore, to assess the benefits of MRR-based UAV deployment for FSO
communications, the MRR-based UAV FSO channel is characterized by taking into
account tracking system errors along with UAV’s orientation fluctuations, link
length, UAV’s height, optical beam divergence angle, effective area of MRR,
atmospheric turbulence and optical channel loss in the double-pass channels.
To enable effective performance analysis, tractable and closed-form
expressions are derived for probability density function of end-to-end signal
to noise ratio, outage probability and bit error rate of the considered system
under both weak-to-moderate and moderate-to-strong atmospheric turbulence
conditions. The accuracy of the analytical expressions is verified by
extensive simulations. Analytical results are then used to study the
relationship between the optimal system design and tracking system errors.

###### Index Terms:

Angle of arrival (AoA) fluctuations, FSO communications, UAV, modulating
retro-reflector (MRR).

##  I Introduction

Unmanned aerial vehicles (UAVs) have been introduced to overcome many of the
shortcomings of the current terrestrial infrastructure by operating as aerial
communication nodes and providing robust line-of-sight (LoS) connectivity to
ground devices [1, 2]. UAVs deployed as flying communication nodes using radio
frequency (RF) will interfere with ground devices, hence degrading the
performance of the ground network. Employing UAVs equipped with free space
optical (FSO) technology is a promising method for future ultra dense wireless
networks [3]. FSO communication systems use very small beam divergence which
is physically inaccessible to RF technologies, thus, making FSO links
extremely secure. Moreover, FSO system is much faster, easier to deploy, more
compact, and cheaper than RF [4, 5, 6]. However, vulnerability to signal
blockage is one of the fundamental limitations of FSO links that essentially
confines the receiver to be placed within the line-of-sight (LoS) of the
transmitter. A potential application of UAV-assisted FSO systems is in dense
cities with tall buildings where a UAV can act as an aerial relay to connect
source and destination nodes where the LoS between ground optical nodes is
interrupted by tall buildings [7, 8].

###  I-A Literature Review and Statement of the Problem

Recently, significant research works have been focused on communication
problems of UAV-based FSO systems and how to utilize their vulnerabilities [9,
10, 11, 12, 13, 14, 15, 16, 17, 18]. For instance, two possible scenarios are
proposed in [9] for the integration of UAVs as buffer-aided moving relays into
the conventional relay-assisted FSO systems. In [10], the 3D deployment and
resource allocation of a UAV Base Station with FSO-based backhaul is studied
in a given hotspot area. In [11, 12], the authors design an efficient
algorithm for FSO-based UAVs relay network topology to achieve a high network
reliability. The trajectory optimization of a fixed-wing UAV using FSO
communication is addressed in [13]. In particular, the authors focus on
maximizing the flight time of the UAV by considering practical constraints
including limited propulsion energy and required data rates. The very small
beam divergence which inherently increases the secrecy and capacity of an FSO
link makes it very sensitive to beam misalignment. Unlike an stable ground
node, the position and orientation of UAVs fluctuate due to independent random
effects such as wind speed, changes in the air pressure, propeller rotation,
engine operation, attitude control system faults and platform stability error
[19]. However, the results of [9, 10, 11, 12, 13] are obtained by neglecting
the effect of UAVs’ random fluctuations which are valid for larger UAVs
equipped by fast and precise stabilizers.

More recently, the authors in [14, 15, 16, 17, 18], studied the problem of
UAV-based FSO links by including the effects of UAV’s orientation and position
fluctuations. In [14], the authors analyzed the end-to-end system performance
of mixed RF-FSO networks employing UAVs as buffer-aided and non-buffer-aided
relays in terms of the ergodic sum rate by taking into account the impacts of
UAV’s fluctuations. A novel FSO channel model between a hovering UAV and a
central unit is developed in [15] by quantifying the corresponding geometric
and misalignment losses, while taking into account the non-orthogonality of
the laser beam. Considering the joint effects of UAVs’ fluctuations as well as
atmospheric turbulence, a novel channel model was proposed in [16, 17] that is
suitable for hovering UAV-based FSO links with zero boresight angle. In [18],
the authors completed the results of [16, 17] and provided a novel UAV-based
FSO channel model that takes into account the effect of nonzero boresight
pointing errors. In [15], the considered value for standard deviation (SD) of
orientation fluctuations is $\sigma_{\theta_{o}}=0.2-0.8$ mrad (or
equivalently is $\sigma_{\theta_{o}}=0.0115^{o}-0.0458^{o}$ in degrees), in
[17] is $\sigma_{\theta_{o}}=5$ mrad, and in [18] is $\sigma_{\theta_{o}}=1-5$
mrad. However, due to the flight time, payload, and power consumption
limitations of lightweight multi-rotor UAVs, reaching such values for SD of
orientation fluctuations may not be always possible. Another challenge is an
asymmetric nature of the ground-to-UAV and the UAV-to-ground links. Although
the ground-to-UAV link can withstand larger orientation fluctuations, the
results of [16, 7] show that severe fluctuations of UAVs greatly reduce the
quality of the UAV-to-ground link. The ground station (GS) has more power and
payload capacity than the UAV node and employing a fast and precise tracking
system in GS node is feasible. Therefore, it can well track the UAV node. In a
ground-to-UAV link, the transmitter (Tx) is mounted on an stable GS node with
precise tracking system, thus, by increasing UAV’s orientation fluctuations,
only the SD of angle-of-arrival (AoA) increases in a ground-to-UAV link and
the receiver’s (Rx’s) field-of-view (FoV) can be increased to relax this
degrading effect [20]. However, in a UAV-to-ground link, the Tx is mounted on
an unstable UAV. For better understanding, consider a 1000 m UAV-to-ground
link. When the orientation of Tx deviates more than 10 mrad, the received
optical beam deviates more than 10 m from the center of Rx’s aperture, and
this leads to an unreliable communication. Moreover, for compensation of UAV’s
orientation fluctuations in a UAV-to-ground link, the Tx needs a power
amplifier to increase the transmitted power which increases the payload and
power consumption of UAV and limits the maneuverability and flight time of
UAV.

Modulating retro-reflectors (MRR)111An MRR consists of an optical
retroreflector with a modulator to first modulate the incoming optical signal
and then reflect it toward the transmitter. This feature makes it possible to
act as an optical communications device without sending its own optical power
as graphically depicted in Fig. 2. A number of technologies have been
considered for the modulation component, including electro-optic modulators,
liquid crystal modulators, multiple quantum well devices, and actuated
micromirrors. The modulator tries to block the reflected signal intensity for
a bit ”0” and tries to pass all the reflected signals for a bit ”1”. In other
words, the modulator changes the intensity of the reflected optical signal in
proportion to the transmitted On–off keying (OOK) signal sequence, which is
called switching speed or switching rate. is potentially attractive in
asymmetric situations such as small UAV platforms which are too small to carry
a conventional FSO terminal [21].222Here, the GS has more power and payload
capacity than the UAV node, and GS sends light towards the remote Rx mounted
on UAV. The UAV is equipped with a small MRR which upon sensing the incoming
interrogating beam, modulates and reflects it directly back to the GS [22].
Retro-reflector links are used in limited duplex communication where Rx have
low power to support full transceiver operations. MRRs are also used to reduce
the pointing and tracking requirements by directly reflecting the incoming
light to the GS independent of orientation of the retro-reflector [22]. In the
MRR-based topology, complexity and tracking equipment are transferred from the
UAV node to the GS. Even though MRR-based optical wireless communications have
been well studied in the context of underwater optical wireless communications
[22], this subject is restricted to few works in the context of UAV-based FSO
communications [23, 24, 25, 26, 27, 28, 29, 30]. In [23, 24], the authors
design and implement a real-time localization and tracking system for a UAV-
to-ground FSO link. In [25], a general geometrical model of the corner cube
reflector (CCR) is established based on the ray tracing method and then, the
authors used the Wave Optics simulations to investigate the double-pass
channel in the MRR-based FSO systems. The probability density function (PDF)
and cumulative distribution function (CDF) of the double-pass MRR-based FSO
systems are derived in [26] under weak turbulence conditions and in [27] under
strong turbulence conditions. Impacts of shape and size of MRR cell along with
turbulence condition on the beam spot are investigated in [28] by the
theoretical analysis and wave-optics simulation. In [29, 30], the effects of
the parameters such as atmospheric turbulence conditions, link length, Rx’s
aperture diameter and the average signal-to-noise ratio (SNR) are studied on
the performance of MRR-based FSO links. However, the results of [25, 26, 27,
28, 29, 30] are obtained by neglecting the effect of geometrical pointing
errors which are only valid when the GS node perfectly pointed optical beam
towards the aperture of MMR mounted on UAV.

Even though a larger value of MRR’s aperture improves the link budget, in a
practical implementation, we are not allowed to use large values for MRR’s
aperture because the switching rate of MRR modulator is inversely proportional
to MRR’s aperture. This leads to a large geometrical loss in an MRR-based FSO
system with respect to the conventional FSO systems with much larger aperture
area. To compensate this problem, the beamwidth must be chosen much smaller
than the values of the beamwidth used in conventional FSO systems. This makes
the MRR-based FSO system very sensitive to tracking errors. Depending on link
length, we will show that any tracking angle error in the order of $\micro$rad
can significantly affect the performance of MRR-based FSO systems. Therefore,
to assess the benefits of MRR-based UAV deployment for FSO communications,
performance analyses of the considered system under tracking system errors is
very important and necessary. To the best of authors’ knowledge, there is no
prior work in the literature that models and analyzes MRR-based FSO systems
for UAVs under tracking system errors by taking into account the effects of
MRR’s orientation fluctuations as well as atmospheric turbulence conditions.

(a)

(b)

Figure 1: The schematic of the considered communication link between an
optical GS and an hovering UAV equipped with an MRR. The $z$ axis refers to
the direction that extends from GS toward UAV node and the GS tries to adjust
optical beam in the direction of $z$-axis. MRR mounted on UAV modulates and
reflects the incoming beam directly back to the GS. (b) Our GS-to-UAV-to-GS
channel model takes into account seven impairments.

###  I-B Major Contributions and Novelty

In this paper, we consider the MRR-based FSO system for UAV-to-ground
communication where the GS sends an unmodulated continuous laser beam towards
the MRR mounted on UAV. The incident beam is modulated by the MRR and is
directly reflected back to the GS. The main contribution of this paper is the
performance analysis and system design of the considered UAV-based FSO system
when UAV is equipped with MRR under tracking errors by taking into account the
effects of MRR’s orientation fluctuations as well as atmospheric turbulence
conditions. In summary, our key contributions include:

  * •

We develop channel models for the UAV-to-ground MRR-based FSO communication
system for both weak-to-moderate and moderate-to-strong atmospheric turbulence
conditions, by taking into account tracking errors, UAV’s orientation
fluctuations, link length, UAV’s height, optical beam divergence angle,
effective area of MRR, atmospheric turbulence and optical channel loss in the
double-pass channels.

  * •

Based on these models, we derive closed-form analytical expressions for PDF
under both weak-to-moderate and moderate-to-strong atmospheric turbulence
conditions. Then, through Monte Carlo simulations, the accuracy of the derived
statistical distributions is verified.

  * •

We also derive the closed-form expressions for the PDF and CDF of end-to-end
signal-to-noise ratio (SNR), outage probability and bit error rate (BER) of
the MRR-based UAV FSO system. The accuracy of the analytical expressions is
verified by using simulations. Analytical results are then used to study the
impact of the system parameters on the performance of MRR-based UAV FSO system
under different conditions, e.g., a wide range of tracking errors, atmospheric
turbulence strengths, different levels of UAV’s instability, different link
lengths, etc. Our results reveal that unlike conventional FSO systems, optimal
design of a UAV-based MRR FSO system is very important and any change in the
parameters (such as link length, SD of tracking errors, target BER, desired
data rate, etc.) affect the optimal values of other parameters.

###  I-C Organization

We list the main notations in Table I. The organization of the rest of the
paper is as follows. In Section II, we characterize the actual channel models
of MRR-based UAV FSO system. Then, in Section III, we provide the analytical
channel models. Next, in Section IV, we provide the simulation results to
verify the derived analytical channel models and study the link performance
and system parameter optimization. Finally, conclusions are drawn in Section
V.

TABLE I: The list of main notations. Parameter | Description
---|---
$[x;y;z]$ | Cartesian coordinate system that $z$ axis refers to the
| direction that extends from GS toward UAV node
$[x^{\prime};y^{\prime};z^{\prime}]$ | Cartesian coordinate system that indicates the
| coordinates of three perpendicular mirrors of MRR
$R$ | PD responsivity
$\Upsilon_{\textrm{th}}$ | SNR threshold
$r_{g}$ | Radius of GS aperture
$A_{r}$ | MRR effective area
$Z$ | Link length
$\sigma_{\theta_{o}}$ | SD of UAV orientation fluctuations
$\sigma_{\theta_{e}}$ | SD of tracking angle errors
$\lambda$ | Wavelength
$\theta_{\textrm{div}}$ | Divergence angle
$w_{z}$ | Beamwidth at the Rx
$P_{t}$ | Transmit power
$\sigma_{n}^{2}$ | Noise variance
$N$ | Number of sectors
$h_{l_{gu}}$ | Channel Loss of GS to UAV link
$h_{l_{ug}}$ | Channel Loss UAV to GS link
$h_{a_{gu}}$ | Atmospheric turbulence coefficient of GS to UAV link
$h_{a_{ug}}$ | Atmospheric turbulence coefficient of UAV to GS link
$h_{p_{u}}$ | Attenuation due to pointing errors at the MRR aperture
$h_{p_{g}}$ | Attenuation due to geometric loss in GS
$h_{MRR}$ | The ratio of direct reflected power by MRR
$C_{n}^{2}$ | Refractive-index structure
$f_{x}(x)$ | The PDF of RV $x$
$F_{x}(x)$ | The CDF of RV $x$
$Q(\cdot)$ | The Q-function defined in [31]
$\textrm{erf}(\cdot)$ | The error function defined in [31]
$\textrm{erfc}(\cdot)$ | The complementary error function defined in [31]
$G^{m,n}_{p,q}\\!\\!\left(\\!z\Big{|}\\!\\!\\!\\!\begin{array}[]{c}-\\\ -\end{array}\\!\\!\\!\right)$ | The Meijer’s G-function defined in [32]

##  II System Model

Figure 2: An MRR using a transmissive device is illustrated where the
interrogation beam is modulated and directly reflected toward the incoming
direction.

Figure 1a shows MRR-based FSO link between an optical GS and an hovering UAV
equipped with an MRR. We assume that the GS is located at $[0;0;0]$ (in
Cartesian coordinate system $[x;y;z]$) and directs a continuous laser
interrogator beam towards the UAV (located at $[0;0;Z]$), in which $Z$ is the
link length between GS and UAV. As illustrated in Fig. 2, MRR modulates and
reflects the incoming beam back to the GS. MRR is suitable for the small UAVs
with low power and payload limitations. It is used to reduce the pointing and
tracking requirements by retro-reflecting the modulated light back to the
interrogating source. On the other hand, the GS is equipped with a precise
laser tracking system that points a continuous laser interrogation beam
towards the UAV. The accuracy of a tracking system is evaluated by the SD of
angle errors. As depicted in Fig. 3, let
$\theta_{ex}\sim\mathcal{N}(0,\sigma_{\theta_{e}}^{2})$ and
$\theta_{ey}\sim\mathcal{N}(0,\sigma_{\theta_{e}}^{2})$ denote tracking system
angle errors in the directions of $x$ and $y$ axes, respectively, with
$\sigma_{\theta_{e}}$ being the SD of the angle errors. These tracking errors
cause a radial distance between the received beam center and the MRR aperture
center as $d_{p}=\sqrt{d_{px}^{2}+d_{py}^{2}}$, where $d_{px}$ and $d_{py}$
are the distance in the directions of $x$ and $y$ axes, respectively. Based on
Fig. 3, any tracking error in the direction of $x$ and $y$ is formulated as

$\displaystyle\sin\left(\theta_{ex}\right)=\frac{d_{px}}{\hat{Z}+z_{e}},~{}~{}~{}\&~{}~{}~{}\sin\left(\theta_{ey}\right)=\frac{d_{py}}{\hat{Z}+z_{e}}.$
(1)

where $z_{e}=Z-\hat{Z}$, and $\hat{Z}$ is the estimated value of $Z$. In
practice, $Z_{e}$ is in the order of a few tens of cm, however, $Z$ is in the
order of a few hundred meters to a few km ($Z_{e}<<Z$), and thus, with a good
accuracy, we can approximate (1) as

$\displaystyle\sin\left(\theta_{ex}\right)=\frac{d_{px}}{Z}\simeq\frac{d_{px}}{\hat{Z}},~{}~{}~{}\&~{}~{}~{}\sin\left(\theta_{ey}\right)=\frac{d_{py}}{Z}\simeq\frac{d_{py}}{\hat{Z}}.$
(2)

We consider a Gaussian beam at the GS, for which the normalized spatial
distribution of the received intensity at distance $Z$, is given by [33]

$\displaystyle I_{r}(d,Z)=\frac{2}{\pi
w_{z}^{2}}\exp\left(-\frac{2(x^{2}+y^{2})}{w_{z}^{2}}\right),$ (3)

where $d=[x,y]$ is the radial distance vector from the beam center. Also,
$w_{z}$ is the beamwidth at distance $Z$ and can be approximated as
$w_{z}\simeq\theta_{div}Z$, where $\theta_{\textrm{div}}$ is the optical beam
divergence angle [34].

Figure 3: Graphical illustration of a ground-to-UAV FSO link wherein the
ground station tracks the UAV mobility with the angle error
$\sigma_{\theta_{e}}$.

The considered MRR consists of three perpendicular triangular mirrors which
has an aperture with a triangular effective area as depicted in Figs. 4 and
5a. For considered MRR-based topology with effective area $A_{r}$, the
effective channel coefficient due to geometric spread with pointing error
$d_{p}=\sqrt{d_{px}^{2}+d_{py}^{2}}$ is obtained as [34]

$\displaystyle h_{p_{u}}$ $\displaystyle=\frac{2}{\pi
w_{z}^{2}}\int\int_{p_{A_{r}}(x,y)}$
$\displaystyle~{}~{}~{}\times\exp\left(-\frac{2((x-d_{px})^{2}+(y-d_{py})^{2})}{w_{z}^{2}}\right)\textrm{d}x\textrm{d}y,$
(4)

where $p_{A_{r}}(x,y)$ is the position of effective aperture area of MRR in
$x-y$ plane. It can be shown that $h_{p_{u}}$ is proportional to the MRR
aperture. This dependence suggests the use of large $A_{r}$ to increase the
link budget. However, the switching rate of MRR modulator is inversely
proportional to $A_{r}$. Accordingly, in practice, the active area of MRR is
usually less than 1 $\textrm{cm}^{2}$ which is much smaller than the beamwidth
$w_{z}$. From this point, the beam is approximately plane in the aperture of
MRR and thus, Eq. (II) can be well approximated as

$\displaystyle h_{p_{u}}$ $\displaystyle=\frac{2A_{r}}{\pi
w_{z}^{2}}\exp\left(-\frac{2(d_{px}^{2}+d_{py}^{2})}{w_{z}^{2}}\right).$ (5)

The vector perpendicular to the aperture area is denoted by $r_{z}$. Here, we
assume that the UAVs orient themselves and/or use a simple servo motor (which
has much lower weight and price than a stabilizer) to set the vector $r_{z}$
in the direction of $z$ axis.333It is assumed that the the positions of GS and
UAV are known to UAV, which can be realized through periodic data exchange
between them and/or by using the global navigation satellite systems. However,
due to the inherent orientation fluctuations of UAVs, the instantaneous
orientation of the vector $r_{z}$ is deviated from $z$ axis. Let $\theta$
denote the instantaneous orientation deviation between vector $r_{z}$ and $z$
axes. As depicted in Fig. 5, a portion of the power collected by the MRR
aperture is scattered depending on $\theta$. The ratio of the power reflected
directly to all of the collected optical power by the MRR is denoted by
$h_{MRR}$. The PDF of $h_{MRR}$ will be derived in the next section. Then, the
reflected optical power is directly back to the GS. It is further assumed that
the GS uses a circular aperture for transmission and reception. We neglect the
effect of beam wandering. This assumption is valid for FSO communication with
link lengths up to several kilometers [35]. Therefore, the beam is received at
the GS with a negligible deviation (due to negligible beam wandering) and thus
the center of the received beam is located approximately at the center of the
GS aperture. Hence, the attenuation due to geometric loss in GS aperture is

$\displaystyle
h_{p_{g}}=\frac{2r_{g}^{2}}{w_{zg}^{2}}\simeq\frac{2r_{g}^{2}}{Z^{2}\theta_{\textrm{div}}^{2}},$
(6)

where $r_{g}$ is the radius of GS aperture, and $w_{zg}$ is the optical
beamwidth in the GS.

In addition to the aforementioned factors, the optical power launched from the
GS is also affected by atmospheric turbulence induced fading and atmospheric
loss before arriving back at the GS. In this study, a bistatic channel is
assumed where the transmitter and the receiver of GS are separated in space by
more than a Fresnel zone, such that, the round trip channel is modeled as a
product of two independent turbulence channels444MRR-based FSO communication
systems can be implemented in monostatic or bistatic configurations. For a
monostatic system, the transmitter and receiver are colocated, whereby, there
is a high correlation between the instantaneous atmospheric coefficients of
the forward and backward passes, which decreases the performance,
significantly. For instance, from the results of [27, Fig. 4], there is a SNR
gap of more than 20 dB between a bistatic configuration with a correlation
coefficient near to zero and a monostatic configuration with a correlation
coefficient equal to 0.6 for a target BER of $10^{-8}$. On the other hand, the
Fresnel zone of FSO links is in the order of one centimeter which is typically
much smaller than optical half power (3 dB) beamwidth at the receiver, and
thus, we can separate the transmitter and receiver of GS more than the Fresnel
zone of FSO link to have a bistatic channel which causes an additional
geometrical loss much less than 3 dB. However, bistatic configuration offers
the advantage of preserving independent forward and backward paths [36]. More
importantly, beam-tracking system impairs for correlated FSO fading channels
[37]. In this paper, we will show that the performance of UAV-based MRR FSO
links significantly depends on the accuracy of the ground beam-tracking
system. Therefore, bistatic is a preferable configuration in most of the work
in the context of MMR-based FSO communications (for instance see [29, 38, 39,
40, 41, 42, 43, 44]). In practice, the interrogator with a bistatic
configuration is also adopted by Naval Research Laboratory to conduct a
communication experiment over 16 km FSO links from the Chesapeake Bay to a
modulated RR array [45, 46].[36]. Therefore, the considered system model
consists of two paths and thus we have two independent atmospheric turbulence
induced fading effects and two independent atmospheric attenuation effects.
Let $h_{a_{gu}}$ and $h_{a_{ug}}$ denote the instantaneous atmospheric
turbulence coefficients of GS-to-UAV and UAV-to-GS, respectively. Further, let
$h_{l_{gu}}$ and $h_{l_{ug}}$ denote the atmospheric attenuation of GS-to-UAV
and UAV-to-GS, respectively. The atmospheric attenuation is typically modeled
by the Beer-Lambert law as [34]

$\displaystyle h_{l_{gu}}=h_{l_{ug}}=\exp(-Z\zeta),$ (7)

where $\zeta$ is the scattering coefficient and is a function of visibility.
Notice that $h_{l_{gu}}$ and $h_{l_{ug}}$ are equal since the parameters
$\zeta$ and $Z$ are the same for both links. The log-normal (LN) and Gamma-
Gamma (GG) distributions are good candidates to efficiently model weak to
moderate and moderate to strong ranges of atmospheric turbulence conditions
[34, 36] Under weak to moderate turbulence conditions, we use LN distribution
to model $h_{a_{gu}}$ and $h_{a_{ug}}$ as [36]

$\displaystyle
f_{L}\left(h_{a_{n}}\right)=\frac{1}{2h_{a_{n}}\sqrt{2\pi\sigma_{L}^{2}}}\exp\left(-\frac{\left(\ln(h_{a_{n}})+2\sigma_{L}^{2}\right)^{2}}{8\sigma_{L}^{2}}\right),$
(8)

where $n\in\\{ub,bu\\}$, $\sigma_{L}^{2}\simeq\sigma_{R}^{2}/4$ is the SD of
the log-normal distribution and $\sigma_{R}^{2}$ is the Rytov variance which
can be obtained for two nodes with different heights as [36, p. 509]

$\displaystyle\sigma_{R}^{2}=9\left({2\pi}/{\lambda}\right)^{7/6}\left(Z/Z_{h_{d}}\right)^{11/6}$
(9)
$\displaystyle~{}~{}~{}\times\int_{Z_{h_{g}}}^{Z_{h_{u}}}C_{n}^{2}(Z_{h})\left(1-\frac{Z_{h}-Z_{h_{g}}}{Z_{h_{d}}}\right)^{5/6}(Z_{h}-Z_{h_{g}})^{5/6}\textrm{d}Z_{h}.$

In (9), $Z_{h_{u}}$ and $Z_{h_{g}}$ denote the height of UAV and GS,
respectively, and $Z_{h_{d}}=Z_{h_{u}}-Z_{h_{g}}$ denote the difference height
between the UAV and GS, and

$\displaystyle
C_{n}^{2}(Z_{h})=0.00594(\mathcal{V}/27)^{2}\left(10^{-5}Z_{h}\right)^{10}\exp\left(-\frac{Z_{h}}{1000}\right)$
$\displaystyle~{}~{}~{}+2.7\times
10^{-16}\exp\left(-\frac{Z_{h}}{1500}\right)+C_{n}^{2}(0)\exp\left(-\frac{Z_{h}}{100}\right),$

is the refractive-index structure parameter at height $Z_{h}$ which
characterizes the atmospheric turbulence, $\mathcal{V}$ (in m/s) is the speed
of strong wind and $C_{n}^{2}(0)$ (in $\textrm{m}^{-2/3}$) is a strong nominal
ground turbulence level. Under moderate to strong turbulence conditions, we
use GG distribution to model the random variables (RVs) $h_{a_{gu}}$ and
$h_{a_{ug}}$ as [36]

$\displaystyle
f_{G}\left(h_{a_{n}}\right)=\frac{2(\alpha\beta)^{\frac{\alpha+\beta}{2}}}{\Gamma(\alpha)\Gamma(\beta)}h_{a_{n}}^{{}^{\frac{\alpha+\beta}{2}-1}}k_{\alpha-\beta}\left(2\sqrt{\alpha\beta
h_{a_{n}}}\right),$ (10)

where $\Gamma(\cdot)$ is the Gamma function and $k_{m}(\cdot)$ is the modified
Bessel function of the second kind of order $m$. Also, $\alpha$ and $\beta$
are respectively the effective number of large-scale and small-scale eddies,
which depend on Rytov variance $\sigma_{R}^{2}$ [36].

From the aforementioned results and as graphically depicted in Fig. 1b, the
instantaneous GS-to-UAV-to-GS channel coefficient is formulated as

$\displaystyle
h=h_{l_{ug}}h_{l_{gu}}h_{a_{ug}}h_{a_{gu}}h_{p_{u}}h_{p_{g}}h_{MRR}.$ (11)

Finally, the instantaneous end-to-end SNR is obtained as [34]

$\displaystyle\Upsilon=\frac{2R^{2}P_{t}^{2}h^{2}}{\sigma_{n}^{2}},$ (12)

where $P_{t}$ denotes the transmitted optical power, $R$ denotes the photo-
detector responsivity, and $\sigma_{n}^{2}$ is the variance of additive
thermal noise.

Figure 4: Gaussian beam footprint at the MRR with the effective triangular
aperture area. The tracking errors cause a radial distance between the
received beam center and the MRR aperture center as
$d_{p}=\sqrt{d_{px}^{2}+d_{py}^{2}}$, where $d_{px}$ and $d_{py}$ are the
distance in the directions of $x$ and $y$ axes, respectively.

##  III Channel Modeling

The considered MRR consists of three perpendicular triangular mirrors that are
located in Coordinate planes $x^{\prime}-y^{\prime}$, $x^{\prime}-z^{\prime}$,
and $y^{\prime}-z^{\prime}$, as depicted in Fig. 5a where
$[x^{\prime};y^{\prime};z^{\prime}]$ is a new Cartesian coordinate system. Let
the vectors $r_{xy}$, $r_{xz}$, and $r_{yz}$ denote the vectors perpendicular
to the aperture area in $x^{\prime}-y^{\prime}$, $x^{\prime}-z^{\prime}$, and
$y^{\prime}-z^{\prime}$, respectively. Also, the angles between the incident
laser beam and the vectors $r_{xy}$, $r_{xz}$, and $r_{yz}$ in
$x^{\prime}-y^{\prime}$, $x^{\prime}-z^{\prime}$, and $y^{\prime}-z^{\prime}$
are denoted respectively by
$\theta_{xy}\sim\mathcal{N}\left(0,\sigma_{\theta_{o}}^{2}\right)$,
$\theta_{xz}\sim\mathcal{N}\left(0,\sigma_{\theta_{o}}^{2}\right)$, and
$\theta_{yz}\sim\mathcal{N}\left(0,\sigma_{\theta_{o}}^{2}\right)$ as depicted
in Fig. 5b. A fraction of the collected laser power by aperture is directly
reflected to the GS which is a function of RVs $\theta_{xy}$, $\theta_{xz}$,
and $\theta_{yz}$. Due to the independence of RVs $\theta_{xy}$,
$\theta_{xz}$, and $\theta_{yz}$, we have

$\displaystyle
h_{MRR}=h_{MRR_{xy}}(\theta_{xy})\,h_{MRR_{xz}}(\theta_{xz})\,h_{MRR_{yz}}(\theta_{yz}),$
(13)

where $h_{MRR_{xy}}$, $h_{MRR_{xz}}$, and $h_{MRR_{yz}}$ are the ratio of
direct reflected optical power to all of the collected optical power by MRR in
planes $x^{\prime}-y^{\prime}$, $x^{\prime}-z^{\prime}$, and
$y^{\prime}-z^{\prime}$, respectively. For a better understanding, and for a
given $\theta_{xz}$ in $x^{\prime}-z^{\prime}$ plane, we have graphically
shown in Figs. 5b and 5c which category of input optical signal is directly
reflected (denoted by $L_{a,xz}$) and which category is scattered (denoted by
$L_{s,xz}$). From this, we have

$\displaystyle
h_{MRR_{n}}(\theta_{n})=\frac{L_{a,n}}{L_{a,n}+L_{s,n}}=1-\tan(\theta_{n}).$
(14)

for $n\in\\{xy,xz,yz\\}$. From (13) and (14), finding a closed-form expression
for the PDF of $h_{MRR}$ is very difficult if not impossible. However, for
system analysis, we need to know the PDF of $h_{MRR}$ denoted by
$f_{h_{MRR}}(h_{MRR})$. Accordingly, in the sequel, $f_{h_{MRR}}(h_{MRR})$ is
approximated by two simple models based on the mean and variance of $h_{MRR}$.
In Fig. 6, the distribution of $h_{MRR}$ is obtained by using simulations for
different values of $\sigma_{\theta_{o}}$. As shown in Fig. 6, the mean and
variance of $h_{MRR}$ (denoted by $\mu_{MRR}$ and $\sigma_{MRR}^{2}$,
respectively) are functions of the angular instability of the UAV
characterized by $\sigma_{\theta_{o}}$. In Table II, $\mu_{MRR}$ and
$\sigma_{MRR}^{2}$ of $h_{MRR}$ are obtained by using simulations for
different values of $\sigma_{\theta_{o}}$. For the rest of
$\sigma_{\theta_{o}}$ values, the corresponding mean and variance values can
be obtained by interpolation from the given values. For instance, from the
results of Table II, for $\sigma_{t_{o}}=1^{o}$ and $\sigma_{t_{o}}=3^{o}$ we
have $\mu_{MRR}=0.95$ and $\mu_{MRR}=0.85$, respectively. Now, for
$\sigma_{t_{o}}=2^{o}$, by using interpolation, we obtain
$\mu_{MRR}\simeq\frac{0.95+0.85}{2}=0.9$.

(a)

(b)

(c)

Figure 5: 3D illustration of MRR configuration: (a) the MRR consists of three perpendicular triangular mirrors that are located in Coordinate planes $x^{\prime}-y^{\prime}$, $x^{\prime}-z^{\prime}$, and $y^{\prime}-z^{\prime}$; (b) the angles between the incident laser beam and the vectors perpendicular to the aperture area in $x^{\prime}-z^{\prime}$ plane; (c) A fraction of collected laser power by aperture is directly reflected to the GS which is a function of RVs $\theta_{xy}$, $\theta_{xz}$, and $\theta_{yz}$.  Figure 6: The distribution of $h_{MRR}$ obtained using simulation for different values of $\sigma_{\theta_{o}}$.  Table II: Mean and SD values of RV $h_{MRR}$ obtained using simulation for different values of $\sigma_{\theta_{o}}$.  $\sigma_{\theta_{o}}$ | $1^{o}$ | $2^{o}$ | $3^{o}$ | $4^{o}$ | $5^{o}$ | $6^{o}$ | $7^{o}$ | $8^{o}$ | $9^{o}$ | $10^{o}$ | $11^{o}$
---|---|---|---|---|---|---|---|---|---|---|---
$\mu_{MRR}$ | 0.96 | 0.93 | 0.89 | 0.86 | 0.83 | 0.8 | 0.76 | 0.73 | 0.70 | 0.66 | 0.62
$\sigma_{MRR}$ | 0.0178 | .035 | 0.052 | 0.066 | 0.083 | 0.094 | 0.11 | 0.12 | 0.13 | 0.145 | 0.158

###  III-A Weak Turbulence Conditions

We use the log-normal distribution to model atmospheric turbulence induced
fading that is valid for weak to moderate turbulence conditions.

Theorem 1. The distribution of $h$ under weak to moderate atmospheric
turbulence conditions is derived as

$\displaystyle
f_{h}(h)=C_{4}h^{K-1}Q\left(\frac{\ln(h)+C_{5}}{\sqrt{C_{1}}}\right),$ (15)

where

$\displaystyle\left\\{\begin{array}[]{ll}&\\!\\!\\!\\!\\!\\!\\!C_{1}=\ln\left(1+\frac{\sigma_{MRR}^{2}}{\mu_{MRR}^{2}}\right)+8\sigma_{L}^{2},\\\
&\\!\\!\\!\\!\\!\\!\\!C_{2}=\ln\left(\frac{\sqrt{\mu_{MRR}^{2}+\sigma_{MRR}^{2}}}{\mu_{MRR}^{2}}\right)+4\sigma_{L}^{2},\\\
&\\!\\!\\!\\!\\!\\!\\!C_{3}=\frac{\pi w_{z}^{2}}{2A_{r}h_{c}},\\\
&\\!\\!\\!\\!\\!\\!\\!C_{4}=KC_{3}^{K}\exp\left(\frac{C_{1}K^{2}-2KC_{2}}{2}\right),\\\
&\\!\\!\\!\\!\\!\\!\\!C_{5}=\ln(C_{3})+C_{1}K+C_{2},\\\
&\\!\\!\\!\\!\\!\\!\\!K=\frac{w_{z}^{2}}{Z^{2}\sigma_{\theta_{e}}^{2}},~{}~{}~{}\&~{}~{}~{}h_{c}=h_{l_{ug}}h_{l_{gu}}h_{p_{g}},\end{array}\right.$
(22)

and $Q(\cdot)$ is the well-known Q-function.

###### Proof:

Please refer to Appendix A. ∎

As we will see in the next section, the accuracy of the derived analytical
channel model will be validated by employing Monte Carlo simulations.
Moreover, the channel model provided in (15) is very simple and tractable
which allows us to easily analyze the performance of the considered MRR-based
FSO system without performing time consuming simulations. Next, the BER and
outage probability of the considered system under weak to moderate atmospheric
turbulence conditions are derived.

Proposition 1. The CDF of $h$ under weak to moderate atmospheric turbulence
conditions is derived as

$\displaystyle F_{h}(h)$
$\displaystyle=\frac{C_{4}}{K}e^{-KC_{5}}\left[e^{K\ln(h)+KC_{5}}Q\left(\frac{\ln(h)+C_{5}}{\sqrt{C_{1}}}\right)\right.$
$\displaystyle~{}~{}~{}\left.+~{}e^{K^{2}C_{1}/2}Q\left(\frac{KC_{1}-\ln(h)-C_{5}}{\sqrt{C_{1}}}\right)\right].$
(23)

###### Proof:

The CDF of $h$ is defined as

$\displaystyle F_{h}(h)=\int_{0}^{h}f_{h}(x)dx.$ (24)

Using (15), we have

$\displaystyle
F_{h}(h)=C_{4}\int_{0}^{h}x^{K-1}Q\left(\frac{\ln(x)+C_{5}}{\sqrt{C_{1}}}\right)dx.$
(25)

Applying a change of variable $y=\ln(x)+C_{5}$ and using (25) and [47, eq.
(06.27.21.0011.01)], the closed-form expression for CDF of $h$ is derived in
(III-A). ∎

Proposition 2. The PDF and CDF of end-to-end SNR under weak to moderate
atmospheric turbulence conditions are derived as

$\displaystyle
f_{\Upsilon}(\Upsilon)=\frac{C_{4}}{2\Upsilon_{1}^{K/2}}\Upsilon^{K/2-1}Q\left(\\!\frac{\ln(\Upsilon)-\ln(\Upsilon_{1})\\!+\\!2C_{5}}{2\sqrt{C_{1}}}\\!\right),$
(26)

and

$\displaystyle F_{\Upsilon}(\Upsilon)$
$\displaystyle=\frac{C_{4}}{K}e^{-KC_{5}}\left[Q\left(\frac{\ln(\Upsilon)-\ln(\Upsilon_{1})+2C_{5}}{2\sqrt{C_{1}}}\right)\right.$
(27)
$\displaystyle~{}\times\exp\left(\frac{K\ln(\Upsilon)-K\ln(\Upsilon_{1})+2KC_{5}}{2}\right)$
$\displaystyle~{}\left.+~{}e^{K^{2}C_{1}/2}Q\left(\frac{2KC_{1}-\ln(\Upsilon)-\ln(\Upsilon_{1})-2C_{5}}{2\sqrt{C_{1}}}\right)\right],$

where $\Upsilon_{1}=\frac{2R^{2}P_{t}^{2}}{\sigma_{n}^{2}}$.

###### Proof:

Based on (12), we have

$\displaystyle
F_{\Upsilon}(\Upsilon)=\textrm{Prob}\\!\left\\{\\!\frac{2R^{2}P_{t}^{2}}{\sigma_{n}^{2}}h^{2}<\Upsilon\\!\right\\}=\textrm{Prob}\\!\left\\{\\!\Upsilon_{1}h^{2}<\\!\Upsilon\\!\right\\}.$
(28)

Using (III-A) and (28), the CDF of $\Upsilon$ is derived in (27). Now, by
differentiating (27) with respect to $\Upsilon$ and after some manipulations,
the closed-form analytical expression for the PDF of $\Upsilon$ is derived in
(26). ∎

In slow fading channels, which is a valid assumption for channel condition in
FSO links [5], the outage probability (i.e., the probability that the
instantaneous end-to-end SNR falls below a threshold $\Upsilon_{\textrm{th}}$)
is the most appropriate performance metric. It is formulated as
$\mathbb{P}_{\textrm{out}}=F_{\Upsilon}(\Upsilon_{\textrm{th}})$. Accordingly,
the results of Proposition 2 let us to easily compute outage probability by
substituting $\Upsilon_{\textrm{th}}$ instead of $\Upsilon$ in (27).

Proposition 3. The closed-form expression for BER of the considered MRR-based
FSO communication for On–off keying (OOK) modulation over weak to moderate
turbulence conditions is derived as

$\displaystyle\mathbb{P}_{e}\simeq
L_{1}L_{2}\Bigg{[}\frac{1}{b}\left(e^{\frac{b^{2}}{2a^{2}}}\textrm{erfc}\left(\frac{b}{2a}-aL_{3}\right)+e^{bL_{3}}\textrm{erfc}(aL_{3})\right)$
$\displaystyle-\sum_{m=0}^{M}\frac{L_{1m}}{b_{m}}\left(e^{\frac{b_{m}^{2}}{2a^{2}}}\textrm{erfc}\left(\frac{b_{m}}{2a}-aL_{3}\right)+e^{b_{m}L_{3}}\textrm{erfc}(aL_{3})\right)\Bigg{]},$
(29)

$\displaystyle\left\\{\begin{array}[]{ll}&\\!\\!\\!\\!\\!\\!\\!L_{1}=\frac{C_{4}}{4\Upsilon_{1}^{K/2}},\\\
&\\!\\!\\!\\!\\!\\!\\!L_{2}=\frac{1}{2}\exp\left(\frac{K(\ln(\Upsilon_{1})-2C_{5})}{2}\right),\\\
&\\!\\!\\!\\!\\!\\!\\!L_{3}=\ln(\Upsilon_{\textrm{max}})-\ln(\Upsilon_{1})+2C_{5},\\\
&\\!\\!\\!\\!\\!\\!\\!L_{1m}=\frac{(-1)^{m}}{m!\sqrt{\pi}(2m+1)2^{m-\frac{1}{2}}}\exp\left(\frac{(2m+1)(\ln(\Upsilon_{1})-2C_{5})}{2}\right),\\\
&\\!\\!\\!\\!\\!\\!\\!a=\frac{1}{2\sqrt{2C_{1}}},~{}~{}~{}\&~{}~{}~{}b=\frac{K}{2},~{}~{}~{}\&~{}~{}~{}b_{m}=\frac{K+2m+1}{2}.\end{array}\right.$
(35)

###### Proof:

Please refer to Appendix B. ∎

The accuracy of the approximated BER provided in (III-A) increases by
increasing $M$ and $\Upsilon_{\textrm{max}}$. Notice that the term of
$e^{\frac{b_{m}^{2}}{2a^{2}}}$ in (III-A) becomes uncountably large by
increasing $M$. for large values of $M$ and $\Upsilon_{\textrm{max}}$.
Therefore, due to the computational limitations, we will limit our simulations
to moderate values of $M$ and $\Upsilon_{\textrm{max}}$. After an exhaustive
search, we found that the BER is obtained with good accuracy for $M=20$ and
$\Upsilon_{\textrm{max}}=4$.

###  III-B Moderate to Strong Turbulence Conditions

We use the GG distribution to model atmospheric turbulence induced fading that
is valid for moderate to strong turbulence conditions.

Theorem 2. The PDF and CDF of $h$ under moderate to strong atmospheric
turbulence conditions are derived as

$\displaystyle f_{h}(h)=B_{s}B_{n}\sum_{n=1}^{N}$ (36)
$\displaystyle\left[G^{6,0}_{2,6}\left(B_{n}^{\prime}h\bigg{|}\begin{array}[]{c}K,1\\\
0,\alpha-1,\beta-1,K-1,\alpha-1,\beta-1\end{array}\right)\right.$ (39)
$\displaystyle-\left.G^{6,0}_{2,6}\left(B_{n}^{\prime\prime}h\bigg{|}\begin{array}[]{c}K,1\\\
0,\alpha-1,\beta-1,K-1,\alpha-1,\beta-1\end{array}\right)\right],$ (42)

and

$\displaystyle F_{h}(h)=B_{s}B_{n}h\sum_{n=1}^{N}$ (43)
$\displaystyle\left[G^{6,1}_{3,7}\left(B_{n}^{\prime}h\bigg{|}\begin{array}[]{c}0,K,1\\\
0,\alpha-1,\beta-1,K-1,\alpha-1,\beta-1,-1\end{array}\right)\right.$ (46)
$\displaystyle-\left.G^{6,1}_{3,7}\left(B_{n}^{\prime\prime}h\bigg{|}\begin{array}[]{c}0,K,1\\\
0,\alpha\\!-\\!1,\beta\\!-\\!1,K-1,\alpha-1,\beta-1,-1\end{array}\right)\right],$
(49)

where

$V_{n}=\left\\{\begin{array}[]{ll}&\\!\\!\\!\\!\\!\\!\\!2\mu_{MRR}-1,~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}n=1,\\\
&\\!\\!\\!\\!\\!\\!\\!V_{n-1}+\frac{(2-2\mu_{MRR})n}{N},~{}~{}~{}n\in\\{2,...,N+1\\},\end{array}\right.$

and $U(x)=\left\\{\begin{array}[]{ll}&\\!\\!\\!\\!\\!\\!\\!1,~{}~{}~{}x>0,\\\
&\\!\\!\\!\\!\\!\\!\\!0,~{}~{}~{}x<0,\end{array}\right.$ is the Heaviside step
function, $B_{n}^{\prime}=\frac{\pi
w_{z}^{2}\alpha^{2}\beta^{2}}{2A_{r}h_{c}V_{n+1}}$,
$B_{n}^{\prime\prime}=\frac{\pi
w_{z}^{2}\alpha^{2}\beta^{2}}{2A_{r}h_{c}V_{n}}$, and $B_{s}=\frac{K\left(\pi
w_{z}^{2}\alpha^{2}\beta^{2}\right)}{(\Gamma(\alpha)\Gamma(\beta))^{2}\left(2A_{r}h_{c}\right)}$.

###### Proof:

Please refer to Appendix C. ∎

The coefficients $B_{n}$ in (36) depend on the parameters
$\sigma_{\theta_{o}}$ and the number of sectors denoted by $N$. Obviously, the
accuracy of the proposed model in (36) directly depends on the number of
sectors $N$, and for sufficiently large values of $N$, an exact match between
simulations and analysis can be achieved at the cost of higher complexity.
Hence, choosing an optimal value for $N$ involves a trade-off between
complexity and accuracy. We will show, via simulations, that $N=8$ achieves
sufficient accuracy. In Table III, the coefficients $B_{n}$ obtained using
simulation for different values of $\sigma_{\theta_{o}}$ and $N=8$. For the
rest of $\sigma_{\theta_{o}}$, the corresponding coefficients $B_{n}$ can be
obtained by interpolation.

Table III: Coefficient $B_{n}$ for $N=8$ obtained using simulation for different values of $\sigma_{\theta_{o}}$.  $\sigma_{\theta_{o}}$ | $1^{o}$ | $3^{o}$ | $5^{o}$ | $7^{o}$ | $9^{o}$ | $11^{o}$
---|---|---|---|---|---|---
$B_{1}$ | 2.63 | 0.85 | 0.47 | 0.29 | 0.19 | 0.1
$B_{2}$ | 5.74 | 1.99 | 1.24 | 0.91 | 0.72 | 0.58
$B_{3}$ | 10.37 | 3.73 | 2.42 | 1.87 | 1.58 | 1.42
$B_{4}$ | 15.2 | 5.49 | 3.56 | 2.75 | 2.3 | 2.07
$B_{5}$ | 17.8 | 6.19 | 3.9 | 2.94 | 2.4 | 2.06
$B_{6}$ | 14.7 | 4.99 | 3.03 | 2.2 | 1.73 | 1.42
$B_{7}$ | 7.05 | 2.45 | 1.44 | 1 | 0.76 | 0.6
$B_{8}$ | 1.26 | 0.4 | 0.23 | 0.15 | 0.11 | 0.08

Proposition 4. The PDF and CDF of end-to-end SNR under moderate to strong
atmospheric turbulence conditions are derived as

$\displaystyle
f_{\Upsilon}(\Upsilon)=\frac{B_{s}B_{n}}{2\sqrt{\Upsilon_{1}\Upsilon}}\sum_{n=1}^{N}$
(50)
$\displaystyle\left[G^{6,0}_{2,6}\left(B_{n}^{\prime}\sqrt{\frac{\Upsilon}{\Upsilon_{1}}}\bigg{|}\begin{array}[]{c}K,1\\\
0,\alpha-1,\beta-1,K-1,\alpha-1,\beta-1\end{array}\right)\right.$ (53)
$\displaystyle-\left.G^{6,0}_{2,6}\left(B_{n}^{\prime\prime}\sqrt{\frac{\Upsilon}{\Upsilon_{1}}}\bigg{|}\begin{array}[]{c}K,1\\\
0,\alpha\\!-\\!1,\beta\\!-\\!1,K-1,\alpha-1,\beta-1\end{array}\right)\right],$
(56)

and

$\displaystyle
F_{\Upsilon}(\Upsilon)=B_{s}B_{n}\sqrt{\frac{\Upsilon}{\Upsilon_{1}}}\sum_{n=1}^{N}$
(57)
$\displaystyle\left[G^{6,1}_{3,7}\left(B_{n}^{\prime}\sqrt{\frac{\Upsilon}{\Upsilon_{1}}}\bigg{|}\begin{array}[]{c}0,K,1\\\
0,\alpha-1,\beta-1,K-1,\alpha-1,\beta-1,-1\end{array}\right)\right.$ (60)
$\displaystyle-\\!\left.G^{6,1}_{3,7}\left(\\!B_{n}^{\prime\prime}\sqrt{\frac{\Upsilon}{\Upsilon_{1}}}\bigg{|}\begin{array}[]{c}0,K,1\\\
0,\alpha\\!-\\!1,\beta\\!-\\!1,K\\!-\\!1,\alpha\\!-\\!1,\beta\\!-\\!1,-\\!1\end{array}\\!\\!\right)\\!\right],$
(63)

###### Proof:

Please refer to Appendix C. ∎

(a)

(b)

(c)

(d)

Figure 7: Comparison of the accuracy of the channel models under weak
turbulence conditions, when $C_{n}^{2}=5\times 10^{-15}$, $Z=1000$ m, and
$\sigma_{\theta_{e}}=100~{}\micro$rad and for (a) channel PDF given in Theorem
1, (b) CDF of channel coefficients given in Proposition 1, (c) PDF of SNR
given in Proposition 2, and (d) CDF of SNR given in Proposition 2.

Proposition 5. For small values of $\sigma_{\theta_{o}}$, (36) and (43) can be
simplified respectively as

$\displaystyle f_{h}(h)\simeq\frac{\pi
w_{z}^{2}\alpha^{2}\beta^{2}K}{2A_{r}(\Gamma(\alpha)\Gamma(\beta))^{2}}\times$
(64) $\displaystyle G^{5,0}_{1,5}\left(\frac{\pi
w_{z}^{2}\alpha^{2}\beta^{2}}{2A_{r}}h\bigg{|}\begin{array}[]{c}K\\\
\alpha-1,\beta-1,K-1,\alpha-1,\beta-1\end{array}\right),$ (67)

and

$\displaystyle F_{h}(h)\simeq\frac{\pi
w_{z}^{2}\alpha^{2}\beta^{2}K}{2A_{r}(\Gamma(\alpha)\Gamma(\beta))^{2}}h\times$
(68) $\displaystyle G^{5,1}_{2,6}\left(\\!\\!\frac{\pi
w_{z}^{2}\alpha^{2}\beta^{2}}{2A_{r}}h\bigg{|}\\!\\!\begin{array}[]{c}0,K\\\
\alpha-1,\beta-1,K-1,\alpha-1,\beta-1,-1\end{array}\\!\\!\right).$ (71)

Also, for small values of $\sigma_{\theta_{o}}$, (50) and (57) can be
simplified respectively as

$\displaystyle f_{\Upsilon}(\Upsilon)\simeq\frac{\pi
w_{z}^{2}\alpha^{2}\beta^{2}K}{2A_{r}(\Gamma(\alpha)\Gamma(\beta))^{2}}\frac{1}{2\sqrt{\Upsilon_{1}\Upsilon}}\times$
(72) $\displaystyle G^{5,0}_{1,5}\left(\frac{\pi
w_{z}^{2}\alpha^{2}\beta^{2}}{2A_{r}}\sqrt{\frac{\Upsilon}{\Upsilon_{1}}}\bigg{|}\begin{array}[]{c}K\\\
\alpha\\!-\\!1,\beta\\!-\\!1,K\\!-\\!1,\alpha-1,\beta-1\end{array}\right),$
(75)

and

$\displaystyle F_{\Upsilon}(\Upsilon)\simeq\frac{\pi
w_{z}^{2}\alpha^{2}\beta^{2}K}{2A_{r}(\Gamma(\alpha)\Gamma(\beta))^{2}}\sqrt{\frac{\Upsilon}{\Upsilon_{1}}}\times$
(76) $\displaystyle G^{5,1}_{2,6}\left(\\!\\!\frac{\pi
w_{z}^{2}\alpha^{2}\beta^{2}}{2A_{r}}\sqrt{\frac{\Upsilon}{\Upsilon_{1}}}\bigg{|}\\!\\!\begin{array}[]{c}0,K\\\
\alpha-1,\beta-1,K-1,\alpha-1,\beta-1,-1\end{array}\\!\\!\right).$ (79)

###### Proof:

Please refer to Appendix C. ∎

In the next section, we show that for UAVs with higher stability, i.e.,
$\sigma_{\theta_{o}}<2^{o}$, the results of Proposition 5 is accurate.

Proposition 6. The closed-form expression for BER of the considered MRR-based
FSO communication for OOK modulation over moderate to strong turbulence
conditions is derived as

$\displaystyle\mathbb{P}_{e}=\frac{B_{s}}{4\sqrt{\pi\Upsilon_{1}\Upsilon}}\sum_{n=1}^{N}B_{n}\frac{2^{3\alpha+2\beta}}{256\pi^{2}}\left[G^{12,2}_{6,13}\left(\frac{(B_{n}^{\prime})^{2}}{32\Upsilon_{1}}\bigg{|}\begin{array}[]{c}<\\!a_{i}\\!>\\\
<\\!b_{j}\\!>\end{array}\right)\right.$ (82)
$\displaystyle-\left.G^{12,2}_{6,13}\left(\frac{(B_{n}^{\prime\prime})^{2}}{32\Upsilon_{1}}\bigg{|}\begin{array}[]{c}<\\!a_{i}\\!>\\\
<\\!b_{j}\\!>\end{array}\right)\right]$ (85)

where

$\displaystyle\left\\{\begin{array}[]{ll}&\\!\\!\\!\\!\\!\\!\\!<\\!a_{i}\\!>=0,\frac{1}{2},\frac{K}{2},\frac{K+1}{2},\frac{1}{2},1\\\
&\\!\\!\\!\\!\\!\\!\\!<\\!b_{j}\\!>=0,\frac{1}{2},\frac{\alpha-1}{2},\frac{\alpha}{2},\frac{\beta-1}{2},\frac{\beta}{2},\frac{K-1}{2},\frac{K}{2},\frac{\alpha-1}{2},\frac{\alpha}{2},\frac{\beta-1}{2},\frac{\beta}{2},-1.\end{array}\right.$
(87)

###### Proof:

Please refer to Appendix D. ∎

TABLE IV: System Parameters for Simulations. Description | Parameter | Setting
---|---|---
PD responsivity | $R$ | 0.8 A/W
SNR threshold | $\Upsilon_{\textrm{th}}$ | 5 dB
Radius of GS aperture | $r_{g}$ | 8 cm
MRR effective area | $A_{r}$ | 0.5-4 $\textrm{cm}^{2}$
Link length | $Z$ | 500-1500 m
SD of UAV orientation fluctuations | $\sigma_{\theta_{o}}$ | 1-8 in degree
SD of tracking angle errors | $\sigma_{\theta_{e}}$ | 50-400 $\micro$rad
Wavelength | $\lambda$ |  $1550$ nm
Divergence angle | $\theta_{\textrm{div}}$ | 0.1-2 mrad
Beamwidth at the Rx | $w_{z}$ | 0.1-2 m
Transmit power | $P_{t}$ | 0-30 dBm
Noise variance | $\sigma_{n}^{2}$ | -11 dBm
Number of sectors | $N$ | 8
Channel Loss | $h_{l_{gu}}~{}\&~{}h_{l_{ug}}$ | 0.7
Refractive-index structure | $C_{n}^{2}$ | $10^{-14}-10^{-13}$
for moderate-to-strong turbulence |  |
Refractive-index structure | $C_{n}^{2}$ | $10^{-15}-10^{-14}$
for weak-to-moderate turbulence |  |
Figure 8:  Comparison of the accuracy of the channel models provided in
Theorem 2 and Proposition 5 under moderate to strong turbulence conditions,
when $Z=1000$ m, and $\sigma_{\theta_{e}}=90~{}\micro$rad. Figure 9:  Outage
probability versus $P_{t}$ for $Z=1000$, $\sigma_{\theta_{o}}=6^{o}$,
$\sigma_{\theta_{e}}=100~{}\micro$rad, and three different values of
$C_{n}^{2}=10^{-14}$, $5\times 10^{-14}$, and $10^{-13}$.  Figure 10:  Outage
probability versus $P_{t}$ for $Z=1000$, $\sigma_{\theta_{o}}=6^{o}$,
$\sigma_{\theta_{e}}=100~{}\micro$rad, $w_{z}=40$ cm, and different values of
effective aperture area $A_{r}$.

##  IV Numerical Results

In this section, first, we utilize computer simulations to verify the accuracy
of our proposed analytical channel models for MRR-based FSO system. Second,
the performance of the considered system is studied in terms of outage
probability and BER. The main considered parameters for the simulation results
are summarized in Table IV, mostly adopted from the standard values for system
parameters in [34].

For evaluation of analytical channel models provided in Section III, we
perform Monte-Carlo simulations. The details of the simulation process are
described as follows. For a given $\sigma_{\theta_{o}}$, we generate $5\times
10^{7}$ independent RVs $\theta_{xz}$, $\theta_{xy}$ and $\theta_{yz}$. Then,
based on (14), we generate $5\times 10^{7}$ independent coefficients of
$h_{MRR_{n}}$ for $n\in\\{xy,xz,yz\\}$. Now, using (13), $5\times 10^{7}$
independent coefficients of $h_{MRR}$ are generated. Moreover, for a given
$\sigma_{\theta_{e}}$, we generate $5\times 10^{7}$ independent RVs
$\theta_{ex}$ and $\theta_{ey}$. Then, using the generated RVs $\theta_{ex}$
and $\theta_{ey}$, we generate $5\times 10^{7}$ independent coefficients of
$h_{pu}$ from (2) and (II). For a given $C_{n}^{2}<10^{-14}$, we also generate
$5\times 10^{7}$ independent coefficients of $h_{a_{gu}}$ and $h_{a_{ug}}$
which have log-normal distribution as given in (8). For a given
$C_{n}^{2}>10^{-14}$, we generate $5\times 10^{7}$ independent coefficients of
$h_{a}$ which have GG distribution as given in (10). We then obtain $5\times
10^{7}$ independent values of UAV-based optical channel coefficients based on
(11). Finally, we find the channel distribution diagrams.

The accuracy of the proposed channel models under weak turbulence conditions
is evaluated in Fig. 7 for two different values of the UAV’s instability
parameters $\sigma_{\theta_{o}}=2^{o}$ and $8^{o}$. In particular, we
corroborate the accuracy of PDF of instantaneous channel coefficients and end-
to-end SNR (respectively provided in Theorem 1 and Proposition 2) in Figs. 7a
and 7c, respectively, and we corroborate the accuracy of CDF of instantaneous
channel coefficients and end-to-end SNR (respectively provided in Propositions
1 and 2) in Figs. 7b and 7d, respectively. The simulation results clearly
confirm the accuracy of analytical models under weak turbulence conditions. In
Fig. 8, we also investigate the accuracy of PDF of instantaneous channel
coefficients under strong turbulence conditions for two different values of
the UAV’s instability parameters $\sigma_{\theta_{o}}=2^{o}$ and $8^{o}$. From
Theorem 2, the parameter $N$ impacts on the validity of channel PDF. The
variable $N$ is used for approximating the distribution of RV $h_{MRR}$. The
optimal value for $N$ is its minimum value that satisfies a predefined
accuracy. The results of Fig. 8 clearly show that the analytical channel model
derived in Theorem 2 with $N=8$ is valid for all conditions. In Proposition 5,
we also propose a more tractable closed-form channel model. As we observe, the
analytical channel model derived in proposition 5 is accurate for more stable
UAVs with $\sigma_{\theta_{o}}<2^{o}$.

The effect of atmospheric turbulence conditions on optical link performance
can be characterized by the index of refraction structure parameter
$C_{n}^{2}$. In Fig. 9, outage probability is plotted for $Z=1000$,
$\sigma_{\theta_{o}}=6^{o}$, $\sigma_{\theta_{e}}=100~{}\micro$rad, and a wide
range of $C_{n}^{2}$ i.e., $C_{n}^{2}=10^{-14}$, $5\times 10^{-14}$, and
$10^{-13}$. The results of Fig. 9 clearly show the impairments caused by
atmospheric turbulence on the performance of MRR-based FSO system. The effect
of atmospheric turbulence on the considered MRR-based FSO system is more
severe than the conventional FSO system because in MRR FSO system, the
atmospheric turbulence affects both GS-to-UAV and UAV-to-GS links
characterized by $h_{a_{gu}}$ and $h_{a_{ug}}$, respectively.

The effective area of MRR denoted by $A_{r}$ is another important parameter
that plays a key role in the link budget. In Fig. 10, we investigate the
effect of $A_{r}$ on the performance of MRR FSO system. A larger value of
$A_{r}$ lets the MRR to collect more optical power and thus, improves the link
budget which leads to a lower outage probability. However, in a practical
implementation, we are not allowed to use large values for $A_{r}$ because the
switching rate of MRR modulator is inversely proportional to $A_{r}$.
Accordingly, in practice, the active area of MRR is usually selected by a
trade-off between link budget and desired data rate.

Figure 11:  Average BER of MRR-based FSO system for $w_{z}=30$ cm, $Z=1000$ m,
and different values of $\sigma_{\theta_{e}}$ and $\sigma_{\theta_{o}}$.

Other main parameters for a UAV-based MRR system are the SD of tracking system
errors and the SD of UAV orientation fluctuations which are denoted by
$\sigma_{\theta_{e}}$ and $\sigma_{\theta_{o}}$, respectively. These
parameters have a significant impact on system performance and related link
budget. To get a better insight, in Fig. 11, the BER of the considered system
is plotted for different values of $\sigma_{\theta_{e}}$ and
$\sigma_{\theta_{o}}$. From the results of Fig. 11, at the target BER
$\mathbb{P}_{e}=10^{-6}$ and $\sigma_{\theta_{e}}=100$ $\micro$rad, the
considered system requires approximately 3 dB more transmit power to
compensate the degrading effect of increasing UAV’s instability from
$\sigma_{\theta_{o}}=2^{o}$ to $6^{o}$. However, for
$\sigma_{\theta_{o}}=2^{o}$, by increasing tracking system error from
$\sigma_{\theta_{e}}=100$ to 200 $\micro$rad, the performance significantly
degrades. As a result, for an MRR based FSO system, the parameter
$\sigma_{\theta_{e}}$ has a greater impact on performance of MRR-based FSO
system with respect to $\sigma_{\theta_{o}}$. From the results of [7, 48], for
a short link UAV-based FSO system, the tolerable UAVs orientation fluctuations
is in the order of several mrads and for a long link, is in the order of
$\micro$rad. However, for the small UAVs with low power and payload
limitations (that it is not possible to employ stabilizers) it is not possible
to achieve such a UAV’s angular stability. One of the main advantages of
employing MRR is to compensate of UAV’s orientation fluctuations, especially,
for the small UAVs. From the results of Fig. 11, the considered MRR-based
system has an acceptable performance even for high UAV’s angular fluctuations
equal to $\sigma_{\theta_{o}}=6^{o}=104.7$ mrad. However, due to the limited
data rate, we are forced to use a smaller $A_{r}$ that makes this MRR-based
system more sensitive to tracking system errors with respect to the
conventional FSO systems. The results of Fig. 11 clearly confirm the
aforementioned points.

Figure 12:  Average BER versus $P_{t}$ for $\theta_{div}=0.4$ mrad,
$A_{r}=1\,\textrm{cm}^{2}$, and different link lengths $Z=800$, $1000$,
$1200$, and $1400$.

Now, the impact of link length is evaluated on the performance of the
considered MRR-based FSO system in Figs. 12 and 13. In Fig. 12, the BER
performance is evaluated for different values of link lengths $Z=800$, $1000$,
$1200$, and $1400$. As expected, by increasing link length, the performance
degrades, significantly. However, it should be noted that the beamwidth at the
MRR changes by varying link length which changes the distribution of RV
$h_{pu}$ and end-to-end SNR. Accordingly, for any given link length, we must
find an optimal beamwidth to achieve minimum outage probability and/or BER.
The beamwidth is tuned by divergence angle at the GS node. Thus, for any given
link length, in the considered MRR-based FSO system, finding and tuning an
optimal value for divergence angle is very important. To get a better insight,
in Fig. 13, the outage probability is plotted versus $\theta_{div}$ for
different values of link length. The results of Fig. 13 clearly shows by
varying link length, the optimal value for $\theta_{div}$ changes and confirm
the importance of finding optimal value for $\theta_{div}$ when the link
length is varied.

In addition to the link length, any changing in other parameters such as
$\sigma_{\theta_{e}}$ can change the optimal value for beamwidth. In Fig. 14,
outage probability is depicted versus $\sigma_{\theta_{e}}$ and $w_{z}$. The
results of Fig. 14 shows that any increase in the SD of tracking system errors
causes an increase in the optimal value for beamwidth. This can be justified
since by increasing $\sigma_{\theta_{e}}$, the beamwidth must be increased to
compensate the fluctuations of Gaussian beam footprint at the MRR. However,
any increase in beamwidth increases the geometrical loss. Therefore, for any
given $\sigma_{\theta_{e}}$, the optimal value for $w_{z}$ can be obtained by
a trade-off between the strength of Gaussian beam footprint fluctuations and
geometrical loss.

Figure 13:  Outage probability versus $\theta_{div}$ for $P_{t}=20$ dBm,
$\sigma_{\theta_{o}}=5^{o}$, $\sigma_{\theta_{e}}=100$ $\micro$rad, and
different values of link length.  Figure 14:  Outage probability versus
$\sigma_{\theta_{e}}$ and $w_{z}$ for $P_{t}=25$ dBm.

##  V Conclusion and Future Road Map

In this paper, we have studied the performance of UAV-based FSO link when UAV
is equipped with MRR. Accordingly, we have characterized the MRR-based UAV FSO
channel by taking into account tracking system errors along with UAV’s
orientation fluctuations, link length, UAV’s height, optical beam divergence
angle, effective area of MRR, atmospheric turbulence and optical channel loss
in the double-pass channels. To enable effective performance analysis, we have
derived the tractable and closed-form expressions for PDF of end-to-end SNR,
outage probability and BER of the considered system under both weak-to-
moderate and moderate-to-strong atmospheric turbulence conditions. We have
then verified the accuracy of analytical models by employing Monte Carlo
simulations. Our results reveal that any change in the parameters (such as
link length, SD of tracking system errors, target BER, desired data rate,
etc.) affect the optimal values of other parameters. For MRR-based FSO
deployments, the proposed analytical methods will assist researchers to easily
analyze and design of such systems without performing any time-consuming
simulations.

Thanks to their high internal gain, avalanche photo-detector (APD) can improve
SNR capability, as compared with PIN-based receivers. However, in such APD-
based receivers, shot noise is mostly dominant. The variance of shot noise
depends on the received optical signal intensity and thus, system analysis
becomes more complex and can be considered as a future work.

Figure 15: Comparison of the accuracy of the provided PDF for the RV $h_{MRR}$
in (89) with actual distribution under under different SD of UAV orientations.

##  Appendix A The Proof of Theorem 1

We reformulate (11) as

$\displaystyle h=h_{c}h_{L}h_{p_{u}},$ (88)

where $h_{c}=h_{l_{ug}}h_{l_{gu}}h_{p_{g}}$ is a constant coefficient and
$h_{L}=h_{a_{ug}}h_{a_{gu}}h_{MRR}$. After an exhaustive search over the
distribution of $h_{MRR}$ denoted by $f_{h_{MRR}}(h_{MRR})$, we obtain an
approximate mathematical function for $f_{h_{MRR}}(h_{MRR})$ as

$\displaystyle
f_{h_{MRR}}(h_{MRR})\simeq\frac{1}{h_{MRR}\sqrt{2\pi\ln\left(1+\frac{\sigma_{MRR}^{2}}{\mu_{MRR}^{2}}\right)}}$
(89)
$\displaystyle~{}~{}~{}\times\exp\left(-\frac{\left[\ln\left(\frac{\sqrt{\mu_{MRR}^{2}+\sigma_{MRR}^{2}}}{\mu_{MRR}^{2}}h_{MRR}\right)\right]^{2}}{2\ln\left(1+\frac{\sigma_{MRR}^{2}}{\mu_{MRR}^{2}}\right)}\right).$

To find a better view about the accuracy of the approximated PDF in 32, the
accuracy of the approximated PDF for $h_{MRR}$ is compared with actual
distribution in Fig. 15 under different SD of UAV’s orientations. As shown in
Fig. 15, the results obtained from (89) is very close to the actual values of
the PDF of $h_{MRR}$ for $\sigma_{\theta_{o}}=1^{o}$. As $\sigma_{\theta_{o}}$
increases, it is observed that the results of (89) deviate slightly from the
actual values, and for $\sigma_{\theta_{o}}=10^{o}$ this deviation increases.
As will be shown in the simulation section, for the interval
$\sigma_{\theta_{o}}<8^{o}$, the end-to-end channel distribution function is
well close to the value obtained from the Monte-Carlo simulations, and for
interval $\sigma_{\theta_{o}}>8^{o}$, we will see an offset between the
analytical and the Monte-Carlo results. This offset is mainly caused by (89).
From (8) and (89), after some derivations, we obtain

$\displaystyle
f_{h_{L}}(h_{L})=\frac{1}{h_{L}\sqrt{2\pi\left[\ln\left(1+\frac{\sigma_{MRR}^{2}}{\mu_{MRR}^{2}}\right)+8\sigma_{L}^{2}\right]}}$
(90)
$\displaystyle~{}~{}~{}\times\exp\left(-\frac{\left[\ln\left(\frac{\sqrt{\mu_{MRR}^{2}+\sigma_{MRR}^{2}}}{\mu_{MRR}^{2}}h_{L}\right)+4\sigma_{L}^{2}\right]^{2}}{2\left[\ln\left(1+\frac{\sigma_{MRR}^{2}}{\mu_{MRR}^{2}}\right)+8\sigma_{L}^{2}\right]}\right).$

As shown in simulation results, for a reliable communication, the SD of
tracking system errors characterized by $\theta_{ex}$ small (less than a few
mrad). Under such tracking errors, (2) can be simplified as

$\displaystyle d_{px}\simeq
Z\theta_{ex},~{}~{}~{}~{}\&~{}~{}~{}~{}d_{py}\simeq Z\theta_{ey}.$ (91)

From (91), the distribution of $d_{p}=\sqrt{d_{px}^{2}+d_{py}^{2}}$ well
approximated as

$\displaystyle
f_{d_{p}}(d_{p})=\frac{d_{p}}{Z^{2}\sigma_{\theta_{e}}^{2}}\exp\left(-\frac{d_{p}^{2}}{2Z^{2}\sigma_{\theta_{e}}^{2}}\right),~{}~{}~{}~{}~{}~{}d_{p}\geq
0.$ (92)

From (5) and (92), we obtain

$\displaystyle f_{h_{p_{u}}}(h_{p_{u}})=K\left(\frac{\pi
w_{z}^{2}}{2A_{r}}\right)^{K}h_{p_{u}}^{K-1},~{}~{}~{}0<h_{p_{u}}\leq\frac{2A_{r}}{\pi
w_{z}^{2}},$ (93)

where $K=\frac{w_{z}^{2}}{Z^{2}\sigma_{\theta_{e}}^{2}}$. From (88), (90), and
(93), we have

$\displaystyle f_{h}(h)=\frac{K(2A_{r}h_{c})^{-K}(\pi
w_{z}^{2})^{K}h^{K-1}}{\sqrt{2\pi\left[\ln\left(1+\frac{\sigma_{MRR}^{2}}{\mu_{MRR}^{2}}\right)+8\sigma_{L}^{2}\right]}}\int_{\frac{\pi
w_{z}^{2}h}{2A_{r}h_{c}}}^{\infty}h_{L}^{-K-1}$
$\displaystyle\times\exp\left(-\frac{\left[\ln\left(\frac{\sqrt{\mu_{MRR}^{2}+\sigma_{MRR}^{2}}}{\mu_{MRR}^{2}}h_{L}\right)+4\sigma_{L}^{2}\right]^{2}}{2\left[\ln\left(1+\frac{\sigma_{MRR}^{2}}{\mu_{MRR}^{2}}\right)+8\sigma_{L}^{2}\right]}\right)dh_{L}.$
(94)

Finally, applying a change of variable
$x=\ln\left(\frac{\sqrt{\mu_{MRR}^{2}+\sigma_{MRR}^{2}}}{\mu_{MRR}^{2}}h_{L}\right)$
and using [49, eq. (2.33)], the closed-form channel distribution is derived in
(15).

##  Appendix B The Proof of Proposition 3

Let $p_{0}$ and $p_{1}$ denote the a priori probability of transmission bits
“0” and “1”, respectively. The BER of intensity modulated direct detection
with on-off keying (OOK) signaling is given by
$\mathbb{P}_{e}=p_{0}\mathbb{P}_{e|0}+p_{1}\mathbb{P}_{e|1}$ where
$\mathbb{P}_{e|0}$ and $\mathbb{P}_{e|1}$ denote the conditional bit error
probabilities when the transmitted bit is “0” and “1”, respectively.
Considering also that $p_{0}=p_{1}$ and $\mathbb{P}_{e|0}=\mathbb{P}_{e|1}$,
the BER is derived as [50]

$\displaystyle\mathbb{P}_{e}=\int_{0}^{\infty}Q\left(\sqrt{\Upsilon}\right)f_{\Upsilon}(\Upsilon)~{}\textrm{d}\Upsilon.$
(95)

Substituting (15) in (95) and using a series expansion [47, eq.
(06.25.06.0002.01)], $\mathbb{P}_{e}$ is well approximated as

$\displaystyle\mathbb{P}_{e}$
$\displaystyle\simeq\frac{C_{4}}{4\Upsilon_{1}^{K/2}}\left(\int_{0}^{\Upsilon_{\textrm{max}}}\Upsilon^{\frac{K-2}{2}}Q\left(\frac{\ln(\Upsilon)-\ln(\Upsilon_{1})+2C_{5}}{2\sqrt{C_{1}}}\right)\textrm{d}\Upsilon\right.$
$\displaystyle~{}~{}~{}-\frac{2}{\sqrt{\pi}}\sum_{m=0}^{M}\frac{(-1)^{m}}{m!(2m+1)2^{m+\frac{1}{2}}}\int_{0}^{\Upsilon_{\textrm{max}}}\Upsilon^{\frac{2m+K-1}{2}}$
$\displaystyle\left.{\color[rgb]{1,1,1}\int_{0}^{\Upsilon_{\textrm{max}}}}\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\times
Q\left(\frac{\ln(\Upsilon)-\ln(\Upsilon_{1})+2C_{5}}{2\sqrt{C_{1}}}\right)\textrm{d}\Upsilon\right).$
(96)

where $M=20$ and $\Upsilon_{\textrm{max}}=4$. In the following derivation, we
use an integral identity [47, eq. (06.27.21.0011.01)]

$\displaystyle\int
e^{bx}\textrm{erfc}(ax)=\frac{1}{b}\left(e^{bx}\textrm{erfc}(ax)-e^{\frac{b^{2}}{2a^{2}}}\textrm{erf}\left(\frac{b}{2a}-ax\right)\right),$
(97)

where $\textrm{erf}(\cdot)$ is the error function and $\textrm{erfc}(\cdot)$
is the complementary error function [31]. Applying a change of variables
$x=[\ln(\Upsilon)-\ln(\Upsilon_{1})+2C_{5}]$, and given the fact that
$Q(x)=\frac{1}{2}\textrm{erfc}\left(\frac{x}{\sqrt{2}}\right)$ and
$Q(x)=\frac{1}{2}-\frac{1}{2}\textrm{erf}\left(\frac{x}{\sqrt{2}}\right)$, the
closed form expression for (B) is derived in (III-A).

##  Appendix C The Proof of Theorem 2

We reformulate (11) as

$\displaystyle
h=\underbrace{h_{l_{ug}}h_{l_{gu}}h_{p_{g}}}_{h_{c}}\underbrace{h_{a_{gu}}\overbrace{h_{a_{ug}}h_{p_{u}}}^{h^{\prime}}}_{h^{\prime\prime}}h_{MRR}.$
(98)

In the sequel, we have used the PDF of product of two RVs, which is generally
given as [51]

$\displaystyle
f_{z}(z)=\int_{-\infty}^{\infty}\frac{1}{|x|}f_{x}(x)f_{y}\left(\frac{z}{x}\right)\textrm{d}x,$
(99)

where $z$ is the product of RVs $x$ and $y$ ($z=xy$). We express the
$k_{\nu}(.)$ in terms of the Meijer’s G-function as
$k_{\nu}(x)=G^{2,0}_{0,2}\left(\frac{x^{2}}{4}\bigg{|}\begin{array}[]{c}-\\\
(\nu/2),-(\nu/2)\end{array}\right)$ [31]. Substituting (10) and (93) in (99)
and using [47, eqs. (07.34.21.0002.01) and (07.34.17.0007.01)], the PDF of RV
$h^{\prime}$ is obtained as

$\displaystyle f_{h^{\prime}}(h^{\prime})$ $\displaystyle=\frac{\pi
w_{z}^{2}\alpha\beta K}{2A_{r}\Gamma(\alpha)\Gamma(\beta)}$ (100)
$\displaystyle~{}~{}~{}\times G^{3,0}_{1,3}\left(\frac{\pi
w_{z}^{2}\alpha\beta h^{\prime}}{2A_{r}}\bigg{|}\begin{array}[]{c}K\\\
K-1,\alpha-1,\beta-1\end{array}\right).$ (103)

Now, substituting (10) and (100) in (99) and using [52, eq. (21)], the PDFof
RV $h^{\prime\prime}=h_{a_{gu}}h^{\prime}$ is derived in (106).

$\displaystyle f_{h^{\prime\prime}}(h^{\prime\prime})=\frac{K\left(\pi
w_{z}^{2}\alpha^{2}\beta^{2}\right)^{\frac{\alpha+\beta}{2}}}{(\Gamma(\alpha)\Gamma(\beta))^{2}\left(2A_{r}\right)^{\frac{\alpha+\beta}{2}}}h^{\prime\prime\frac{\alpha+\beta}{2}-1}G^{0,5}_{5,1}\left(\frac{2A_{r}}{\pi
w_{z}^{2}\alpha^{2}\beta^{2}h^{\prime\prime}}\bigg{|}\begin{array}[]{c}\frac{2-\alpha+\beta}{2},\frac{2-\beta+\alpha}{2},\frac{2+\alpha+\beta}{2}-K,\frac{2+\beta-\alpha}{2},\frac{2+\alpha-\beta}{2}\\\
\frac{\alpha+\beta}{2}-K\end{array}\right).$ (106)

Using [47, eqs. (07.34.16.0001.01) and (07.34.16.0002.01)] the (106) can be
further simplified as (107).

$\displaystyle f_{h^{\prime\prime}}(h^{\prime\prime})=\frac{\pi
w_{z}^{2}\alpha^{2}\beta^{2}K}{2A_{r}(\Gamma(\alpha)\Gamma(\beta))^{2}}\times$
(107) $\displaystyle G^{5,0}_{1,5}\left(\frac{\pi
w_{z}^{2}\alpha^{2}\beta^{2}h^{\prime\prime}}{2A_{r}}\bigg{|}\begin{array}[]{c}K\\\
\alpha-1,\beta-1,K-1,\alpha-1,\beta-1\end{array}\right).$ (110)

We propose an approximate sectorized model for $f_{h_{MRR}}(h_{MRR})$ as

$\displaystyle f_{h_{MRR}}$ $\displaystyle(h_{MRR})\simeq$ (111)
$\displaystyle\sum_{n=1}^{N}B_{n}\left[U(h_{MRR}-V_{n})-U(h_{MRR}-V_{n+1})\right],$

where

$V_{n}=\left\\{\begin{array}[]{ll}&\\!\\!\\!\\!\\!\\!\\!2\mu_{MRR}-1,~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}n=1,\\\
&\\!\\!\\!\\!\\!\\!\\!V_{n-1}+\frac{(2-2\mu_{MRR})n}{N},~{}~{}~{}n\in\\{2,...,N+1\\},\end{array}\right.$

and $U(x)=\left\\{\begin{array}[]{ll}&\\!\\!\\!\\!\\!\\!\\!1,~{}~{}~{}x>0,\\\
&\\!\\!\\!\\!\\!\\!\\!0,~{}~{}~{}x<0,\end{array}\right.$ is the Heaviside step
function. Also, the coefficients $B_{n}$ in (111) depend on the parameters
$\sigma_{\theta_{o}}$ and the number of sectors denoted by $N$. In Fig. 16, we
compare the proposed sectorized model with respect to the distribution of
$h_{MRR}$ obtained using simulation for two different values of
$\sigma_{\theta_{o}}$ and $N=12$.

Figure 16:  Sectorized model of the PDF of $h_{MRR}$ for $N=12$ and comparison
with the distribution of $h_{MRR}$ obtained using simulation for two different
values of $\sigma_{\theta_{o}}$.

Obviously, the accuracy of the proposed sectorized model directly depends on
the number of sectors $N$, and for sufficiently large values of $N$, an exact
match between simulations and analysis can be achieved at the cost of higher
complexity. Next, we use the proposed model to find a closed-form analytical
model for $h$ under moderate to strong turbulence conditions. In the section
of simulation results, it is shown that $N=8$ achieves a sufficient accuracy.
In Table III, the coefficients $B_{n}$ obtained using simulation for different
values of $\sigma_{\theta_{o}}$ and $N=8$. For the rest of
$\sigma_{\theta_{o}}$, the related coefficients $B_{n}$ can be obtained by
interpolation from the given values. Substituting (107) and (111) in (99), we
obtain (C).

$\displaystyle f_{h}(h)=\frac{K\left(\pi
w_{z}^{2}\alpha^{2}\beta^{2}\right)^{\frac{\alpha+\beta}{2}}}{(\Gamma(\alpha)\Gamma(\beta))^{2}\left(2A_{r}h_{c}\right)^{\frac{\alpha+\beta}{2}}}h^{\frac{\alpha+\beta}{2}-1}\sum_{n=1}^{N}B_{n}$
$\displaystyle\left[\int_{0}^{V_{n+1}}h_{MRR}^{-\frac{\alpha+\beta}{2}}G^{0,5}_{5,1}\left(\frac{2A_{r}h_{c}}{\pi
w_{z}^{2}\alpha^{2}\beta^{2}h}h_{MRR}\bigg{|}\begin{array}[]{c}\frac{2-\alpha+\beta}{2},\frac{2-\beta+\alpha}{2},\frac{2+\alpha+\beta}{2}-K,\frac{2+\beta-\alpha}{2},\frac{2+\alpha-\beta}{2}\\\
\frac{\alpha+\beta}{2}-K\end{array}\right)\textrm{d}h_{MRR}\right.$ (114)
$\displaystyle-\left.\int_{0}^{V_{n}}h_{MRR}^{-\frac{\alpha+\beta}{2}}G^{0,5}_{5,1}\left(\frac{2A_{r}h_{c}}{\pi
w_{z}^{2}\alpha^{2}\beta^{2}h}h_{MRR}\bigg{|}\begin{array}[]{c}\frac{2-\alpha+\beta}{2},\frac{2-\beta+\alpha}{2},\frac{2+\alpha+\beta}{2}-K,\frac{2+\beta-\alpha}{2},\frac{2+\alpha-\beta}{2}\\\
\frac{\alpha+\beta}{2}-K\end{array}\right)\textrm{d}h_{MRR}\right].$ (117)

Finally using [47, 07.34.21.0084.01] and after some manipulations, the closed-
form expression for $h$ is derived in (36). Substituting (36) in (24) and
using [47, eq. (07.34.21.0084.01)], the CDF of $h$ is derived in (64). Based
on (28), we obtain

$\displaystyle
F_{\Upsilon}(\upsilon)=F_{h}\left(\sqrt{\frac{\Upsilon}{\Upsilon_{1}}}\right),$
(118)

and then substitute (43) into (118), the CDF of $\Upsilon$ can be derived in
(57).

Furthermore, from (111), for lower values of $\sigma_{\theta_{o}}$, the (36)
and (50) can be simplified as (64) and (72), respectively. Now, substituting
(64) and (72) in (24) and using [47, eq. (07.34.21.0084.01)], the CDF of $h$
and $\Upsilon$ are approximated as (68) and (76), respectively.

Also, for small values of $\sigma_{t_{o}}$ (lower than $0.05^{o}$),
$f_{h_{MRR}}=({h_{MRR}})$ can be well approximated by a Dirac delta function
as $f_{h_{MRR}}({h_{MRR}})\simeq\delta(h_{MRR}-1)$. From this, and by
following the method used for obtaining (36) and (43), the PDF and CDF of $h$
can be approximated as (64) and (68), respectively.

##  Appendix D The Proof of Proposition 6

Lets to rewrite the Gaussian Q-function as complementary error function by
$Q(x)=2\text{erfc}\left(\sqrt{2}x\right)$. Also, using [47, eq.
(06.27.26.0006.01)], we can rewrite $Q(x)$ as

$\displaystyle
Q(x)=\frac{1}{2\sqrt{\pi}}G^{2,0}_{1,2}\left(\frac{x^{2}}{2}\bigg{|}\begin{array}[]{c}1\\\
0,1/2\end{array}\right).$ (121)

Using (50) and (95), the BER of the considered system can be obtained as

$\displaystyle\mathbb{P}_{e}=\frac{B_{s}}{4\sqrt{\pi\Upsilon_{1}\Upsilon}}\sum_{n=1}^{N}B_{n}\int_{0}^{\infty}G^{2,0}_{1,2}\left(\frac{\Upsilon}{2}\bigg{|}\begin{array}[]{c}1\\\
0,1/2\end{array}\right)$ (124)
$\displaystyle\left[G^{6,0}_{2,6}\left(B_{n}^{\prime}\sqrt{\frac{\Upsilon}{\Upsilon_{1}}}\bigg{|}\begin{array}[]{c}K,1\\\
0,\alpha-1,\beta-1,K-1,\alpha-1,\beta-1\end{array}\right)-\right.$ (127)
$\displaystyle\left.G^{6,0}_{2,6}\left(B_{n}^{\prime\prime}\sqrt{\frac{\Upsilon}{\Upsilon_{1}}}\bigg{|}\begin{array}[]{c}K,1\\\
0,\alpha\\!-\\!1,\beta\\!-\\!1,K-1,\alpha-1,\beta-1\end{array}\right)\right]~{}\textrm{d}\Upsilon.$
(130)

Finally, using (124) and [52, eq. (21)], after some manipulations, the closed-
form expressions for BER is derived in (82).

## References

  * [1] M. Mozaffari, W. Saad, M. Bennis, Y.-H. Nam, and M. Debbah, “A tutorial on UAVs for wireless networks: Applications, challenges, and open problems,” _IEEE Commun. Surveys Tuts._ , vol. 21, no. 3, pp. 2334–2360, 2019.
  * [2] M. T. Dabiri, H. Safi, S. Parsaeefard, and W. Saad, “Analytical channel models for millimeter wave UAV networks under hovering fluctuations,” _IEEE Trans. Wireless Commun._ , vol. 19, no. 4, pp. 2868–2883, 2020.
  * [3] M. Alzenad, M. Z. Shakir, H. Yanikomeroglu, and M.-S. Alouini, “FSO-based vertical backhaul/fronthaul framework for 5G+ wireless networks,” _IEEE Commun. Mag._ , vol. 56, no. 1, pp. 218–224, 2018.
  * [4] “Mynaric: The technology behind laser communication.” [Online]. Available: https://mynaric.com/technology/overview
  * [5] M. A. Khalighi and M. Uysal, “Survey on free space optical communication: A communication theory perspective,” _IEEE Commun. Surveys Tuts_ , vol. 16, no. 4, pp. 2231–2258, 2014.
  * [6] H. Kaushal and G. Kaddoum, “Optical communication in space: Challenges and mitigation techniques,” _IEEE Commun. Surveys Tuts._ , vol. 19, no. 1, pp. 57–96, 2016.
  * [7] M. T. Dabiri and S. M. S. Sadough, “Optimal placement of UAV-assisted free-space optical communication systems with DF relaying,” _IEEE Comm. Let._ , vol. 24, no. 1, pp. 155–158, 2019.
  * [8] M. T. Dabiri, S. Khankalantary, M. J. Piran, I. S. Ansari, M. Uysal, W. Saad, and C. S. Hong, “UAV-assisted free space optical communication system with amplify-and-forward relaying,” _IEEE Trans. Veh. Technol._ , vol. 70, no. 9, pp. 8926–8936, 2021.
  * [9] W. Fawaz, C. Abou-Rjeily, and C. Assi, “UAV-aided cooperation for FSO communication systems,” _IEEE Commun. Mag._ , vol. 56, no. 1, pp. 70–75, 2018.
  * [10] S. Zhang and N. Ansari, “3D drone base station placement and resource allocation with FSO-based backhaul in hotspots,” _IEEE Trans. Veh. Technol._ , vol. 69, no. 3, pp. 3322–3329, 2020.
  * [11] Z. Gu, J. Zhang, and Y. Ji, “Topology Optimizing in FSO-based UAVs Relay Networks for Resilience Enhancement,” _Mobile Networks and Applications_ , vol. 25, no. 1, pp. 350–362, 2020.
  * [12] Z. Gu, J. Zhang, X. Sun, and Y. Ji, “Optimizing networked flying platform deployment and access point association in FSO-based fronthaul networks,” _IEEE Wireless Commun. Let._ , 2020.
  * [13] J.-H. Lee, K.-H. Park, Y.-C. Ko, and M.-S. Alouini, “A UAV-mounted free space optical communication: Trajectory optimization for flight time,” _IEEE Trans. Wireless Commun._ , vol. 19, no. 3, pp. 1610–1621, 2019.
  * [14] H. Ajam, M. Najafi, V. Jamali, and R. Schober, “Ergodic sum rate analysis of UAV-based relay networks with mixed RF-FSO channels,” _IEEE Open J. Commun. Soc._ , vol. 1, pp. 164–178, 2020.
  * [15] M. Najafi, H. Ajam, V. Jamali, P. D. Diamantoulakis, G. K. Karagiannidis, and R. Schober, “Statistical Modeling of the FSO Fronthaul Channel for UAV-based Communications,” _IEEE Trans. Commun._ , 2020.
  * [16] M. T. Dabiri, S. M. S. Sadough, and M. A. Khalighi, “Channel modeling and parameter optimization for hovering UAV-based free-space optical links,” _IEEE J. Sel. Areas Commun._ , vol. 36, no. 9, pp. 2104–2113, 2018.
  * [17] M. T. Dabiri, S. M. S. Sadough, and I. S. Ansari, “Tractable optical channel modeling between UAVs,” _IEEE Trans. Veh. Technol._ , vol. 68, no. 12, pp. 11 543–11 550, 2019.
  * [18] M. T. Dabiri, M. Rezaee, I. S. Ansari, and V. Yazdanian, “Channel Modeling for UAV-based Optical Wireless Links with Nonzero Boresight Pointing Errors,” _IEEE Trans. Veh. Technol._ , 2020.
  * [19] M. T. Dabiri, M. Rezaee, V. Yazdanian, B. Maham, W. Saad, and C. S. Hong, “3D channel characterization and performance analysis of UAV-assisted millimeter wave links,” _IEEE Trans. Wireless Commun._ , 2020.
  * [20] H. Safi, A. Dargahi, J. Cheng, and M. Safari, “Analytical Channel Model and Link Design Optimization for Ground-to-HAP Free-Space Optical Communications,” _J. Lightw. Technol._ , 2020.
  * [21] P. G. Goetz, W. S. Rabinovich, R. Mahon, J. L. Murphy, M. S. Ferraro, M. R. Suite, W. R. Smith, H. R. Burris, C. I. Moore, W. W. Schultz _et al._ , “Modulating retro-reflector lasercom systems for small unmanned vehicles,” _IEEE J. Sel. Areas Commun._ , vol. 30, no. 5, pp. 986–992, 2012.
  * [22] H. Kaushal and G. Kaddoum, “Underwater optical wireless communication,” _IEEE access_ , vol. 4, pp. 1518–1547, 2016.
  * [23] C. Quintana, G. Erry, A. Gomez, Y. Thueux, G. Faulkner, and D. O’Brien, “Design of a holographic tracking module for long-range retroreflector free-space systems,” _Applied Optics_ , vol. 55, no. 25, pp. 7173–7178, 2016\.
  * [24] C. Quintana, A. Gomez, G. Faulkner, G. R. Erry, Y. Thueux, and D. O’Brien, “Novel non-mechanical fine tracking module for retroreflective free space optics,” in _Unmanned/Unattended Sensors and Sensor Networks X_ , vol. 9248\. International Society for Optics and Photonics, 2014, p. 92480S.
  * [25] G. Yang, S. You, M. Bi, B. Fan, Y. Lu, X. Zhou, J. Li, H. Geng, and T. Wang, “Wave-optics simulation of the double-pass beam propagation in modulating retro-reflector FSO systems using a corner cube reflector,” _Applied optics_ , vol. 56, no. 26, pp. 7474–7483, 2017.
  * [26] G. Yang, Z. Li, M. Bi, X. Zhou, R. Zeng, T. Wang, and J. Li, “Channel modeling and performance analysis of modulating retroreflector FSO systems under weak turbulence conditions,” _IEEE Photon. J._ , vol. 9, no. 2, pp. 1–10, 2017.
  * [27] G. Yang, C. Li, J. Li, H. Geng, M. Bi, B. Fan, and T. Wang, “Performance analysis of full duplex modulating retro-reflector free-space optical communications over single and double gamma-gamma fading channels,” _IEEE Trans. Commun._ , vol. 66, no. 8, pp. 3597–3609, 2018.
  * [28] G. Yang, J. Zhang, J. Zhang, M. Bi, T. Chen, S. You, X. Zhou, T. Wang, J. Li, and H. Geng, “Wavefront compensation with the micro corner-cube reflector array in modulating retroreflector free-space optical channels,” _J. Lightw. Technol._ , 2020.
  * [29] X. Li, X. Zhao, P. Zhang, W. Yang, T. Wang, and H. Jiang, “Probability density function of turbulence fading in MRR free space optical link and its applications in MRR free space optical communications,” _IET Commun._ , vol. 11, no. 16, pp. 2476–2481, 2017.
  * [30] X. Li, X. Zhao, and P. Zhang, “Bit error rate analysis for modulating retro-reflector free space optical communications with adaptive threshold over correlated gamma gamma fading channels,” _IEEE Commun. Let._ , vol. 23, no. 12, pp. 2275–2278, 2019.
  * [31] I. S. Gradshteyn and I. M. Ryzhik, _Table of integrals, series, and products_. 7th ed. Academic press, 2007\.
  * [32] E. W. Weisstein, “Meijer G-Function.” [Online]. Available: https://mathworld.wolfram.com/MeijerG-Function.html
  * [33] B. E. Saleh and M. C. Teich, _Fundamentals of photonics_. john Wiley & sons, 2019.
  * [34] Z. Ghassemlooy, W. Popoola, and S. Rajbhandari, _Optical wireless communications: system and channel modelling with Matlab®_. CRC press, 2019\.
  * [35] S. Huang and M. Safari, “Free-space optical communication impaired by angular fluctuations,” _IEEE Trans. Wireless Commun._ , vol. 16, no. 11, pp. 7475–7487, 2017.
  * [36] L. C. Andrews and R. L. Phillips, “Laser beam propagation through random media.” SPIE, 2005.
  * [37] K. Kiasaleh, “Beam-tracking in FSO links impaired by correlated fading,” in _Free-Space Laser Communications VI_ , vol. 6304. International Society for Optics and Photonics, 2006, p. 63041I.
  * [38] J. Yanmei, L. Congmin, S. Pengfei, and L. Lu, “Modulated retro-reflector-based physical-layer network coding for space optical communications,” _IEEE Access_ , vol. 9, pp. 44 868–44 880, 2021.
  * [39] W. Nelson, J. Palastro, C. Wu, and C. Davis, “Enhanced backscatter of optical beams reflected in turbulent air,” _JOSA A_ , vol. 32, no. 7, pp. 1371–1378, 2015.
  * [40] S. Althunibat, Z. Altarawneh, and R. Mesleh, “Performance analysis of free space optical–based wireless sensor networks using corner cube retroreflectors,” _Transactions on Emerging Telecommunications Technologies_ , vol. 30, no. 12, p. e3707, 2019.
  * [41] B. M. El Saghir, M. B. El Mashade, and A. M. Aboshosha, “Performance analysis of modulating retro-reflector FSO communication systems over Málaga turbulence channels,” _Optics Communications_ , vol. 474, p. 126160, 2020\.
  * [42] D. Geng, P. Du, W. Wang, G. Gao, T. Wang, and M. Gong, “Single laser free-space duplex communication system with adaptive threshold technique and BER analysis in weak turbulent atmosphere,” _Optics letters_ , vol. 39, no. 13, pp. 3950–3953, 2014.
  * [43] A. K. Majumdar, “Modulating retroreflector-based free-space optical (FSO) communications,” in _Advanced free space optics (FSO)_. Springer, 2015, pp. 243–293.
  * [44] A. M. Scott and K. D. Ridley, “Calculations of bit error rates for retroreflective laser communications systems in the presence of atmospheric turbulence,” in _Advanced Free-Space Optical Communications Techniques and Technologies_ , vol. 5614. International Society for Optics and Photonics, 2004, pp. 31–42.
  * [45] G. G. Peter, S. R. William, R. Mahon, L. M. James, S. F. Mike, R. S. Michele, R. S. Walter, B. X. Ben, R. B. Harris, I. M. Christopher, W. S. Warren, T. F. Wade, S. Frawley, M. Colbert, M. M. Barry, and K. H. S. Reese, “Modulating retro-reflector lasercom systems at the Naval Research Laboratory,” in _2010 - MILCOM 2010 MILITARY COMMUNICATIONS CONFERENCE_ , 2010, pp. 1601–1606.
  * [46] M. L. Plett, W. S. Rabinovich, R. Mahon, M. S. Ferraro, P. G. Goetz, C. I. Moore, and W. Freeman, “Free-space optical communication link across,” _Optical Engineering_ , vol. 47, no. 4, p. 045001, 2008.
  * [47] Wolfram, “The wolfram functions site: http://functions.wolfram.com/,” 2001.
  * [48] S. Khankalantary, M. T. Dabiri, and H. Safi, “BER performance analysis of drone-assisted optical wireless systems with APD receiver,” _Optics Communications_ , vol. 463, p. 125309, 2020.
  * [49] G. Izrail’S and I. M. Ryžik, _Table of integrals, series, and products_. Acad. Press, 1996.
  * [50] M. K. Simon and M.-S. Alouini, _Digital communication over fading channels_. John Wiley & Sons, 2005, vol. 95.
  * [51] A. Papoulis and S. U. Pillai, _Probability, random variables, and stochastic processes_. Tata McGraw-Hill Education, 2002.
  * [52] V. Adamchik and O. Marichev, “The algorithm for calculating integrals of hypergeometric type functions and its realization in reduce system,” in _Proceedings of the international symposium on Symbolic and algebraic computation_ , 1990, pp. 212–224.