Ergodic Sum Rate Analysis of UAV-Based Relay Networks with Mixed RF-FSO Channels

Unmanned aerial vehicle (UAV)-based communications is a promising new technology that can add a wide range of new capabilities to the current network infrastructure. Given the flexibility, cost-efficiency, and convenient use of UAVs, they can be deployed as temporary base stations (BSs) for on-demand situations like BS overloading or natural disasters. In this work, a UAV-based communication system with radio frequency (RF) access links to the mobile users (MUs) and a free-space optical (FSO) backhaul link to the ground station (GS) is considered. In particular, the RF and FSO channels in this network depend on the UAV's positioning and (in)stability. The relative position of the UAV with respect to the MUs impacts the likelihood of a line-of-sight (LOS) connection in the RF link and the instability of the hovering UAV affects the quality of the FSO channel. Thus, taking these effects into account, we analyze the end-to-end system performance of networks employing UAVs as buffer-aided (BA) and non-buffer-aided (non-BA) relays in terms of the ergodic sum rate. Simulation results validate the accuracy of the proposed analytical derivations and reveal the benefits of buffering for compensation of the random fluctuations caused by the UAV's instability. Our simulations also show that the ergodic sum rate of both BA and non-BA UAV-based relays can be enhanced considerably by optimizing the positioning of the UAV. We further study the impact of the MU density and the weather conditions on the end-to-end system performance.


