A Unified Statistical Model for Atmospheric Turbulence-Induced Fading in Orbital Angular Momentum Multiplexed FSO Systems

Space division multiplexing (SDM) technique allows to increase the capacity of a communication system by transmitting several independent data streams in parallel. In optical communications, SDM can be realized by sending different spatial modes in a multi-mode fiber or using many cores in a multi-core fiber [1]. In analogy to optical fiber transmission systems, orbital angular momentum (OAM) is proposed to transmit multiple signals over free space channels [2, 3]. This simultaneous transmission of information on OAM modes is possible thanks to the orthogonality property of OAM beams that allows propagation without interference between signals. A simple scheme of OAM FSO transmission consists in using only the intensity of OAM beams to carry the information. At the receiver, photo-detectors are used to detect light intensity, and the system is operating in intensity-modulation direct-detection (IM/DD). This configuration has a low cost and is more feasible for practical deployment. On the other hand, coherent detection mostly used with optical fibers can also be used for OAM FSO systems. In coherent systems, the spectral efficiency can be further increased by using the two polarizations of light and higher order modulation formats. However, these benefits come at the cost of the high receiver expense and complexity that increases with the number of modes. By using OAM multiplexing as an additional degree of freedom to polarization and wavelength, coherent laboratory experiment has achieved more than 1 Pbit/s capacity using 26 modes [4].
In real life communication scenarios, transmitted OAM beams are subject to atmospheric turbulence induced distortions that if not properly addressed and compensated may deteriorate the entire FSO link [5, 6, 7]. Several experimental and simulation works have studied the effect of atmospheric turbulence and different techniques to emulate turbulence have been proposed. Atmospheric turbulence can be emulated by using a rotatable phase plate that has a pseudo-random distribution that obeys Kolmogorov statistics. This technique allows to control the turbulence strength and was used in [8], at the transmitter for a 2-OAM transmission in a 100-m round trip FSO link. The same process was used in [9] to emulate turbulence in a 4-OAM modes transmission. An alternative to rotatable phase plate is spatial light modulator (SLM). In fact, the distorted phase hologram of a specific OAM mode after propagation can be generated by computer and printed on a SLM. The latter is used at the receiver as a demodulating SLM that contains the effect of atmospheric turbulence [10]. Besides, propagation of OAM modes through the turbulent atmosphere has also been subject to numerical simulation studies. This can be realized using the split-step Fourier propagation method where turbulence is induced along the propagation path [7, 11, 12, 13, 14]. This technique is very accurate and can simulate the propagation of any OAM mode for different turbulence regimes. Nonetheless, it is computationally complex as it involves many multiplication operations of matrices having large dimensions. Therefore it is important to have a statistical model to describe the propagation of OAM beams in turbulent free space.

In classical FSO systems employing only the Gaussian beam for transmission, several statistical models were proposed for irradiance fluctuations. For example, a Log-Normal distribution can be used to model weak turbulence variations [15]. For moderate, to strong turbulence, a Double-Weibull stochastic model was proposed in [16]. Moreover, the Gamma-Gamma distribution is a widely accepted model used for all levels of turbulence [17]. For all the previous models, closed-form expressions for performance metrics such as the average capacity, the outage probability, and the bit-error rate (BER) were derived. In OAM FSO systems, to the best of our knowledge, the only proposed statistical channel models for turbulence-induced fading were given in [18, 19]. In this models, the self-channel irradiance of OAM modes was shown to obey a Johnson $S_{B}$ distribution whereas the crosstalk between OAM modes follows an Exponential distribution. Knowing that the statistical properties of the Johnson $S_{B}$ distribution are analytically intractable, it is not possible to extract closed-form expressions for performance metrics. Moreover, the lack of a unified turbulence model makes it difficult to assess the effect of the interference on the performance of OAM modes. In [20], the interference fading between OAM modes was assumed to be part of the added Gaussian noise. This consideration has allowed deriving a theoretical BER expression but was unfortunately not corresponding to the numerically simulated BER.
In this work, we show that the Generalized Gamma distribution (GGD) can efficiently model the turbulence induced self-channel fading as well as the crosstalk between OAM modes. The GGD gives a perfect fit for the weak, moderate, and strong turbulence regimes. Closed-form expressions for the outage probability, the average capacity, and the BER are derived for single-input single-output (SISO) transmission using only one OAM mode. Afterward, we consider a multiple-input multiple-output (MIMO) system using several OAM modes to provide diversity gain and hence improve the performance. Furthermore, space-time (ST) coding is proposed to achieve diversity and multiplexing gains, and the corresponding performance is analyzed.
The remainder of this paper is organized as follows. In Section II, we describe the OAM modes generation and propagation in turbulence, and then we present the GGD and its main statistical parameters. In Section III, we compare our proposed GGD model to numerical simulations of atmospheric turbulence fading and also to the Johnson $S_{B}$ and Exponential models. In Section IV, we derive closed-form expressions for performance metrics using the GGD model and compare the fitting goodness with numerical simulations. Afterward, spatial diversity using several OAM modes is studied and shown to give better performance than SISO links in section V. In section VI, ST coding is proposed to achieve full diversity and multiplexing gains and the theoretical error probability is derived. Finally, the conclusion is given in Section VII.

