## 2 System and Channel Models
### Channel Model  
In this paper, the channel gain \(h_{m,j},\ m\in\{1,\cdots,M\},j\in\{1,2\}\) is mainly influenced by four factors [19], i.e.,  
\[h_{m,j}=h_{m,j}^{1}h_{m,j}^{\rm a}h_{m,j}^{\rm P}h_{m,j}^{\rm a\infty} \tag{7}\]  
where \(h_{m,j}^{1}\) is atmospheric loss, \(h_{m,j}^{\rm a}\) is atmospheric turbulence, \(h_{m,j}^{\rm P}\) is pointing error, and \(h_{m,j}^{\rm a\infty}\) is AoA fluctuation. According to the exponential Beers-Lambert law, the atmospheric loss can be modeled as [19]  
\[h_{m,j}^{1}=\exp\bigl{(}-Z_{m,j}\xi\bigr{)} \tag{8}\]  
Figure 1: Model of FSO communication system based on parallel UAV relays. where \(Z_{m,j}\) is the length of the \(j\)th hop of the \(m\)th path and \(\xi\) represents the atmospheric attenuation coefficient. For atmospheric turbulence, the Gamma-Gamma distribution is employed since it is suitable for all turbulence conditions [25]. The probability density function (PDF) of \(h_{m,j}^{\mathrm{a}}\) can be expressed as [26]  
\[f_{h_{m,j}^{\mathrm{a}}}\left(h_{m,j}^{\mathrm{a}}\right)=\frac{2(\alpha_{m,j} \beta_{m,j})^{\frac{\alpha_{m,j}+\beta_{m,j}}{2}}}{\Gamma\left(\alpha_{m,j} \right)\Gamma\left(\beta_{m,j}\right)}\left(h_{m,j}^{\mathrm{a}}\right)^{ \frac{\alpha_{m,j}+\beta_{m,j}}{2}-1}K_{\alpha_{m,j}-\beta_{m,j}}\left(2\sqrt{ \alpha_{m,j}\beta_{m,j}h_{m,j}^{\mathrm{a}}}\right),h_{m,j}^{\mathrm{a}}>0 \tag{9}\]  
where \(\alpha_{m,j}\) and \(\beta_{m,j}\) are the effective numbers of large-scale and small-scale eddies in the scattering process and given by  
\[\left\{ \begin{array}{l}\alpha_{m,j}=\left[\exp\left(0.49\sigma_{\mathrm{R},m,j}^{2}\Big{/}\left(1+1.11\sigma_{\mathrm{R},m,j}^{12/5}\right)^{7/6}\right)- 1\right]^{-1}\\ \beta_{m,j}=\left[\exp\left(0.51\sigma_{\mathrm{R},m,j}^{2}\Big{/}\left(1+0.69 \sigma_{\mathrm{R},m,j}^{12/5}\right)^{5/6}\right)-1\right]^{-1}\end{array} \right. \tag{10}\]  
where \(\sigma_{\mathrm{R},m,j}^{2}=1.23C_{n}^{2}k^{7/6}Z_{m,j}^{11/6}\) denotes the Rytov variance, \(C_{n}^{2}\) denotes the index of refraction structure parameter, \(k=2\pi/\lambda\) denotes the number of light waves, and \(\lambda\) denotes the wavelength. Based on [10], we assume that the random variances in position at each UAV are Gaussian-distributed with zero mean and variance \(\sigma_{\mathrm{p},\mathrm{u}}^{2}\), the random variances in the direction deviation of the UAVs are Gaussian-distributed with zero mean and variance \(\sigma_{\mathrm{o}}^{2}\), and the random variances in position deviation at the ground platform are Gaussian-distributed with zero mean and variance \(\sigma_{\mathrm{p},\mathrm{g}}^{2}\). Therefore, the PDF of the pointing error is written as  
\[f_{h_{m,j}^{\mathrm{p}}}\left(h_{m,j}^{\mathrm{p}}\right)=\frac{\xi_{m,j}^{2}} {A_{m,j}^{\mathrm{e}}}\left(h_{m,j}^{\mathrm{p}}\right)^{\xi_{m,j}^{2}-1},0 \leq h_{m,j}^{\mathrm{p}}\leq A_{m,j} \tag{11}\]  
where \(\xi_{m,j}^{2}=\omega_{Z_{\mathrm{eq},m,j}}^{2}\Big{/}\Big{(}4\sigma_{\mathrm{ s},m,j}^{2}\Big{)}\); \(A_{m,j}=\left[\mathrm{erf}(v_{m,j})\right]^{2}\) is the maximal fraction of collected power, with \(v_{m,j}=\sqrt{\pi}a_{m,j}/\Big{(}\sqrt{2}\omega_{Z_{m,j}}\Big{)}\); \(a_{m,j}\) is the aperture radius; \(\omega_{Z_{m,j}}\) is the beam waist radius at distance \(Z_{m,j}\); and \(\omega_{Z_{\mathrm{eq},m,j}}^{2}=\omega_{Z_{m,j}}^{2}\sqrt{\pi}\mathrm{erf} \big{(}v_{m,j}\big{)}\Big{/}\Big{(}2v_{m,j}\exp\left(-v_{m,j}^{2}\right)\Big{)}\) is the equivalent beam waist. Moreover, \(\sigma_{\mathrm{s},m,j}^{2}\) is the variance in the total displacement for different links, which is given by  
\[\sigma_{\mathrm{s},m,j}^{2}=\left\{ \begin{array}{l}\sigma_{\mathrm{p}, \mathrm{g}}^{2}+\sigma_{\mathrm{p},\mathrm{u}}^{2},j=1\\ \sigma_{\mathrm{p},\mathrm{g}}^{2}+\sigma_{\mathrm{p},\mathrm{u}}^{2}+Z_{m,j} ^{2}\sigma_{\mathrm{o}}^{2},j=2\end{array} \right. \tag{12}\]  
From [9; 27], it is known that the link between transmitter and receiver is interrupted when the AoA of the received beam falls outside the FoV of the receiver.