I. INTRODUCTION
Fifth generation (5G) and beyond wireless communication networks are expected to overcome many of the shortcomings of the current infrastructure by offering higher data rates, improving the quality of service (QoS) in crowded areas, and reducing the blind spots of current networks [1]. Among other techniques, unmanned aerial vehicles (UAVs) have been introduced to achieve the aforementioned goals. The unique characteristic of UAVs is their flexible positioning which together with their cost efficiency and easy deployment makes them promising candidates for a wide range of applications. For example, they may be used as relays arXiv:2001.00871v2 [eess.SP] 10 Jan 2020 for coverage enhancement [2], as temporary base stations (BSs) for on-demand situations [3], for adaptive fronthauling/backhauling [4], and for data acquisition for the Internet-of-Things (IoT) [5]. Despite their expected benefits, the backhaul/fronthaul links needed to connect a UAV with a ground station (GS) constitute a major challenge in UAV-based networks. The authors of [6] considered WiFi and satellite links for backhauling. However, for many applications, UAVs have to transfer huge amounts of data to the GS and the backhaul links have to be able to cope with the UAV's mobility and the interference from other UAVs and the mobile users (MUs). To address these issues, the authors of [4] and [7] proposed free-space optical (FSO) systems for the fronthaul/backhaul connections in UAV-based networks. FSO links offer high data rates (up to 10 Gbps) by using the optical range of the frequency spectrum. Moreover, FSO systems are not suceptible to interference owing to their narrow laser beams and are able to communicate over large distances (few kilometers) [8].
Considering the above advantages, the QoS and coverage of the current infrastructure can be improved by UAVs operating as flying relays. Particularly, UAVs can be deployed as mobile relay nodes that forward the huge amounts of data collected via RF access links from MUs to a GS via FSO backhaul links. The viability of UAV-based communication systems has been demonstrated recently in [9] and [10], where end-to-end long term evolution (LTE) connectivity was provided by low altitude UAVs and balloons, respectively. Furthermore, UAVs with FSO backhauling were utilized in the Aquila [11] and Loon [12] projects for providing connectivity to the remote parts of the world. However, the performance of such UAV-based relay networks has not been studied in detail and is expected to be strongly dependent on the UAV's positioning and (in)stability. In particular, the location of the UAV with respect to (w.r.t.) the MUs and GS and its random vibrations in the hovering state affect the quality of the RF and FSO channels. Thus, the impact of the UAV's positioning and instability on the following parameters necessitates a careful study of the performance of the end-to-end network: 1) MU distribution: The MUs are randomly distributed and their data traffic patterns may change over time. The flexible positioning of the UAV above the randomly distributed MUs can reduce path loss and shadowing effects in the RF access links and boost the end-to-end system performance.
2) Line-of-sight (LOS) link: The probability of maintaining a LOS path for the RF access links depends on the elevation angle of the UAV w.r.t. the MUs. When the UAV is at higher altitudes, the LOS path between the UAV and a given MU is less likely to be blocked. Hence, depending on the position of the UAV, the distribution of the RF access channel coefficients can be either Rician in the presence of a LOS path or Rayleigh in the absence of a LOS path. Thus, the positioning of the UAV determines the LOS probability.
3) Quality of the FSO link: In UAV-based FSO communications, tracking errors and the instability of the hovering UAV degrade the intensity of the optical signal received at the photo detector (PD) of the GS. Furthermore, the distance between the UAV and the GS affects the atmospheric loss in the FSO channel.
The above factors have to be taken simultaneously into account for the design of UAV-based networks. For instance, the position of the UAV w.r.t. the MUs affects both the LOS probability of the access links and the atmospheric loss in the backhaul channel. However, in previous works, only some of the above aspects were considered. For example, the authors of [13] considered only the impact of the RF access links and assumed a perfect backhaul connection. They determined the optimal placement of a stationary UAV and the optimal trajectory of a moving UAV in terms of the maximum throughput. Furthermore, the performance of a cluster-based UAV network in terms of coverage probability and energy efficiency was analyzed in [3]. The authors of [14] investigated the optimal positioning in a multi-UAV network with the objective to minimize the total transmit power. In both [3] and [14], the backhaul link was assumed to be ideal. The authors of [15] considered the impact of the positioning of the UAV on both the access and the backhaul links. However, in [15], despite the potentially higher data rates of FSO links, an RF link was considered for backhauling and the position of the UAV and the RF bandwidth shared between the access and backhaul links was optimized. In fact, to the best of the authors' knowledge, the performance of a UAV-based network employing FSO backhaul and RF access links to connect randomly distributed MUs to a GS has not been studied in the literature, yet.
In this paper, we consider a UAV-based network where a hovering UAV acts as a decodeand-forward (DF) relay and connects Poisson distributed MUs via RF access links and an FSO backhaul link to a fixed GS. Because of the mutual orthogonality of the RF and FSO links, the relaying UAV can concurrently transmit and receive. Moreover, depending on whether the data is delay sensitive or not, non-buffer-aided (non-BA) and buffer-aided (BA) relaying UAVs are considered, respectively [16]- [18]. In particular, BA relaying allows the UAV to transmit, receive, or simultaneously transmit and receive depending on the channel conditions [19]. The performance of the considered UAV-based relay network is analyzed in terms of the ergodic sum rate for both BA and non-BA relaying. We validate the proposed analytical results with computer simulations. The main contributions of this paper can be summarized as follows: • End-to-end system model including both access and backhaul links: We investigate the end-to-end performance of UAV-based relay systems and show that the performance is simultaneously dependent on both the access and backhaul links. Thereby, if one link degrades the end-to-end performance, the position of the UAV can be adjusted accordingly to enhance the performance.
• Impact of UAV positioning and (in)stability: The end-to-end ergodic sum rate of the system is analyzed taking into account the impact of the positioning of the UAV w.r.t. the MUs and the GS and the random fluctuations of the position and orientation of the UAV in the hovering state which in turn affect the quality of the access and backhaul links.
We show that these characteristics of UAVs have to be jointly considered for performance analysis of UAV-based communication networks.
• Impact of buffering on performance: We analyze the ergodic sum rate for both BA and non-BA UAV-based relay systems and show that buffering can mitigate the randomness of the FSO link quality induced by the instability of the UAV. Our simulation results reveal that buffering improves the performance of the system at the expense of an increased delay.
We also show that the optimal position of the UAV is in general different for BA and non-BA UAV-based relay systems.
• Impact of weather conditions and MU density: We show that the system performance and the optimal position of the UAV strongly depend on the atmospheric conditions of the backhaul channel and the density of the MUs in the access channel. Our simulation results reveal that by optimal positioning of the UAV the impact of these system and channel parameters can be significantly reduced.
The remainder of this paper is organized as follows. In Section II, the considered system and channel models for UAV-based communication are presented. In Section III, the end-to-end performance of BA and non-BA UAV-based relaying systems is analyzed in terms of the ergodic sum rate, respectively. Simulation results are provided in Section IV, and conclusions are drawn in Section V.
Notations: In this paper, (·) T and (·) H denote the transpose and Hermitian transpose of a matrix, respectively. E{·}, * , and · denote the expectation operator, the convolution operator, and the 2 -norm of a vector, respectively. R + denotes the set of positive real numbers. represents the n × n identity matrix. x ∼ N (µ, Γ) and y ∼ CN (µ, Γ) indicate that x and y are respectively real and complex Gaussian random vectors with mean vector µ and covariance . y ∼ H(q, ω) represents a Hoyt distributed RV y with shape parameter q and spread factor ω. z ∼ lognormal(µ, σ 2 ) is used to indicate that z is a lognormal distributed RV where µ and σ 2 are the mean and variance in dB. Finally, w ∼ Nakagami(m, Ω) indicates that w is a Nakagami distributed RV with shape parameter m and spread factor Ω.

II. SYSTEM AND CHANNEL MODELS
In this section, first we present the system model for the considered UAV-based relay network facilitating uplink communication between multiple MUs and a GS. Subsequently, we introduce the channel models for the MU-to-UAV RF access links and the UAV-to-GS FSO backhaul link.