In this section, we evaluate the suitability of the proposed GGD to fit the simulation data. As we are working in an IM/DD setup, our Monte Carlo simulations of OAM propagation through the turbulent atmosphere are performed at wavelength $\lambda=850~$nm. The beam waist at the transmitter is $\omega_{0}=1.6$ cm. On the other hand, the optical receiver is assumed to be large enough to collect all the power on the OAM modes. This condition is satisfied when the receiver diameter is $d_{Rx}=2w_{z}\sqrt{m_{\text{max}}}$, where $w_{z}$ is the radius of the Gaussian beam and $m_{\text{max}}$ is the highest OAM topological charge used in the system. To create the desired OAM modes, SLMs having $512\times 512$ pixels are used. The propagation distance is set to $z=1~$km, the inner and outer scales of turbulence are set to $l_{0}=5~$mm and $L_{0}=20~$m, respectively. Atmospheric turbulence is emulated by placing $20$ random phase screens each $50~$m along the propagation path. Each phase screen is evaluated as the Fourier transform of a complex random distribution with zero mean and variance equal to $\left(\frac{2\pi}{N\Delta x}\right)^{2}\Phi(\kappa)$, where $N=512$ is the array length and $\Delta x=1.36~$mm is the grid spacing assumed to be equal in both dimensions $x$ and $y$. Propagation through the turbulent atmosphere is simulated using the commonly used split-step Fourier method [11].
Let $I_{mn}^{1},...,I_{mn}^{N}$ be the set of irradiance realizations for a particular $m$ and $n$. $N$ is the number of samples, where for each $I_{mn}$ we have simulated $N=10^{6}$ realizations. We use a maximum-likelihood (ML) estimation in order to estimate the parameters $a$, $b$, and $c$ of $I_{mn}$. The ML estimator computes the loglikelihood function given by [38]:

The estimated parameters $a$, $b$, and $c$ represent the values for which the previous function is maximized and are giving by [38, 39]:

where $\psi_{0}$ denotes the digamma function. In our simulations, we have used the ML implementation given in the wafo toolbox [40] .
To validate the goodness of the proposed GGD fit, we used the mean square error (MSE) test. The latter can give a measure on the accuracy of the GGD fitting to the simulated data. The MSE estimator is defined as:

