## 3 Problem Formulation and Solution  
In this section, we aim to search for the optimal PMF of the discrete constellations under a Gamma-Gamma channel. We compute the AIR of the proposed scheme with bit-metric decoding, for various symbol to noise ratios (SNRs). In the following, we review the basic principle of probabilistic shaping. We mainly focus on the AIR of the BMD, rather than symbol-metric decoding (SMD) due to the higher computational complexity of SMD. In this scheme, the input signal \(X\) is composed of \(m\) bit levels \(B=B_{1},\ldots,B_{m}\). The BMD rate is defined as follows [24]:  
\[ \begin{split}\text{R}_{\text{BMD}}&=\left[\mathbb{ H}(X)-\sum\limits_{i=1}^{m}\mathbb{H}(B_{i}\mid Y,G)\right]^{+}\\ &=\left[\mathbb{H}(X)-\sum\limits_{i=1}^{m}\mathbb{E}[\log_{2}\{1 +\exp[-(1-2B_{i})L_{i}]\}]\right]^{+},\end{split}  \tag{6}\]  
where \([\bullet]^{+}\) means the largest number of \(\bullet\) and \(0\), namely \(\max(\bullet,0)\). \(\mathbb{H}(B_{i}\mid Y,G)\) calculates the uncertainty between \(B_{i}\) at \(i\)-th position and the channel output \(Y\) under the turbulence condition. These log-likelihood ratios (LLRs) are computed with the auxiliary channel as:  
\[L_{i}=\log\frac{p_{B_{i}|Y}(0\mid y)}{p_{B_{i}|Y}(1\mid y)}=\log\frac{\sum \limits_{x\in X_{i}^{0}}q_{Y|X}(y\mid x,g)P_{X}(x)}{\sum\limits_{x\in X_{i}^{ 1}}q_{Y|X}(y\mid x,g)P_{X}(x)}, \tag{7}\]  
where \(X_{i}^{0}\) and \(X_{i}^{1}\) denote the sets of PS-QAM symbols whose \(i\)-th bit is \(0\) and \(1\). In an auxiliary channel, the condition probability \(q_{Y|X}(Y\mid X,g)\) in Equation (7) can be given by:  
\[q_{Y|X}(Y\mid X,g)=\frac{1}{\sqrt{\pi}\sigma}\exp\left(\frac{(y-gx)^{2}}{ \sigma^{2}}\right), \tag{8}\]  
The AIR depends on the distribution of the input signal \(X\), for various SNRs. Hence, the objective function at the specified SNR can be written as:  
\[ \begin{split}& C=\max\limits_{P_{X}(x)}\{\text{R}_{\text{BMD}}\}\\ &=\max\limits_{P_{X}(x)}\left[\mathbb{H}(X)-\sum\limits_{i=1}^{m} \mathbb{E}[\log_{2}\{1+\exp[-(1-2B_{i})L_{i}]\}]\right]^{+},\end{split}  \tag{9a}\]\[ \begin{split}\text{s.t.}\langle P_{X}=[P_{X}(x_{1}),\dots,P_{X}(x_{i} ),\dots,P_{X}(x_{M})]\rangle=1;\\ 0\leq P_{X}(x_{i})\leq 1.\end{split}  \tag{9b}\]  
where \(\langle\bullet\rangle\) represents the summation formula. The value of PMF is a non-negative number less than one with a sum of one. As mentioned above, the objective function \(\text{R}_{\text{BMD}}\) is a lack of the closed-form expression. Furthermore, we have to calculate the BMD rate by Monte Carlo simulations. Thus, finding an efficient algorithm that finds the optimal PMF is difficult. Under these circumstances, we adopt the heuristic algorithm to maximize the BMD rate to derive the optimal PMF of PS-QAM signals. In the following, we will introduce the heuristic algorithm. The first stage of the algorithm is initializing the population. We produce individuals with the size of \(P_{pop}\). In this algorithm, each individual contains two vectors, namely position and velocity. Considering the symmetric property of the 2-dim 16-ary QAM constellation, we can initialize the population with 2-dim searching space instead of 16-dim. That is, each position represents the 2-dim probability \(\mathbf{p}_{\text{2\,dim}}\), which can correspond to the PMF of 16 QAM. Generalized \(\mathbf{p}_{\text{2\,dim}}=[p_{1}\ p_{2}]\), then \(\mathbf{P}_{\mathbf{X}}\) can be expressed as \(\mathbf{P}_{\mathbf{X}}= \begin{bmatrix}p_{1}&p_{2}&p_{1}&p_{2}\end{bmatrix} ^{\prime} \times \begin{bmatrix}p_{1}&p_{2}&p_{1}&p_{2}\end{bmatrix} \), where \(\times\) denotes the matrix product. We take \(\mathbf{p}_{\text{2\,dim}}=[0.25,0.25]\) as an example, then \(\mathbf{P}_{\mathbf{X}}=[0.625,\dots,0.625]\). The velocity is also a 2-dim vector formed by the change of each position \(\mathbf{p}_{\text{2\,dim}}\) in the previous and subsequent iterations. Besides, we can initialize the position of the population individuals near the Gaussian distribution in order to achieve faster convergence.