A. System Model
We consider an uplink UAV-based communication system, cf. Fig. 1, with K single-antenna MUs transmitting data over RF links to a hovering UAV equipped with N RF antennas and a single aperture which relays the received data over an FSO backhaul link to a GS equipped with a single PD. The position of the GS, MUs, and UAV are as follows. The GS is installed at height z GS of a building located in the origin of the Cartesian coordinate system (x, y, z), i.e., the GS has coordinates (0, 0, z GS ). Moreover, the K MUs are randomly distributed in a cell of radius r 0 centered at (x 0 , y 0 , 0). The random position of MU k , k ∈ {1, 2, . . . , K}, is modeled by a homogeneous Poisson point process, Φ, with density λ, and is characterized by polar coordinates (r k , ϕ k ) ∈ Φ, where 0 ≤ r k ≤ r 0 is the radial distance from the cell center and ϕ k ∈ [0, 2π] is the polar angle. Furthermore, the UAV is located at position p d = (x d , y d , z d ) and its beam direction is determined by orientation variables o d = (θ d , φ d ). Consider (x , y , z ) as a translation of the primary coordinate system, (x, y, z), by vector (x d , y d , z d ), cf. Fig. 1. Then, φ d is defined as the angle between the laser beam and the z -axis and θ d is the angle between the x -axis and the projection of the beam onto the x − y plane. Now, the uplink transmission can be divided into two hops, the MUs-to-UAV hop and the UAV-to-GS hop. The MUs are connected via RF access links to the UAV and by assuming frequency division multiple access (FDMA), each MU's signal is assigned to an orthogonal subchannel. Thus, the signal received at the UAV from MU k is given by where h k = [h k,1 , . . . , h k,n , . . . , h k,N ] T , h k,n is the flat fading channel coefficient from MU k to the n-th antenna of the UAV, and x k is the transmit symbol of MU k with power P = E{|x k | 2 }.
Assuming perfect channel state information (CSI) at the UAV, the elements of the received signal vector, y k , are combined by maximum ratio combining (MRC), and the resulting signal is decoded and forwarded to the GS via the FSO link. Furthermore, the UAV can transmit and receive simultaneously due to the mutual orthogonality of the RF and FSO channels. We consider both non-BA and BA relaying at the UAV. In the former case, the UAV immediately forwards the data received from the MUs to the GS, whereas, in the latter case, the UAV can select to either receive and transmit the packets in the same time slot or to store them in its buffer and forward them in a later time slot when the FSO channel conditions are more favorable [19]. Therefore, unlike non-BA relaying, BA relaying relaxes the constraint to transmit and receive according to a predetermined schedule and the UAV can select the best strategy (transmit, receive, or transmit and receive simultaneously) based on the conditions of the MUs-to-UAV and UAV-to-GS links.
Hence, BA relaying can enhance the end-to-end achievable rate at the expense of an increased delay [16], [20].
For the UAV-to-GS hop, the UAV maintains an FSO connection to the GS via a single laser aperture. Next, assuming perfect CSI at the GS and an intensity modulation and direct detection (IM/DD) FSO system, the signal intensity received at the GS's PD can be modeled as wherex ∈ R + is the intensity modulated optical signal, g is the scalar FSO channel coefficient, is the background Gaussian shot noise obtained after removing the ambient noise [21].

B. UAV-Based RF Channel Model
The main difference between UAV-based RF channels and the conventional channel models for satellite and terrestrial mobile communications lies in the characteristics of the LOS component.
In particular, in satellite channels, the LOS path is almost always present, leading to Rician fading, whereas in urban mobile communication, due to the large number of obstacles (e.g. high rise buildings) in the environment, having an LOS path is less probable and Rayleigh fading is expected [22]. On the other hand, UAV-based communication is performed at altitudes that are between these two cases. Thus, the presence and absence of LOS links depends on the position of the UAV. This behavior was modeled in [23] and [24] via the probability of attaining an LOS link between MU k and the UAV which is given as follows where ψ k = 180 π tan −1 ( z d r DM ) is the elevation angle between the UAV and MU k , r DM = (x d − x 0 − r k cos φ k ) 2 + (y d − y 0 − r k sin φ k ) 2 1/2 is the radial distance between the UAV and MU k , and C and B are constants whose values depend on the environment (e.g., rural, urban, and high-rise areas). Eq. (3) suggests that, at higher UAV operating altitudes, the LOS link is more likely to be present, which in turn comes at the expense of a higher path loss. Accordingly, the probability of having a non-LOS (NLOS) connection from the UAV to MU k is given by P NLOS,k = 1 − P LOS,k . Now, based on the LOS probability, the RF channel coefficient, h k,n , may be either LOS or NLOS and is given by [23], [25] where h p k is the free-space path loss, and h r k,n e jΨ k,n and h sr k,n are the Rayleigh and shadowed Rician fading coefficients, respectively. In particular, the path loss is given by where c = f 23.85 and f is the center operating frequency in MHz. The NLOS scenario is characterized by the Rayleigh fading coefficient, h r k,n e jΨ k,n , with power E{ h r k,n 2 } = η 2 and uniformly distributed phase, Ψ ∼ U(0, 2π). The probability density function (pdf) of h r k,n , denoted by f h r k,n (x), is given by Moreover, the distribution of ς = N n=1 |h r k,n | 2 , which characterizes the combined signal after MRC at the UAV, is given by a chi-distribution with 2N degrees of freedom, i.e., ς ∼ χ 2 (2N ). In the LOS case, the channel is affected by both shadowing and small scale fading which is modeled via the Loo model [26], [27], [28]. In this model, the shadowed Rician fading coefficient is is the lognormal shadowed LOS component which is added to the Rayleigh scattering component. Here, the log-normal shadowed component is identical across all RF antennas of the UAV, but the small scale Rayleigh fading is independent across antennas. Unfortunately, the combination of log-normal shadowing and Ricean fading does not lend itself to a closed-form expression for the resulting distribution.
Thus, as a widely accepted approximation [26], [29], the log-normal distribution is fitted to the Nakagami distribution and the shape and spread parameters of the Nakagami pdf are obtained via moment matching to the log-normal pdf as follows where ε = 20 ln (10) and thus, h s k ∼ Nakagami(q, ω). Based on this approximation, the shadowed-Rician pdf is modeled in [26]. In the following Lemma, we characterize the distribution of Lemma 1: For the proposed Nakagami approximation of h s k , the distribution of τ = N n=1 |h sr k,n | 2 is given by where 1 F 1 (·, ·; ·) is the confluent hypergeometric function.
Proof: The proof is given in Appendix A.
Thus, the pdf of h k 2 = N n=1 |h k,n | 2 is given by

