## 2 System and Channel Model
### HAP-to-Terrestrial Link  
In the HAP-to-terrestrial link, the transmitter first transmits the optical signal to the RIS module with \(N_{r}\) elements. Subsequently, each RIS element reflects the optical signal to the receiver. During this process, the optical signal is affected not only by atmospheric turbulence but also by pointing errors due to beam jitter and RIS jitter. Assuming that the receiver can receive the energy of all beams, the instantaneous SNR of the link can be expressed as [28]  
\[\gamma=\sum_{k=1}^{N_{r}}\gamma_{k}=\overline{\gamma}\sum_{k=1}^{N_{r}}l_{k}^{ 2} \tag{6}\]  
where \(\gamma_{k}\) is the instantaneous SNR of the k-th channel, \(\overline{\gamma}\) is the average SNR, \(l_{k}=I_{a}I_{p}\), \(I_{a}\) is the atmospheric turbulence, and \(I_{p}\) denotes the pointing errors. Figure 2 shows the pointing errors \(I_{p}\) caused by beam jitter and RIS jitter [29]. \(\beta_{k}\) is the reflection error angle of RIS, \(\theta_{k}^{{}^{\prime}}=\sqrt{\theta_{x_{k}}^{2}+\theta_{2_{k}}^{2}}\) is the superimposed pointing error angle, \(\theta_{x_{k}}^{\prime 2}=\Big{(}1+\frac{L_{1}}{L_{2}}\Big{)}\theta_{x_{k}}+2 \beta_{x_{k}}\), \(\theta_{y_{k}}^{\prime 2}=\Big{(}1+\frac{L_{1}}{L_{2}}\Big{)}\theta_{y_{k}}+2 \beta_{y_{k}}\), \(\{\theta_{x_{k}},\theta_{y_{k}}\}\sim N\big{(}0,\sigma_{g}^{2}\big{)}\), \(\{\beta_{x_{k}},\beta_{y_{k}}\}\sim N\Big{(}0,\sigma_{g}^{2}\big{)}\), \(L_{1}\) is the distance from HAP to RIS, \(L_{2}\) is the distance from the RIS to the receiver, \(L=L_{1}+L_{2}\), \(\theta_{x_{k}}\) and \(\beta_{y_{k}}\) are the error angle in the horizontal direction, \(\theta_{y_{k}}\) and \(\beta_{y_{k}}\) are the error angle in the vertical direction, and \(N(0,\sigma)\) is a normal distribution with mean 0 and variance \(\sigma\). Figure 2: Schematic diagram of beam jitter and RIS jitter. The PDF of \(h_{p}\) can be expressed as [29]  
\[f_{I_{p}}\big{(}I_{p}\big{)}=\frac{\xi}{A_{0}}\bigg{(}\frac{I_{p}}{A_{0}}\bigg{)}^ {\xi-1} \tag{7}\]  
where \(\xi=\frac{\omega_{zeq}^{2}}{4\sigma_{e}^{2}I_{2}+16\varepsilon\theta_{z}^{2}I_ {2}^{2}}\), \(A_{0}=erfc^{2}(v)\), \(\omega_{zeq}^{2}=\omega_{z}^{2}\frac{\sqrt{\pi}erf(v)}{2\varepsilon v^{-2}}\), \(v=\sqrt{\frac{\pi}{2}}\frac{a}{\omega_{z}}\), \(a\) is the size of the receiving aperture, \(\omega_{z}\) is the beam width, and \(erfc(\bullet)\) is the supplementary error function. Subsequently, the PDF of \(I_{k}\) can be expressed as  
\[f_{I_{k}}(I)=\int_{I|A_{0}}^{\infty}\frac{1}{I_{a}}f_{I_{p}}\bigg{(}\frac{I}{I _{a}}\bigg{)}f_{I_{k}}(I_{a})dI_{a} \tag{8}\]  
Substituting Equations (2) and (7) into Equation (8), and then performing a few transformations gives the PDF of \(I_{k}^{2}\) as follows:  
\[f_{I_{k}^{2}}(I)=\frac{A_{\xi}^{2}}{2}\sum\limits_{k=1}^{\beta}a_{k}B^{-\frac {a+k}{2}}G_{1,3}^{3,0}\Bigg{[}\frac{\sqrt{I}}{BA_{0}}\Big{|}\frac{1+\xi}{|\xi, a,k}\Bigg{]} \tag{9}\]  
The mean and variance of \(I_{k}^{2}\) can be expressed as  
\[\mu=\frac{A_{0\xi}^{2}A}{2}\sum\limits_{k=1}^{\beta}a_{k}B^{2+\frac{a+k}{2}} \frac{\Gamma(2+\xi)\Gamma(2+\alpha)\Gamma(2+k)}{\Gamma(3+\xi)} \tag{10}\]  
\[\sigma^{2}=\frac{A_{0\xi}^{4}A}{2}\sum\limits_{k=1}^{\beta}a_{k}B^{4+\frac{a+ k}{2}}\frac{\Gamma(4+\xi)\Gamma(4+\alpha)\Gamma(4+k)}{\Gamma(5+\xi)}-\mu^{2} \tag{11}\]  
Using the Central Limit Theorem (CLT) [30] and assuming a large number of \(N_{r}\), \(\sum\limits_{k=1}^{N_{r}}I_{k}^{2}\) can be approximated by a Gaussian random variable with mean \(N_{r}\times\mu\) and variance \(N_{r}\times\sigma^{2}\). Therefore, the PDF of \(\gamma_{p-t}\) can be expressed as  
\[f_{p-t}(\gamma)=\frac{1}{\sqrt{2\pi N_{r}\sigma^{2}\gamma}}\exp\Bigg{(}-\frac {(\gamma-\overline{\gamma}N_{r}\mu)^{2}}{2\overline{\gamma}^{2}N_{r}\sigma^{ 2}}\Bigg{)} \tag{12}\]  
Considering \(\int_{u}^{\infty}\exp\left(-\frac{\gamma^{2}}{4\beta}-\gamma x\right)dx= \sqrt{\pi\beta}e^{\beta\gamma^{2}}\bigg{[}1-\Phi\bigg{(}\gamma\sqrt{\beta}+ \frac{u}{2\sqrt{\beta}}\bigg{)}\bigg{]}\), \(\Phi(\bullet)\) is the error function.