*1.3. Contributions*

Despite their relevance, none of the aforementioned studies, e.g., [6,29–33], consider two important aspects: (i) the random nature of the EH, which is inherited from the variability of the optical channel; and (ii) the trade-off between bit error rate (BER) and EH. In this paper, we investigate the influence of this random behavior on the design of EHoptimized links, while fulfilling a target BER. This will be achieved by transmitting a certain fraction, *ζ*, of the peak transmission power, *Pm*, with the logical symbol '0'. As shown in this paper, a proper design of the *ζ* maximizes the harvested energy while fulfilling a desirable reliability in terms of the BER.

More specifically, we aim to respond to the following inquiries: (i) the optimum *ζ* value to maximize the harvested energy while fulfilling a target BER; and (ii) the maximum average harvested energy that can be achieved under realistic link conditions. To address those two issues we derive some closed-form expressions for the two main metrics here considered: the average BER and the average harvested energy. Therefore, the paper's contributions can be summarized as follows:


The rest of the paper is structured as follows. In Section 2 the system model is mathematically described, including the signal as the channel model. Then, the main metrics are presented and derived in Section 3, namely, the expressions for the amount of harvested energy by the system and its BER. In Section 4 numerical results are obtained to assess the benefits of the proposed scheme, which allows to identify both trends and limiting factors. Finally, relevant conclusions are discussed in Section 5.

## **2. System Model**

## *2.1. Proposed Scenario*

We consider a flying backhaul/fronthaul network composed of UAVs that provides reliable connectivity to a 5G+ radio access network (RAN). The UAVs act as networked flying platform (NFP) nodes that can react to changes in weather or traffic conditions [34]. Its ground-to-air links are based on FSO technology with EH. This solution has four remarkable benefits: (i) its FSO links provide the required high data rates for the backhaul or fronthaul connections; (ii) FSO links do not cause interference to the 5G RAN;

(iii) UAVs can adapt their location to the traffic and channel conditions; and (iv) EH extends the UAVs' service time and thus, improves the network performance.

This approach is illustrated with Figure 1 where it shows the representative nodes and connections. The backhaul traffic represents the IP data transmissions between the 5G core network and the BSs in distributed RAN approaches [35]. Accordingly, the BSs fully implement the 5G RAN protocol stack, and thus they are source/destination of IP packets towards/from the core network.

**Figure 1.** Possible 5G + network scenarios where FSO will extend or complement existing deployments.

A recent approach to the classical distributed RAN arises with the cloud RAN (C-RAN) [36], where a base band unit (BBU) centralizes the RAN processing of many small cells. The signal transmission related to each cell is carried out by radio remote heads (RRHs) that are connected with the BBUs through fronthaul links. This latter approach offers two main benefits compared to distributed RAN. First, C-RAN reduces costs since each small cell can be implemented with a RRH, cheaper than a complete BSs. Second, C-RANs ease the implementation of coordinated mechanisms for interference mitigation.

Coming back to Figure 1, it is straightforward to see the different backhaul links between the core network and the BSs through the FSO-UAVs based network. This scheme requires a FSO gateway that creates the link with a mother NFP node, i.e., a mother UAV that forwards the IP traffic towards the NFP node currently connected to the target BS. Analogously, the fronthaul links between the BBU and the RRHs also require a connection between a FSO gateway and the mother NFP node. Certainly, many applications beyond classical mobile networks can benefit from the proposed scenario shown in Figure 1 since it can provide backhaul/fronthaul connectivity to 5G V2X, or even future underwater communications [11–13,37].

#### *2.2. Received Signal Model*

Throughout this paper, we consider a non-coherent intensity modulation with directdetection (IM/DD) scheme since it is a practical and widely extended technique for FSO communications [38]. With this scheme, the intensity of the emitted light is employed to convey the information. Therefore, the PD output current can be expressed as follows [39,40]

$$
\dot{n} = h\mathbf{R}\mathbf{x} + n\_\prime \tag{1}
$$

where *x* is the transmitted intensity in W, *R* is the detector responsivity in W/A, *h* is the channel coefficient and *n* represents the noise term at the receiver side. We assume that the

shot noise induced by the ambient light is the dominant noise source in our analysis. Thus, we model *n* as a signal-independent zero-mean Gaussian noise with variance given by

$$
\sigma\_n^2 = 2qRB\_nP\_{amb\nu} \tag{2}
$$