C. UAV-Based FSO Channel Model
The UAV-based FSO channel differs from conventional FSO channels with fixed transceivers mounted on top of buildings in the following two aspects. First, in contrast to fixed transceivers, the FSO beam of the UAV is not necessarily orthogonal to the PD plane. Second, the instability of the UAV, i.e., the random vibrations of the UAV, introduces a random power loss. Taking these effects into account, the point-to-point UAV-based FSO channel, g, can be modeled as follows [30], [31] where R s , g p , g a , and g g represent the responsivity of the PD, the atmospheric loss, the atmospheric turbulence, and the geometric and misalignment loss (GML), respectively.

1) Atmospheric Loss:
The atmospheric loss, g p , is due to scattering and absorption of the laser beam by atmospheric particles and is given by [30] g p = 10 − κL 10 , where κ is the attenuation factor, whose value depends on the weather conditions (e.g. clear, is the distance between the UAV and the GS. 2) Atmospheric Turbulence: The atmospheric turbulence, g a , is caused by variations of the refractive index in different layers of the atmosphere due to fluctuations in pressure and temperature. As a universal model that considers both small and large scale irradiance fluctuations, the Gamma-Gamma (GG) distribution, g a ∼ GG(α, β), is considered. Here, α and β are the small and large scale turbulence parameters, respectively, and are given by [32] where ι = (k(2a) 2 /4L) 1/2 , σ 2 R = 1.23C 2 n k 7/6 L 11/6 is the Ryotov variance, C 2 n represents the index of refraction structure parameter, k denotes the wave number, and a is the radius of the PD. The Hufnagle-Valley (H-V) model suggests that C 2 n decreases with increasing height of the UAV (up to 3 km) and is given by is the ground level refraction structure parameter [7], [33]. Thus, the Ryotov variance and accordingly the scintillation variance, var{g a } = 1 α + 1 β + 1 αβ , depend on both the distance between the laser aperture and the PD, L, and the UAV operating height, z d . Now, for conventional applications where UAVs are used to build temporary networks, short distances between the UAV and the GS, e.g., L ≤ 600 m, are expected. In this operating range and typical heights of z d > 30 m, the impact of scintillation becomes very weak even in clear weather condition, i.e., var{g a } ≤ 0.15, and hence for our analysis in Section III, we ignore the effect of turbulence, i.e., g a = E{g a } = 1 is assumed. Then, in Section IV, we use simulations to investigate the impact of this simplifying assumption on the end-to-end achievable rate of the system.
3) GML: The GML, g g , comprises the geometric loss due to the beam spread along the propagation path and the misalignment loss due to the random fluctuations of the position and orientation of the UAV. These random fluctuations are caused by different phenomena, including random air fluctuations around the UAV, internal vibrations of the UAV, and tracking errors, see [31] for a detailed discussion. Also, for the UAV-based FSO channel, the positioning of the UAV may lead to non-orthogonality between the laser beam and the PD plane, which in turn introduces an additional geometric loss. Taking the aforementioned effects into account, the GML for UAV-based FSO channels can be modeled as follows [30] where u is the misalignment factor, w is the beam width at distance L, A 0 = erf(ν min )erf(ν max ), dt is the error function. Assuming perfect tracking, i.e., E{u} = 0, the random variations of the position and orientation of the UAV can be characterized by a Hoyt distributed misalignment factor, i.e., u ∼ H(m, Ω), and accordingly, the pdf of the GML is given by [30] where I 0 (·) is the zero-order modified Bessel function of the first kind, , Ω = λ 1 + λ 2 , and λ 1 and λ 2 are the eigenvalues of matrix Here, denotes the variance of the random fluctuations of the UAV along the position and orientation variables, , and c 5 = − cot(φ d ) cos(θ d ) . Next, given (9) and (13), the pdf of the FSO channel, disregarding the atmospheric turbulence, can be modeled as Remark 1: To shed some light on the impact of the various system parameters on the distribution of the misalignment factor u, we consider special cases. Let us assume σ If the UAV flies in the x-z plane, i.e., y d = θ d = 0, then we obtain Under this assumption, we consider the following special cases for u ∼ H(m, Ω): • UAV flies along the z-axis and x d = x 0 : By increasing |z d − z GS |, m reduces and Ω increases.
• UAV flies along the x-axis and z d = z GS : Now, φ d = π/2 and m = 1 and decreasing x d reduces Ω.
• UAV hovers in front of the GS: y d = 0, z d = z GS , θ d = 0, and φ d = π 2 , then we obtain m = 1, and Ω = 2(σ 2 p + x 2 d σ 2 o ). In this case, the FSO beam is orthogonal to the PD plane and the misalignment factor u is Rayleigh distributed.
In summary, the models for both the RF channel and the FSO channel of UAV-based relay networks differ substantially from the corresponding models for conventional relay networks without UAVs. In particular, the positioning of the UAV affects the path loss and LOS characteristics of the RF access links and the atmospheric loss of the FSO backhaul channel. Furthermore, the instability of the UAV impacts the GML of the FSO backhaul channel. In the following, we study the impact of these effects on the end-to-end achievable rate of the system.