where $\hat{F}^{k}(I_{mn})$ is the empirical distribution function of $I_{mn}$ and $F^{k}_{I_{mn}}(I_{mn})$ is the cumulative distribution function of $I_{mn}$ computed using the estimated GGD parameters of $I_{mn}$.
In Figs. 2, 3, and 4, we plot the distributions of $I_{mn}$ for $m\in\mathcal{M}=\left\{+1,+3,+5,+10\right\}$ and $n\in\mathcal{N}=\left\{m-1,m,m+1\right\}$ for different atmospheric turbulence regimes. Also for each distribution, a GGD fitting along with a Johnson $S_{B}$ or Exponential distribution fitting are plotted.
In Fig. 2, a weak turbulence regime given by $C^{2}_{n}=5\times 10^{-15}$ and $\sigma_{R}^{2}=0.2$ is considered. The center column of the figure shows the obtained results for the self-channel irradiance (i.e., the distribution of the power transmitted on OAM mode $m$ and detected on the same OAM mode $n=m$). The obtained distributions have a perfect correspondence with the distributions obtained in [18]. The histograms are left skewed and very close to the value of $1$. This is explained by the fact that almost all power transmitted on a specific OAM mode remains in this mode when atmospheric turbulence is weak. Furthermore, we notice that the GGD gives a perfect fit with the obtained histograms and very low values of the MSE were achieved. In addition, the left and right columns of Fig. 2 shows the distribution of the crosstalk of the adjacent OAM modes to the transmitted OAM mode $m$. We notice the Exponential shape of the crosstalk histograms which are perfectly fitted by the GGD, and achieving an MSE always less than $4\times 10^{-6}$.
In Fig. 3, a moderate to strong atmospheric turbulence regime with $C^{2}_{n}=7\times 10^{-14}$ and $\sigma_{R}^{2}=2.80$ is considered. From the center column of Fig. 3, we notice that the effect of turbulence is remarkable on the distributions of $I_{mm}$ which are no longer tight and close to 1. In this regime of turbulence, the GGD is also capable of fitting the simulated data histograms with an MSE less than $1.6\times 10^{-4}$. For OAM mode $m=+10$, the GGD provides a better fit than the Johnson $S_{B}$. In addition, the crosstalk histograms are perfectly fitted by the GGD and its fitting is also much better than the Exponential fit. Lastly, for the saturation regime, we have considered $C^{2}_{n}=3\times 10^{-13}$ and $\sigma_{R}^{2}=12$. As for the previous two cases, the GGD gives an excellent fit for all the obtained distributions and achieves a very low MSE. Moreover, for the crosstalk distributions, we notice from Fig. 4 that the GGD gives a better fit than the Exponential fit. Consequently, all the previous results indicate that the GGD is a very suitable candidate to model the irradiance distributions of OAM FSO transmission systems in the presence of atmospheric turbulence. Furthermore, we also notice from Figs. 2, 3, and 4 that the parameters $a$, $b$, and $c$ of the GGD depend on the OAM mode as well as the atmospheric turbulence strength. This situation is more complex than the classical FSO systems where the parameters of the irradiance of the Gaussian beam modeled as a Gamma-Gamma distribution could be directly derived from the inner and outer scale of the atmospheric turbulence [17]. Consequently, the computation of the parameters of the GGD in OAM FSO systems as a function of the OAM mode order and the turbulence strength is subject of future works.

In this section, we consider that only one OAM mode is transmitted and detected in a system operating in intensity modulation and direct detection. For notational convenience, we omit the subscript $m$ and refer to $I_{mm}$ by $I$. The instantaneous electrical signal-to-noise ratio (SNR) is defined as $\gamma=\left(\eta I\right)^{2}/N_{0}$. The average electrical SNR is defined as $\mu=\left(\eta\mathbb{E}\left\{I\right\}\right)^{2}/N_{0}$. After a simple power transformation of the variable $I$ in Eq. (8), the PDF of the instantaneous SNR can be written as:

The CDF of the SNR defined as: $F_{\gamma}\left(\gamma\right)=\int_{0}^{\infty}f_{\gamma}\left(\gamma\right)d\gamma$ can be obtained using the Meijer’s G function definition in [35, Eq. (2.9.1)]:

Now that we have defined the PDF and CDF of the SNR, we can proceed to the performance metrics derivation.

