## 2 System Model
Assuming lossless reflection power at all the RISs (i.e., \(\alpha_{k,n}=1,\forall k,n\)), the reflection matrix is given by  
\[\mathbf{\Phi}_{k}=\mathrm{diag}\Big{(}\mathrm{e}^{i\eta_{k,1}},\mathrm{e}^{i\eta_ {k,2}},\cdots,\mathrm{e}^{i\theta_{k,N}}\Big{)},\forall k. \tag{3}\]  
In (1), the signals directly transmitted from the satellite to the UEs are neglected because of severe path loss and shadowing. For a practical air-to-ground (A2G) channel, \(\mathbf{f}_{k}=\sqrt{\xi_{k}}\mathbf{\hat{f}}_{k}\) is generated by correlated Rician fading, which consists of both line-of-sight (LoS) and non-line-of-sight (NLoS) paths, where \(\zeta_{k}\) is the path loss of the RIS-UE link and \(\mathbf{\hat{f}}_{k}\) is Rician fading, as follows:  
\[\zeta_{k,\mathrm{dB}} =-10\log_{10}\Big{(}\frac{4\pi d_{k}^{\mathrm{ru}}}{\lambda} \Big{)}^{2}-\eta_{\mathrm{L}}P_{k}^{\mathrm{L}}-\eta_{\mathrm{NL}}P_{k}^{ \mathrm{NL}}, \tag{4}\] \[\mathbf{\hat{f}}_{k} =\sqrt{\frac{\kappa_{k}}{\kappa_{k}+1}}\mathbf{\hat{f}}_{k}+\sqrt{ \frac{1}{\kappa_{k}+1}}\mathbf{\hat{f}}_{k}. \tag{5}\]  
Here, \(\zeta_{k,\mathrm{dB}}\) is the path loss in dB scale, \(d_{k}^{\mathrm{ru}}\) is the distance from RIS \(k\) to UE \(k\), \(\lambda\) represents the wavelength of a signal, and the average additional losses for LoS and NLoS conditions are denoted as \(\eta_{\mathrm{L}}\) and \(\eta_{\mathrm{NL}}\), respectively. Based on [21], the presence of a LoS path is contingent upon the distance between the RIS and the UE. Then, the probability of LoS is expressed as follows:  
\[P_{k}^{\mathrm{L}}=\frac{1}{1+c_{1}\exp\left[-c_{2}\bigg{(}\arctan\bigg{(} \frac{d_{k}^{\mathrm{ru}}}{d_{k}^{\mathrm{L}}}\bigg{)}\bigg{)}\right]}, \tag{6}\]  
where \(d_{k}^{\mathrm{ru}}\) and \(d_{k}^{\mathrm{h}}\) are, respectively, vertical and horizontal distances from UAV \(k\) and UE \(k\), and \(c_{1}\) and \(c_{2}\) are determined by the particular configuration of the environment [21]. Also, \(P_{k}^{\mathrm{NL}}\) is the probability of NLoS as given by \(P_{k}^{\mathrm{NL}}=1-P_{k}^{\mathrm{L}}\). In (5), \(\kappa_{k}\) is the Rician factor between UE \(k\) and RIS \(k\), and \(\mathbf{\hat{f}}_{k}\in\mathbb{C}^{N\times 1}\), which is composed of deterministic values whose entries are determined depending on the locations of RIS \(k\) and UE \(k\). By defining \(\hat{f}_{k,n}\) as the \(n\)-th element of \(\mathbf{\hat{f}}_{k}\), we have  
\[\mathbf{\hat{f}}_{k,n}=\exp\left(-\mathrm{j}2\pi d_{k,n}^{\mathrm{ru}}/\lambda \right), \tag{7}\]  
where \(d_{k,n}^{\mathrm{ru}}\) is the distance between UE \(k\) and RIS element \(n\) of RIS \(k\)[22]. Considering the spatial correlation in the RIS-UE channel, we define \(\mathbf{\hat{f}}_{k}=\mathbf{G}_{k}^{1/2}\mathbf{\hat{f}}_{k}\), where \(\mathbf{G}_{k}\in\mathbb{C}^{N\times N}\) is the spatial correlation matrix at RIS \(k\) and \(\mathbf{\hat{f}}_{k}\in\mathbb{C}^{N\times 1}\) is composed of NLoS components,with each entry being an i.i.d. complex Gaussian random variable with a mean of zero and a variance of one. For notational simplicity, we define \(\tilde{\beta}_{k}=\sqrt{\zeta_{k}\kappa_{k}/(\kappa_{k}+1)}\) and \(\tilde{\beta}_{k}=\sqrt{\zeta_{k}/(\kappa_{k}+1)}\). In (1), the channel gain vector between the satellite and the RIS, \(\mathcal{g}_{k}\), includes a multibeam pattern, atmospheric fading, and phase rotations, which can be decomposed as follows [23]:  
Figure 1:0  
where \(\mathcal{D}_{k}\in\mathbb{C}^{N\times N}\) is a rain fading matrix and \(\mathcal{b}_{k}\in\mathbb{C}^{N\times 1}\) is the coefficient vector that represents the beam radiation pattern gain with path loss component.