III. END-TO-END ERGODIC SUM RATE ANALYSIS
In this section, the ergodic sum rate, E{C sum }, is adopted as a metric for analyzing the end-toend system performance. Assuming DF relaying at the UAV, the end-to-end achievable sum rate is restricted to the rate of the weaker of the two involved links as a consequence of the max-flow min-cut theorem [34]. In particular, for non-BA relaying, the achievable sum rate depends on the instantaneous fading states of both hops [18]. Thus, the non-BA ergodic sum rate, denoted byC NB sum , is given byC where C RF and C FSO are the instantaneous achievable rates of the RF and FSO links, respectively.
In the BA scenario, the relay is equipped with buffers to store the data received from the MUs and to transmit it when the FSO channel is in a favorable state [17]. Here, for unlimited buffer sizes, the BA ergodic sum rate is given bȳ Remark 2: Although, we assume an unlimited buffer size in (17), in practice, the buffer size is finite and hence, (17) constitutes a performance upper bound for practical BA relaying systems.
However, in [17] and [18], it has been shown that the performance of BA relays with sufficiently large buffer sizes closely approaches the upper bound for unlimited buffer size.
Hence, the ergodic sum rate achieved with BA relaying is an upper bound for the ergodic sum rate of the non-BA case which suggests that buffering data is advantageous for achieving a high ergodic sum rate for the end-to-end system. Nevertheless, this gain comes at the expense of a higher end-to-end delay. Therefore, BA relaying is suitable for delay-tolerant applications.
Next, we analyze the ergodic rate for non-BA (16) and BA (17) relay UAVs for the UAV-based mixed RF-FSO channel model presented in Section II .

A. Ergodic Rate Analysis for BA Relay UAV
The achievable rate for a BA relay UAV depends only on the individual ergodic rates of the RF and FSO channels, i.e., E{C RF } and E{C FSO } in (17). Therefore, in the following, we analyze these ergodic rates separately.

1) Achievable Ergodic Sum Rate of the RF Channel:
Given that the MUs employ orthogonal subchannels, the instantaneous rate of the RF channel, C RF , can be written as a summation of all MUs' rates and is given by where W RF sub is the subchannel bandwidth and R RF k = log 2 (1 + P ζ 2 h k 2 ) is the achievable rate of MU k . Now, the RF ergodic rate is determined by averaging over the random fluctuations of the shadowed Rician and Rayleigh fading in the RF channel and the random MU positions. Taking these effects into account, the ergodic sum rate of the RF channel, denoted byC RF , is given in the following theorem.
Theorem 1: The ergodic sum rate of the RF channel is given as follows where Additionally, the ergodic rates for Rayleigh and shadowed-Rician fading are respectively given as follows where Γ(·, ·) and (x) n denote incomplete Gamma function and the (rising) Pochhammer symbol, respectively [35].
Proof: The proof is given in Appendix B. 2) Achievable Ergodic Rate of the FSO Channel: In the second hop, for an average power constraint,p, the achievable rate of an IM/DD FSO system is given by [21] C FSO = 1 2 where W FSO denotes the FSO bandwidth. Unfortunately, the ergodic rate of the FSO system cannot be computed in closed form for the entire range of signal-to-noise ratios (SNRs). In fact, the expected value of (23) w.r.t. the squared Hoyt variable, u 2 , i.e.,

dx, does not have a closed-form solution and
can only be obtained numerically. Nevertheless, the following theorem presents the ergodic rate for low and high SNRs.

Theorem 2:
The ergodic rate of the FSO system for the low and high SNR regimes, i.e.,γ < 1 andγ 1, respectively, is given bȳ Proof: Please refer to Appendix C.
As will be shown in Section IV, (24a) yields an accurate approximation even if the number Finally, the BA ergodic sum rate in (17) is given bȳ whereC RF ,C FSO low , andC FSO high are given by (20), (24a), and (24b), respectively. Eq. (25) illustrates the inherent trade-off between the ergodic rates of the RF and FSO channels which are dependent on the position and the instability of the hovering UAV. The position of the UAV affects the LOS probability, path loss, and geometric loss of the RF and FSO channels and the instability of the UAV influences the misalignment loss of the FSO channel. Considering these effects, in Section IV, the expression in (25) is employed to optimize the positioning of the BA UAV for maximization of the end-to-end system performance.