where *q* denotes the electron charge, *Bn* is the receiver noise bandwidth, and *Pamb* is the ambient light power. This latter term can be obtained as *Pamb* = *<sup>π</sup>r*<sup>2</sup>*SnBo*Ω*FOV*, where *r* is the receiver aperture radius, *Sn* is the ambient light spectral radiance, *Bopt* is the optical bandwidth in nm and Ω*FOV* is the receiver field of view (FOV) in srad.

The transmitted intensity is taken as symbols drawn with equal probability from the OOK constellation such that *x* ∈ { *P*1, *<sup>P</sup>*0}, where *P*1 and *P*0 represent the power corresponding to the transmission of a '1' and a '0', respectively.

In this respect, we can cite three primary atmospheric phenomena that affect optical wave propagation and that constitute the total channel coefficient *h*: (i) the deterministic path loss, *hl* characterized by absorption and scattering; (ii) the geometric spread and pointing errors *hp*; and (iii) the refractive-index fluctuations (i.e., the atmospheric turbulence), *ha*, leading to irradiance fluctuations. Thus, the total channel attenuation is modeled as the product of these aforementioned channel factors as:

$$h = h\_l h\_p h\_{a\nu} \tag{3}$$

where *hp* and *ha* are random variables (RVs). In the next sections, we describe the underlying distribution that model those RVs. Finally, we define the signal to noise ratio as [41]:

$$\text{SNR} = \frac{\left(RhP\_{av}\right)^2}{\sigma\_{\text{ul}}^2},\tag{4}$$

where *Pav* = ( *P*0 + *<sup>P</sup>*1)/2 represents the average transmit power.

#### *2.3. Atmospheric Turbulence Model*

In order to model the intensity fluctuations caused by the atmospheric turbulence, the statistical Gamma-Gamma distribution is here assumed because of its mathematical tractability and accuracy to characterize a wide variety of scenarios ranging from weak to strong turbulence. Thus, and following [21], the probability density function (pdf) of *ha* is written as

$$f\_{h\_a}(h\_a) = \frac{2(a\beta)^{(a+\beta)/2}}{\Gamma(a)\Gamma(\beta)} h\_a^{(a+\beta)/2-1} K\_{a-\beta} \left( 2(a\beta h\_a)^{1/2} \right),\tag{5}$$

where *Kp*(*x*) is the modified Bessel function of the second kind, and <sup>Γ</sup>(*x*) is the Gamma function, with *α* representing the effective number of large-scale cells of the scattering process and with *β* denoting the effective number of small-scale cells. Namely, from [39], these latter parameters, *α* and *β*, can be obtained as

$$\mathfrak{a} = \left[ \exp \left( 0.49 \sigma\_R^2 (1 + 1.11 \sigma\_R^{12/5})^{-7/6} \right) - 1 \right]^{-1},\tag{6}$$

and

$$\beta = \left[ \exp\left( 0.51 \sigma\_R^2 (1 + 0.69 r\_R^{12/5})^{-5/6} \right) - 1 \right]^{-1},\tag{7}$$

respectively, where *σ*<sup>2</sup> *R* is the Rytov variance which, for uplinks paths, is defined as [39]

$$
\sigma\_R^2 = 2.25k^{7/6} (Z - z\_0)^{5/6} \int\_{z\_0}^Z \mathbb{C}\_n^2(z) \left(\frac{z - z\_0}{Z - z\_0}\right)^{5/6} dz,\tag{8}
$$

