## 2 Statistical Characteristics of \(\alpha\)-\(\mu\)/IGA Distribution
### Bivariate Case
When multipath fading is superimposed on shadowing, \(W_{i}\) slowly varies. Here, we assume \(W_{i}\) follows the IGA distribution in (4) and \(\mathbb{E}\big{[}r_{i}^{2}\big{]}=\Omega_{i}\). By using a simple operation of Equation (8) given in [44], the joint PDF of the bivariate IGA distribution can be derived as  
\[f_{W_{1},W_{2}}(w_{1},w_{2})=\frac{(n-1)^{n+1}(w_{1}w_{2})^{-0.5(n+3)}}{\Gamma (n)(1-\rho_{G})\rho_{G}^{(n-1)/2}}\exp\left[-\frac{(n-1)(w_{1}^{-1}+w_{2}^{-1}) }{1-\rho_{G}}\right]I_{n-1}\left[\frac{2(n-1)\sqrt{\rho_{G}}}{\sqrt{w_{1}w_{2 }}(1-\rho_{G})}\right], \tag{15}\]  
where \(n>1\) is the IGA shaping parameter, \(\rho_{G}\) is the power correlation coefficient between \(w_{1}^{2}\) and \(w_{2}^{2}\), and \(\mathbb{E}[w_{i}]=1\). Thus, by using the same approach shown in [29] (Equation (5)), and performing some basic operations, the joint PDF of the bivariate \(\alpha\)-\(\mu\)/IGA composite distribution is written as  
\[f_{R_{1},R_{2}}(r_{1},r_{2})=\sum\limits_{k=0}^{\infty}\sum\limits_{l=0}^{ \infty}\frac{\Phi p_{N}^{k}p_{G}^{l}r_{1}^{-1}r_{2}^{-1}}{kl!\Gamma(k+\mu) \Gamma(l+n)}{\prod\limits_{i=1}^{2}H_{1,1}^{1,1}[\Xi_{i}r_{i}^{\alpha_{i}}|_{ (\mu+k,1)}^{(1-n-l,0.5\alpha_{i})}]}, \tag{16}\]  
where \(\Phi=\alpha_{1}\alpha_{2}(1-\rho_{G})^{n}(1-\rho_{N})^{\mu}/[\Gamma(\mu)\Gamma (n)]\), \(\Xi_{i}=(1-\rho_{G})\beta_{i})^{0.5\alpha_{i}}/[(1-\rho_{N})(\Omega_{i}(n-1))^ {0.5\alpha_{i}}]\). By setting \(\gamma_{i}=r_{i}^{2}E_{s}/N_{0}\) and \(\overline{\gamma}_{i}=\Omega_{i}E_{s}/N_{0}\), after a simple variable transformation, the joint PDF of \(\gamma_{1}\) and \(\gamma_{2}\) over the correlated \(\alpha\)-\(\mu\)/IGA channels can be obtained as  
\[f_{\gamma_{1},\gamma_{2}}(\gamma_{1},\gamma_{2})=\sum\limits_{k=0}^{\infty} \sum\limits_{l=0}^{\infty}\frac{0.25\Phi p_{N}^{k}p_{G}^{l}\gamma_{1}^{-1} \gamma_{2}^{-1}}{kl!\Gamma(k+\mu)\Gamma(l+n)}{\prod\limits_{i=1}^{2}H_{1,1}^{1,1}[\phi_{i}\gamma_{i}^{0.5\alpha_{i}}|_{(\mu+k,1)}^{(1-n-l,0.5\alpha_{i})}]}, \tag{17}\]  
where \(\phi_{i}=((1-\rho_{G})\beta_{i})^{0.5\alpha_{i}}/[(1-\rho_{N})((n-1)\overline{ \gamma}_{i})^{0.5\alpha_{i}}]\). Hence, with the aid of [37] (Equation (3.194.1)), the corresponding joint CDF of \(\gamma_{1}\) and \(\gamma_{2}\) can be yielded as  
\[ \begin{array}{ll}F_{\gamma_{1},\gamma_{2}}(\gamma_{1},\gamma_{2})&=\int_{0 }^{\gamma_{1}}\int_{0}^{\gamma_{2}}f_{\gamma_{1},\gamma_{2}}(\gamma_{1},\gamma _{2})d\gamma_{1}d\gamma_{2}\\ &=\sum\limits_{k=0}^{\infty}\sum\limits_{l=0}^{\infty}\frac{0.25\Phi p_{N}^{k }p_{G}^{l}}{kl!\Gamma(k+\mu)\Gamma(l+n)}{\prod\limits_{i=1}^{2}H_{2,2}^{1,2}[ \theta_{i}\gamma_{i}^{0.5\alpha_{i}}|_{(\mu+k,1),(0.5\alpha_{i})}^{(1,0.5\alpha _{i}),(1-n-l,0.5\alpha_{i})}]}.\end{array}  \tag{18}\]  
Moreover, based on the definition of the joint MGF in [1], the joint MGF of \(\gamma_{1}\) and \(\gamma_{2}\) can be given by  
\[ \begin{array}{ll}M_{\gamma_{1},\gamma_{2}}(s_{1},s_{2})&=\mathbb{E}[\exp(-s _{1}\gamma_{1}-s_{2}\gamma_{2})]\\ &=\sum\limits_{k=0}^{\infty}\sum\limits_{l=0}^{\infty}\frac{0.25\Phi p_{N}^{k }p_{G}^{l}}{kl!(k+\mu)\Gamma(l+n)}{\prod\limits_{i=1}^{2}H_{2,1}^{1,2}[ \theta_{i}s_{i}^{-0.5\alpha_{i}}|_{(\mu+k,1)}^{(1,0.5\alpha_{i}),(1-n-l,0.5 \alpha_{i})}]}.\end{array}  \tag{19}\]