B. Ergodic Rate Analysis for Non-BA Relay UAV
In this section, we assume that the UAV is not equipped with a buffer and the instantaneous rate in each channel hop determines the system's end-to-end rate. The non-BA ergodic rate in (16) can be written as [36] where x = min C RF , C FSO . Here, the cumulative distribution function (CDF), F x (x), depends on the CDF of both the RF and FSO channels as follows Here, is the summation of the CDFs of the sum of a squared shadowed Rician RV and a squared Rayleigh RV. Thus, F C RF (u) does not lend itself to a closed-form expression. To cope with this issue, the following lemma is proposed.
Lemma 2: Assuming the number of MUs is sufficiently large, lim K→∞ C RF =C RF , then the ergodic rate for the non-BA relay UAV is given by Proof: Please refer to Appendix D.
In Section IV, we provide simulation results to confirm that even for comparatively small numbers of MUs (e.g.,K = E{K} ≥ 50), (28) is an accurate approximation. Based on (28), two terms are needed to analyze the non-BA ergodic rate, namely F C FSO (C RF ) and E g {C FSO |C FSO ≤ C RF }, which are determined in the following. The CDF of the achievable rate of the FSO channel is a function of the CDF of a squared Hoyt distributed RV, u 2 , [37] and is given by where log 2 (γ + 1), and the Rice-Ie function is defined as where and Q(·, ·) denotes the Marcum Q-function. For the special case, where the FSO beam is orthogonal to the PD plane, i.e., m = 1, (29) becomes the CDF of an exponential distribution. This result is in line with [30], where for an orthogonal beam, misalignment factor u was shown to be Rayleigh distributed, which implies that u 2 is exponentially distributed.
Theorem 3: For low and high SNRs, the average FSO rate, conditioned on the FSO channel being the instantaneous bottleneck of the end-to-end achievable rate, is given by where a = 4 kgw 2 + δ, b = (1−m 4 ) 4m 2 Ω , and δ = (1+m 2 ) 2 4m 2 Ω . Proof: The proof is given in Appendix E.
In Section IV, we will show that the infinite sum in (31a) converges to the exact result if only the first 5 terms are used. Moreover, (31a) and (31b) include finite-range integrals that can be calculated numerically. Furthermore, it can be shown that Theorem 2 is a special case of Theorem 3. In particular, if the RF channel always supports a higher achievable rate than the FSO channel, or equivalently if the conditions in the expected values in (31a) and (31b) are always fulfilled by assumingC RF → ∞, then (31) approaches the ergodic rate of the FSO channel in Theorem 2.
In Section IV, we investigate the accuracy of (32) for small numbers of MUs and employ this expression to study the performance of non-BA relaying UAV-based communications networks and to optimize the position of the non-BA relaying UAV.

IV. SIMULATION RESULTS
In this section, first we use simulation results to validate the analytical results in (25) and (32) for BA and non-BA relay UAVs, respectively. Then, we investigate the impact of the positioning of the UAV on the system performance. Finally, we study the inherent trade-offs in the considered network and the impact of different system and channel parameters on the end-to-end ergodic sum rate. For the considered system and channel models, the parameter values provided in Table I  For the RF channel, the two state channel model in (4) with LOS and NLOS states is adopted.
In the NLOS and LOS states, Rayleigh fading with multipath power η 2 and shadowed Rician fading with lognormal shadowing parameters (µ, σ) are assumed, respectively. The LOS probability parameters, (B, C), in (3) are chosen for an urban environment [24]. For the FSO channel, simulations with and without GG atmospheric turbulence were conducted. The atmospheric loss and GML are incorporated according to (9). The UAV is assumed to be able to track the PD with zero mean misalignment factor, u, and the UAV's instability in the hovering state is accounted for by the position and orientation standard deviations (STD) σ o = 0.3 mrad and σ p = 1 cm, respectively. Given the above assumptions and parameters, we obtained the results reported in this section by averaging over 10 6 realizations of the RF and FSO channels as well as of the MU distributions.