with *C*<sup>2</sup> *n*(*z*) being the index of refraction structure parameter at altitude *z*, whereas *k* = 2*π*/*λ* is the optical wave number and *Z* and *z*0 are the UAV and transmitter heights, respectively.

## *2.4. Atmospheric Attenuation*

For optical waves, the effect of the atmospheric attenuation suffered by the light propagating through the atmosphere is mainly caused by either absorption as scattering by air molecules in addition to both absorption and scattering by solid or liquid particles suspended in the air which, as a last resort, indicates the effect of weather conditions on the transmitted laser beam. Absorption and scattering are often grouped together under the topic of extinction, defined as the reduction or attenuation in the amount of radiation passing through the atmosphere. Mathematically speaking, such attenuation is incorporated using the well-known Beer–Lambert law [42–44] given by:

$$h\_l = \exp\left(-a(\lambda)L\right),\tag{9}$$

where *a*(*λ*) is the wavelength dependent attenuation coefficient (extinction coefficient) and *L* = *Z* − *z*0 is the propagation path length from the transmitter. As commented, the coefficient *a* depends on weather conditions and can be obtained from the atmospheric visibility. Since the attenuation coefficient hardly changes over long periods of time, we have assumed *hl*as a deterministic coefficient in our analysis.

#### *2.5. Geometric Spread and Pointing Error Model*

In addition to attenuation and atmospheric turbulence, geometric beam spread and pointing accuracy also affect the performance of these systems. The geometric loss is caused by the divergence of the transmit beam when propagating through the atmosphere, as *ωz* = *θT* · *L*, with *ωz* being the received beam waist, with *θT* denoting the transmitter divergence angle, whereas *L* is the propagation path length. Since the received beam width is usually wider than the lens aperture size, part of the transmitted power cannot be collected, leading to loss. Thus, this geometric loss depends mainly on the ratio between the received beam waist, *ωz*, and the receiver aperture radius, *r*. However, during the designing of a FSO link, it is possible to control the beamwidth produced at a certain distance, and therefore the ratio *<sup>ω</sup>z*/*<sup>r</sup>*, by adjusting properly the laser parameters.

On the other hand, imperfections in the pointing, acquisition, and tracking process between the ground stations and the UAVs can also cause loss. Note that both phenomena are interrelated and, thus, accurate modeling frameworks have been proposed to account for both effects. To this end, in this work, we consider the general model proposed in [44,45]. According to this model, a Gaussian beam profile is assumed with a beam waist, *ωz*, on the receiver plane and a circular aperture receiver of radius *r*; then, the attenuation due to the geometric spread with pointing error can be approximated as the Gaussian form

$$
\hbar\_p \approx A\_o \exp\left(-\frac{2\rho^2}{\omega\_{zcq}^2}\right),
\tag{10}
$$

where *ρ*, is the radial pointing error, *Ao* is the fraction of collected power without pointing error, i.e., only due to geometric spread, and *<sup>ω</sup>*2*zeq* is the equivalent beam width. Here, *Ao* and *<sup>ω</sup>*2*zeq* are given by *Ao* = erf(*μ*)<sup>2</sup> and *<sup>ω</sup>*2*zeq* = *<sup>ω</sup>*2*z*√*π*erf(*μ*)/[<sup>2</sup>*<sup>μ</sup>* exp(−*μ*<sup>2</sup>)], respectively, being erf(·) the error function and *μ* = √*<sup>π</sup>r*/(√<sup>2</sup>*ωz*). Moreover, considering independent identical Gaussian distribution for the horizontal *x* and *y* displacement in the receiver plane, the radial error *ρ* = *x*<sup>2</sup> + *y*2 is modeled as Rayleigh distribution with a jitter variance at the receiver *σ*2*s* . Under these assumptions, the channel attenuation, *hp*, can be seen as a function of the radial displacement, *ρ*, which is a RV. Hence, the pdf of *hp* can be seen as a random variable transformation problem, which leads to the following expression [45]:

$$f\_{\hbar\_p}(\hbar\_p) = \frac{\gamma^2}{A\_o^{\gamma^2}} \hbar\_p^{\gamma^2 - 1},\tag{11}$$

where *γ* = *<sup>ω</sup>zeq*/(<sup>2</sup>*σs*) denotes the ratio between the equivalent beam radius at the receiver and the pointing error displacement standard deviation.

#### *2.6. Composite Channel Model*

Once the three factors included in Equation (3) have been individually discussed, the statistical characterization of the composite channel *h* = *hahlhp* can be achieved. Thus, from [45],

$$f\_h(h) = \frac{2\gamma^2 (a\beta)^{\frac{a+\beta}{2}}}{(A\_o h\_l)^{\gamma^2} \Gamma(a) \Gamma(\beta)} h^{\gamma^2 - 1} \int\_{\frac{h}{A\_o h\_l}}^{\infty} h\_a^{\frac{a+\beta}{2} - 1 - \gamma^2} K\_{a-\beta}(2\sqrt{a\beta h\_a}) \, d\mathbf{h}\_a. \tag{12}$$

where *α* and *β* parameters include the information on the strength of the turbulence, with *γ* containing the severity of the pointing error, where *Ao* denotes the geometric spread attenuation and with *hl* being the path loss.