In this section, we propose to enhance the error performance of OAM FSO systems by using multiple OAM beams to transmit the same data. Typically, FSO systems have used multiple transmitters and receivers having apertures that are separated enough to bring diversity gain [44, 45]. However, this strategy comes at the cost of more hardware complexity and apertures spacing limitations which cannot always be satisfied. Spatial diversity can be used at the receiver in the form of a single-input multiple-output (SIMO) system as in [27]. At the receiver, combining techniques such as equal gain combining (EGC), selective combining (SC) or maximum ratio combining (MRC) can be used. Spatial diversity can also be set at the transmitter in the form of a multiple-input single-output (MISO) system [46]. In this case, beamforming can be used thanks to channel state information to improve the performance [47]. Moreover OAM mode selection can also be used to achieve better rates and lower BER [14, 27]. In the absence of channel state information, space-time coding can also be used as an alternative solution to increase the performance [14]. Spatial diversity can be used at both the transmitter and the receiver in the form of a MIMO system [46]. We consider a square MIMO system $\left(\mathcal{M}=\mathcal{N}\right)$ where $M=\left|\mathcal{M}\right|$ OAM modes are transmitted and received. By using a MRC at the receiver, the electrical SNR is given by [46, Eq. (14)]:

where the factor $1/M^{2}$ in the previous equation allows to compare the SISO and MIMO systems at the same transmitted optical power.
The derivation of the exact PDF of $\gamma_{\text{M}}$ is complex as it involves sums of GGDs that are non identically distributed. To overcome this issue and compare the performance of SISO and MIMO configurations, we approximate $I_{\text{M}}=\frac{1}{M}\sqrt{\sum\limits_{n\in\mathcal{N}}\left(\sum\limits_{m\in\mathcal{M}}I_{mn}\right)^{2}}$ by a GGD and estimate its parameters. Therefore, the PDF and CDF of $\gamma_{\text{M}}=\left(\eta I_{\text{M}}\right)^{2}/N_{0}$ can be directly given from Eq. (18) and Eq. (19), and the outage probability and average BER can be theoretically computed using Eq. (22) and Eq. (24).
Before performance evaluation of the MIMO system, it is important to choose the appropriate OAM modes in order to fully exploit the benefit of the spatial diversity. When an OAM mode of order $m$ is used for transmission, it is normal to detect the power conserved in that mode first. In fact, the power $I_{mm}$ is much higher than the power spread into other modes given by $I_{mn}$ as shown from Figs. 2, 3, and 4. When using spatial diversity, the detected signal on OAM mode $n$ contains the power transmitted from OAM mode $n$ and also crosstalk powers from other transmitted OAM modes as given by Eq. (5). Furthermore, the performance achieved by the spatial diversity is limited by the correlation between spatial paths [44]. Therefore, to achieve the maximum spatial diversity gain, the correlation between the self-channel irradiance $I_{mm}$ and the crosstalk $I_{mn}$ has to be minimal. To have an insight on the correlation between OAM modes, we propose to compute the correlation coefficient between $I_{mm}$ and $I_{mn}$ given by:

where $\text{cov}\left(\cdot,\cdot\right)$ is the covariance and $\sigma_{(\cdot)}$ is the standard deviation.

Consider that OAM mode $m=+1$ is transmitted, we compute the correlation coefficients between $I_{+1+1}$ and $I_{+1n}$ for $n\in\mathcal{N}=\left\{-10,-9,...,+9,+10\right\}$. In Fig. 8, we plot the correlation coefficients for the weak (resp. strong) turbulence regimes given by $C^{2}_{n}=5\times 10^{-15}$ and $\sigma_{R}^{2}=0.2$ (resp. $C^{2}_{n}=7\times 10^{-14}$ and $\sigma_{R}^{2}=2.8$). From Fig. 8, we notice that adjacent modes to the transmitted mode $m=+1$ have inverse correlation coefficient. The latter increases as the mode order gets further from $m=+1$. Therefore, closer OAM mode orders are preferred to achieve a maximum diversity [27]. In addition to the correlation, the power spread from the transmitted mode due to atmospheric turbulence is more important for adjacent OAM modes than further modes. Therefore, higher SNR gain will be provided by using modes with closer OAM order.
In Fig. 9, the performance of OAM FSO system using spatial diversity is presented for the moderate to strong turbulence regime given by $C^{2}_{n}=7\times 10^{-14}$ and $\sigma_{R}^{2}=2.80$. For both the outage probability and the average BER, we notice that the set $\mathcal{M}=\left\{+1,+2\right\}$ achieves better performance than the set $\mathcal{M}=\left\{+1,+3\right\}$ as expected. Furthermore, increasing the number of OAM modes enhances the performance as shown from the curve with the set $\mathcal{M}=\left\{+1,+2,+3\right\}$. In addition, we have also approximated $I_{\text{M}}$ with a GGD and estimated its parameters which allowed to compute the theoretical outage probability and average BER. The plotted curves based on the GGD approximation give an excellent match with the simulated data.