A. Validation of Analytical Results
In this subsection, we investigate the accuracy of the following assumptions and approximations made in our analysis: 1) ignorance of the atmospheric turbulence, g a , in the FSO channel, 2) accuracy ofC NB sum in (28) for finite K, 3) fitting of the lognormal shadowing to Nakagami fading in (7), 4) Gaussian-Chebyshev Quadrature (GCQ) numerical approximation in (20), 5) Taylor series expansions in (25a) and (32a). Furthermore, we confirm our analytical results in (25) and (32) for low and high SNR scenarios. In Section II-C2, we argued that the atmospheric turbulence factor, g a , can be ignored when the UAV flies at typical operating altitudes and for small distances from the GS. Fig. 2 investigates the accuracy of this approximation by comparing the ergodic rates of the FSO channel with and without GG fading as functions of the distance between the UAV and the GS. Here, strong and moderate turbulence conditions 1 , which have different ground level refraction structure parameters, C 2 n (0), and different UAV operating altitudes, z d , are considered. Fig. 2 suggests that at a distance of 1 km, the gap between the curves with and without GG fading at altitudes of 10 m and 30 m is respectively about 7% and 5% for strong turbulence. The smaller gap for the higher altitude is due to the H-V model since C 2 n (z d ) decreases if z d increases. Furthermore, Fig. 2 shows that the gap between the ergodic rates with and without GG fading vanishes for small distances. Considering the dependence of parameters α and β on the distance in (11), the impact of atmospheric turbulence decreases for shorter distances. Since the practical operating 1 We note that the severity of the turbulence affects only the variance of the atmospheric turbulence factor ga, while its mean is always one. In fact, the normalized variance of ga, i.e., the scintillation index var{ga} E{ga} 2 , varies for different operating distances L and different C 2 n (z) [38], and depending on the operating scenario, stronger turbulence can lead to larger/smaller variances than moderate turbulence. This means that, as can be observed in Fig. 2, the ergodic rate of the FSO channel in strong turbulence is not always lower than that in moderate turbulence. range of UAVs is expected to be within one kilometer of the GS and its typical operating altitude is expected to be above 30 m, the impact of atmospheric turbulence can be savely ignored.  Fig. 3 confirms that not taking into account g a does not affect the ergodic rate of the FSO channel and yields accurate results for BA relaying. For non-BA relaying, there is a small gap of about 1% between the simulation results with and without GG fading. Furthermore, in the non-BA case, the simulation results are upper bounded by the analytical results obtained from (32). However, for sufficiently large numbers of MUs or equivalently for sufficiently high MU densities, the instantaneous sum rate of the RF channel approaches the ergodic sum rate of the RF channel (see Lemma 2) and hence, the gap between numerical and analytical results vanishes also for non-BA relaying. Fig. 3 suggests that even for a relatively small number of MUs, this gap is small. For example, for λ = 0.008 or equivalentlyK = 63.3 MU, the gap is only about 2% and 3.4% for w 0 = 0.25  (7) is justified. Overall, we conclude that the analytical ergodic sum rate expressions for BA and non-BA relaying UAVs in (25b) and (32b), respectively, are accurate.

B. Impact of Positioning of the UAV
Next, we investigate the impact of the placement of the UAV on the end-to-end ergodic rate.  height for BA relaying, z BA d , depends on the altitude at which the ergodic sum rate curves of both links intersect, i.e.,C RF =C FSO and is marked by in Fig. 6. On the other hand, for non-BA relaying, the gap between the analytical and simulation results is only 2% and the altitudes that maximize the respective curves, z NB d , are only 8 m apart. Because of the benefits of buffering, BA relaying yields a higher ergodic sum rate than non-BA relaying.
In Fig. 7, the UAV operates at a height of 30 m and moves along the x-axis from the cell center towards the GS. By reducing the distance to the GS, the ergodic rate of the FSO channel drastically increases, due to the exponential reduction of the atmospheric loss. On the other hand, the ergodic rate of the RF channel decreases due to the larger path loss and the smaller LOS probability caused by the smaller elevation angle in (3). Consequently, farther from the cell center, the RF channel is the performance bottleneck and the ergodic rates of both BA and non-BA relaying approach the ergodic sum rate of the RF channel. On the other hand, when the UAV is farther from the GS, the FSO channel limits the performance of both types of relaying.
For BA relaying, the intersection of the RF and FSO ergodic sum rate curves corresponds to the optimal position of the UAV at 11 m, and analysis and simulations yield the same value. On the other hand, comparing the analytical and simulated ergodic rates for non-BA relaying reveals a  confirm that the optimal positioning of the relay depends on the parameters of the RF and FSO channels as well as the type of relaying.

C. Impact of System Parameters
In this subsection, we investigate the impact of the weather-dependent FSO attenuation factor, κ, and the density of the MUs, λ, on the end-to-end system performance. Here, for clarity of presentation, we only show simulated ergodic sum rates. preferable in order to compensate for the reduced ergodic rate in the FSO link. Furthermore, for higher densities λ, the bandwidth available for each MU decreases, and accordingly, the Although, a larger number of MUs does not affect the ergodic rate of the FSO channel, the FSO backhaul has to support the increased rate of the RF channel. Therefore, for larger MU densities, for both BA and non-BA relaying the UAV benefits from moving towards the GS (i.e., increasing |x 0 − x d |) to improve the quality of its backhaul channel.

V. CONCLUSIONS
In this paper, the end-to-end performance of a UAV-based communication system, where a relay UAV connects MUs via RF access links to a GS via an FSO backhaul link, was analyzed in terms of the ergodic sum rate. The UAV's characteristics including its relative position w.r.t. the GS and the MUs and its (in)stability in the hovering state were taken into account in the adopted RF and FSO channel models. In particular, the probability of attaining an LOS path in the RF channel and the GML in the FSO channel were accounted for. Furthermore, to address the application dependent sensitivity to delay, both BA and non-BA relay UAVs were investigated.
Exact and approximate expressions for the ergodic sum rate were derived for BA and non-BA relay UAVs, respectively. We validated the accuracy of the obtained analytical result and investigated the trade-offs affecting the optimum position of BA and non-BA relay UAVs. Our results revealed that the impact of atmospheric turbulence on the quality of the FSO channel can be ignored for practical UAV-GS distances of less than 1 km. Furthermore, the derived approximate analytical expression for the ergodic sum rate for non-BA relays was shown to approach the corresponding simulation results for sufficiently large numbers of MUs. Moreover, our simulations revealed that the random variations of the FSO channel caused by the UAV's instability can be mitigated by BA relaying which results in larger achievable ergodic sum rate compared to non-BA relaying at the expense of introducing an additional delay into the system. Our results also show that when the weather conditions get worse and accordingly the atmospheric loss in the backhaul channel increases, the UAV prefers a position closer to the GS to improve the quality of the backhaul channel. Furthermore, for higher MU densities and the resulting larger amounts of data received via the RF access channel, the end-to-end performance can be improved if the UAV moves closer to the GS in order to enhance the backhaul link quality.
Considering our simulation and analytical results, we conclude that the specific properties of both the FSO backhaul and RF access channels have to be simultaneously taken into account for performance evaluation and optimization of UAV-based communication networks employing mixed RF-FSO channels.
APPENDIX A: PROOF OF LEMMA 1 Given that h sr k,n = h s k + h r k,n e jΨ k,n is a Rician distributed RV with Nakagami distributed parameter h s k ∼ Nakagami(q, ω), the conditional distribution of √ τ is given by a non-central chi distribution as follows Since h s k ∼ Nakagami (q, ω), thus, (h s k ) 2 is a Gamma distributed RV and the pdf of κ 2 is given by Then, we can obtain the unconditional distribution of where 0 F 1 (·, ·) is the confluent hypergeometric function. Thus, the pdf of √ τ is given by Then, using the relation f τ (x) = 1 2 √ , the pdf of τ is obtained as in (7) and this concludes the proof.
APPENDIX B: PROOF OF THEOREM 1 First, the ergodic rate corresponding to (19) can be simplified to a summation of the ergodic rates for the LOS and NLOS states as follows where for equality (a), we exploited the linearity of expectation.
Then, using Campbell's law, which states that E homogeneous Poisson process with intensity λ and f is any nonnegative function [40], we obtain The above integrals can be solved only numerically. To obtain a suitable numerical approximation, we first change variables r k and φ k to x = 2r k r 0 − 1 and y = φ k π − 1, respectively, which yields Then, we note that any definite integral can be transformed to a weighted summation using Gaussian-Chebyshev Quadrature (GCQ) as f (x i ) +Ẽ H , where = π H and x i = cos 2i−1 2H π [41]. Thus, letting f (x, y) = x (1 − x 2 ) × (1 − y 2 ) E hr {R RF k }P NLOS,k +E hsr {R RF k }P LOS,k , (38) can be written as where in (a) and (b) the GCQ is applied to the integrals over x and y, respectively.
Next, we determine the ergodic rate for shadowed Rician fading, denoted by E hsr {R RF k } = E hsr log 2 (1 + γ h k 2 ) as follows where f τ (x) is given in Lemma 1. The above integral can be solved using the identity 1 F 1 (a, b; x) = ∞ n=0 (a)n n!(b)n x n [42] and [39,Eq. (78)]. This leads to the ergodic sum rate for shadowed Rician fading in (22). Then, using [42,Eq. (40)] for the ergodic sum rate for Rayleigh fading leads to (21) and this concludes the proof.

APPENDIX C: PROOF OF THEOREM 2
Based on (23), the ergodic rate of the FSO channel is given by By substituting g = R s g p g g and g g from (12), we obtain 29 Then, for low SNRs, we use the Taylor series ln(1 + x) = ∞ =1 (−1) x (for |x| ≤ 1) to obtain The expectation over the Hoyt-squared variable, u 2 , can be solved by [35, 6.611-4]. This leads to (24a).
where equality holds if C RF = E h,Φ {C RF }. For K → ∞, the instantaneous rate in (19) and the ergodic rate of the RF channel in (20) become identical since where (a) exploits the definition of the ergodic rate. Given this relation, equality holds in (45).
Denoting E h,Φ {C RF } byC RF , we obtain, where in (a) we apply min(c, x) = xPr(x ≤c) +c Pr(c < x), wherec is a constant. In (b), we substitute the definition of the CDF F C FSO (x) = 1 − Pr(C FSO ≥ x). This concludes the proof. 30 APPENDIX E: PROOF OF THEOREM 3 For low SNRs, we can use a conditional version of (43) as follows where in (a), the pdf of a squared Hoyt RV, i.e., f u 2 (x) = 1+m 2 2mΩ exp − (1+m 2 ) 2 x