In the previous section, we have shown that spatial diversity can improve the error performance by sending multiple replicas of the same signal on different OAM modes. Nonetheless, the achieved diversity performance was obtained at the cost of a reduced multiplexing gain. To achieve full diversity and multiplexing gains, we consider space-time coding. ST coding was originally designed for wireless radio-frequency (RF) systems [48], and recently investigated for few-mode fibers [49, 50]. The Golden code [51] was designed for $2\times 2$ wireless MIMO systems and achieves a full diversity with the best coding gain for Rayleigh fading channels. The Alamouti code [52], is a half-rate orthogonal ST code that also achieves a full diversity and has the benefit of a very low decoding complexity that reduces to a simple channel inversion. ST coding was also demonstrated with heterodyne multiple apertures FSO systems [53]. However, due to the high cost of heterodyne implementation and complex signal processing required for ST decoding, it is not easily deployed in practice.

In addition, ST coding was also investigated for multiple apertures IM/DD FSO systems. In [54], a theoretical error probability was derived, and ST coding was shown to achieve full diversity. Furthermore, modified versions of known ST codes from RF were constructed and adapted to IM/DD systems. In [55, 56] modified versions of the Alamouti code were proposed. The new constructed Alamouti code preserves all the useful properties of the original code. The codeword matrix of the Alamouti code is given by [55]:

where the symbol $\sim s_{i}$ refers to the bit-wise not of the transmitted bit $s_{i}$. The main drawback of the Alamouti code is its half-rate code that only allows sending two bits $s_{1},s_{2}$ during two time slots and using two OAM modes.
Furthermore, Mroueh has proposed the Golden-Light (GL) code [57] which is an extended version of the Golden code that was adapted to IM/DD systems. The proposed code was shown to conserve an optimal and non-vanishing determinant. The GL codeword is given by [57]:

where $z^{*}_{1}=\frac{1}{2}\left(1+\sqrt{5}\right)$ and $z^{*}_{2}=2+\sqrt{5}$ are specific algebraic numbers related to the code construction. In addition to a full diversity, The GL code has a full-rate, which means that during two time slots and using two OAM modes, the GL code allows transmitting four bits $s_{1},s_{2},s_{3},s_{4}$.
To evaluate the performance of ST coding in OAM FSO systems affected by atmospheric turbulence, we consider that $M$ OAM modes are transmitted and detected. The resulting MIMO transmission system is given by:

where $\mathbf{X}\in\mathbb{R}^{M\times M}$ and $\mathbf{Y}\in\mathbb{R}^{M\times M}$ are the transmitted and received codewords. $\mathbf{N}\in\mathbb{R}^{M\times M}$ represents an additive white Gaussian noise with zero mean and a variance $N_{0}$ per complex dimension. $\mathbf{I}$ represents the channel matrix with entries $I_{mn}$. In the next section, we derive a theoretical expression of the error probability of ST coded OAM FSO systems.

In this paper, we have proposed a new statistical model for atmospheric turbulence fading in OAM FSO transmission systems. The proposed GGD can model the self-channel irradiance of OAM modes as well as the crosstalk between OAM modes for all atmospheric turbulence regimes. Hence, the proposed model can be used to overcome the computationally complex Monte Carlo simulations to simulate the propagation of OAM beams. Afterward, we have derived analytical expressions for the ergodic capacity, the outage probability, and the average BER. The obtained results show perfect match with simulations. Moreover, we have proposed to enhance the performance of OAM FSO systems by considering spatial diversity where several OAM modes are used to transmit the same data. From the obtained outage probability and average BER, significant performance improvements were achieved. Finally, ST coding was proposed to achieve full diversity and multiplexing gains. A theoretical derivation of the error probability was conducted, and numerical simulations of different ST codes were presented.