In this paper, orbital angular momentum (OAM) modes transmission in the presence of atmosphere turbulence is studied via a coupled mode theory. The Laguerre-Gauss (LG) beams with OAM topological charges are emitted into free space and undergo interactions due to the random index variations in the atmosphere. The coupling between the LG beams can be characterized by a set of coupled average power equation, which resembles the Marcuse' coupled power equation (CPE) originally proposed for the optical waveguides. The coupling coefficients and the modal radiation losses for the equation can be evaluated analytically. The accurate solution and the first order approximate solution to the CPE match the published data and the Mont-Carlos simulation results with good accuracy. The CPE and its approximate analytical solution can work as powerful tools for the analysis of the OAM beam evolution with the presence of the atmosphere turbulence.
© 2015 Optical Society of America
Optical Orbital angular momentum (OAM) multiplexing is one of the most promising techniques to realize space division multiplexing (SDM) in the optical communication systems [1–4 ]. It uses OAM modes with different topological charges as the information carriers and can greatly improve the system capacity.
In the ideal radial symmetrical environment, OAM mode will maintain its topological charge. The orthogonality between the OAM modes enables simple multiplexing/de-multiplexing among the channels. Since most of the current free space OAM transmission experiments have been demonstrated within a short distance in free space [1–3 ], the orthogonality maintains perfectly and therefore, crosstalk between the OAM modes during transmission is negligible. The main penalty comes from the OAM multiplexing/de-multiplexing components.
However, with the increase of the transmission distance, the refractive index variations caused by inhomogeneous temperature and pressure distributions in the atmosphere become non-negligible. It is known that such atmosphere turbulence will bring crosstalk during OAM modes propagation and many studies have been conducted to quantify its influence [5–9 ]. C. Paterson  and Glenn A. Tyler  studied the impact of random index perturbation analytically under the weak turbulence assumption, which models the impact of the index variation by introducing phase change on the beam. B. Rodenburg  verified the phase perturbation statistical theory via experimental demonstrations. Y. Zhang et. al. presented similar results for Whittaker-Gaussian (WG) beams propagating in the turbulent atmosphere . As pointed out in these pioneering studies , the phase perturbation model might not be valid in long distance point-to-point communications , since the random scattering becomes quite complex in such an environment and the phase perturbation model becomes inaccurate. To study the phenomenon in a more general sense, J. A. Anguita used split step Fourier method (SSFM) to simulate the propagation process . It can accurately provide the characterization for the long distance beam propagation in the media with large index variations. However, the method is a little time consuming because the random nature of the atmosphere turbulence requires Mont-Carlo simulations, which are composed of numerous rounds of simulations to find the average value. Therefore, despite the extensive analytical and numerical studies on this topic [5–9 ], there is a great demand for an efficient analytical tool that could account for the strong atmosphere turbulence during the relatively long distance propagation in free space.
In this paper, we present a set of coupled power equation (CPE) for the OAM modes transmission in the presence of the atmosphere turbulence, which is similar to the Marcuse' coupled power theory [10–12 ] for the optical waveguides in the presence of index perturbations. It is presented that the coupling coefficients and the mode radiation losses in the CPE can be calculated analytically. The CPE is a set of ordinary differential equations (ODEs), and it can be solved easily by the numerical integration techniques such as the Runge-Kutta method. Since the numerical integration is accomplished only once, the computational cost is significantly reduced in comparison with the Mont Carlo simulations. Furthermore, since the LG modes vary little when the propagation distance does not exceed the Rayleigh range, the coupling coefficients and the radiation losses can be regarded constants. Therefore, the first order approximate solution to the CPE, which assumes constant coupling coefficients and losses, can provide a good approximation. Since the first order approximation can be calculated analytically, it becomes a powerful tool in combination with the CPE for the analysis of the long distance OAM modes transmission in the presence of the atmosphere turbulence. The results from the proposed CPE are verified by the comparison with the published data in  and the Mont Carlo simulation results.
2. The coupled power equations
To derive the CPE similar to the Marcuse' theory [10–12 ], we start with the derivation of the coupled mode equations for the OAM modes.
It is generally accepted that the scalar wave equation will be sufficient to describe the beam propagation in the atmosphere [5–9,13 ]
where k is the wave number in the media. In the presence of the turbulence, the wave number k can be expressed by 
where k 0 is the free space wave number, δn the refractive index turbulence. To facilitate the analysis, the refractive index variation can be modified as a zero-mean random process in the three dimensional space by subtracting the average value of δn, i. e.
Therefore, one has
Without the atmosphere turbulence, Laguerre-Gauss (LG) beams are the exact solution to the Helmholtz equation in free space. And the hence, if LG beams are used for OAM multiplexing and the atmosphere turbulence are not considered, the field can be expressed by
where m stands for the mth order LG beam. The LG beams u m fulfill the scalar Helmholtz equation in free space, i. e.
and they have the form of 
where w 0 is the initial beam waist, λ the wavelength, Lpm the LG polynomial. In our case we have p = 0, Lpm = 1, and the OAM modes are ring shape beams except for the m = 0 case .
It should be noted that these LG modes have been normalized and they are orthogonal to each other
where * denotes the complex conjugate. It should be noted that the normalization expression adopted in  and this work is different from the normalization expression in [10–12 ]. The normalization expression in [10–12 ] was proposed for the traditional coupled mode theory in optical waveguides, and an additional factor appears, where βm is the propagation constant of the mth waveguide mode. Since in the case of LG modes, the propagation constants are very close to each other, the factor mentioned above becomes a constant for all of the LG modes. Therefore, the normalization used in Eq. (9) does not impact the final results as compared with the traditional coupled mode theory.
In the presence of turbulence, the coefficient am in Eq. (6) becomes z dependent. By assuming the waves are solely formed by the LG beams, we can substitute Eq. (6-7) into Eq. (5). Neglecting the higher order derivatives with respect to z for am and um, one has
The OAM mode um varies very slowly while the phase term changes very fast with respect to the propagation distance z, and hence we have . Therefore, can be neglected because it is much smaller than and we have
It should be noted that although the waves are mainly composed of LG beams, it consists of a portion of other beams, since the LG beam does not form a complete orthogonal basis. The contribution of these other beams can be accounted for as radiation losses of the LG beams. In the derivation, the radiation losses are ignored and added in the final CPE as Marcuse did in his papers [10–12 ].
The turbulence induced index fluctuation has the correlation function as
The Fourier transform of the correlation function R should be the Kolmogorov spectrum ,
where FT denotes the Fourier transform, kx ky and kz the spatial frequency in the x, y, and z directions, l0 and L0 the constants related the spectrum distribution , Cn 2 the refractive-index structure parameter. It should be noted that the expression in this work differs from the expression in  by a factor of (2π)3, because the general definition of the Fourier transform  is used in this work.
Since the propagation distance z is much larger than the beam size in the x and y directions, the random process can be regarded as two independent processes in the x-y direction and in the z direction. Therefore, the correlation function can be approximately expressed by
when FT −1 denotes the inverse Fourier transform. The power of each LG beam can be expressed by
With the coupling coefficient as
Equation (16) has not considered any coupling between the LG beams and the radiation modes. By taking radiation modes into account, the exact coupled power equation should be
where αm stands for the radiation loss of each LG beam.
The coupling coefficient can be evaluated analytically as follows. Equation (17) can be rewritten as
where denotes two dimensional convolution. Equation (19) can be calculated by the property of Fourier transform
The radiation loss αm can be evaluated by summing up (or integrating) the coupling coefficient between the LG beams and the radiation modes.
Since it is quite hard to obtain the expression of these radiation modes, we can use the plane waves to approximate the radiation modes as it is done in the case of optical waveguides . Therefore, the loss can be calculated as
Strictly speaking, the plane approximation for the radiation modes is only valid under the paraxial approximation during the beam propagation, i. e. when kx and ky and small. Since the LG modes are 'large' in comparison with the wavelength, the value of kx and ky are very small in comparison with k, the paraxial approximation is valid. For example, the largest value of kx and ky is determined by the grid size in the numerical simulation, and the grid size is usually in the order of 10−4 m, which is already small enough to depict the variation of the LG beams. Hence, the largest value of kx and ky is two orders of magnitude smaller than the propagation constant k. Therefore, Eq. (21) can be regarded valid for all kx and ky, and the integration intervals of kx and ky in Eq. (21) can be considered as (-∞, + ∞).
It should be noted that the radiation loss calculated in Eq. (21) has included the loss contributed by the LG modes (see Eq. (18)), so the true radiation loss αm should subtract these terms, i. e., from Eq. (21).
From Eq. (21), one can draw a very interesting conclusion. The loss of each LG mode (including the contribution from the other LG modes) is only impacted by the atmosphere turbulence. Therefore, the parameters of the beam will only affect the coupling between the LG modes, not the losses of them.
The CPE, i. e. Equation (18), can be solved via numerical integration methods such as the Runge Kutta method. Its solution can be conveniently expressed by rewriting it in the matrix form,
And its solution is
Because the coupling coefficients almost do not vary much during the propagation if the distance does not exceed the Rayleigh distance ZR, Eq. (23) can be simplified by assuming
Therefore, the first order approximate solution to Eq. (22) can be represented in the analytical form
3. Results and discussions
In order to verify the proposed coupled power model, a numerical example in  is used. The turbulence parameters L 0 = 20m, and l 0 = 5mm. The signal wavelength is 850nm and the propagation distance is 1000m. These parameters are exactly the same as the ones in .
The numerical parameters are as follows: the simulation is conducted in the area of 140cmX140cm with the grid of 1024X1024 elements, which is enlarged in comparison with  (In , it is 70cmX70cm with the grid of 512X512 elements). In the z direction, random phase plates are inserted every 50m with randomly generated phases . Split-step Fourier method (SSFM) is used to simulate the propagation . For the details of the SSFM, the readers are encouraged to refer to .
The validity of the SSFM Mont Carlo simulation results is firstly checked by calculating the intensity and phase distribution of the OAM modes when the input is purely OAM state + 5. Cn 2 is 10−14 m-2/3. The initial beam waist is 3cm and the propagation distance is 1000m. The intensity and phase distribution of the output field during one of the random realization is shown in Fig. 1 .
It can be seen from the figure that the corresponding intensity/phase distribution matches the results in , which were also obtained by the SSFM.
Next, the scintillation index is evaluated. The initial waist is changed to 1.6cm. The refractive-index structure parameter Cn 2 is 10−14 m-2/3, and the propagation distance varies from 100m to 1000m (the Rytov variance is about 0.04 to 0.4). The input beam is assumed to be OAM state 0. Its self scintillation index is evaluated and plotted in Fig. 2 . When the Rytov variance increases, the scintillation index increases accordingly. This is in accordance with the theoretical and experimental results . When z = 1000m, the scintillation index is about 0.14. This is a reasonable value for the Gaussian beam as indicated by . In , when the Rytov variance is 0.4, the measured scintillation index of the narrow collimated Gaussian beam is about 0.2.
The numerical model can also be used to investigate the correlation between the modes, which is of great interest to investigate the MIMO channel capacity. In this case, the initial waist is 1.6cm. The refractive-index structure parameter Cn 2 is 10−14 m-2/3, and the propagation distance is 1000m. OAM modes −9 to + 9 are emitted at the input and detected at the output. The signal on each mode is assumed to be independent of each other and with equal power, which is reasonable during information transmission. The correlations between OAM state 0 and other OAM states at the output are plotted in Fig. 3 . It can be seen from the figure that correlation is the strongest between the neighboring channels. The correlation varies from 0.3 to 0.01, which are reasonable values for the channel correlation in the presence of atmosphere turbulence .
The above calculations indicate the numerical model is quite valid to simulate the atmosphere turbulence. With the above simulation parameters, the final power transfer matrix is calculated. The power transfer matrix M is defined as
It can either be calculated by Eqs. (23) and (25) , or by the Mont-Carlo simulations.
The first verification of the proposed coupled mode power model is accomplished by comparing the detailed data in  and the theoretical predictions of our models. The beams in  are LG beams with the initial beam waist as 1.6cm. The refractive-index structure parameter Cn 2 is 10−14 m-2/3. (It is worth noting that our results based on the Mont-Carlo simulations are very close to the ones in  if the simulation area and the grid elements number are reduced to those in . The data in  is also used for comparison to verify the accuracy of the proposed CPE model.)
The diagonal terms of the power transfer matrix stand for the channel efficiency, i.e. the amount of power which is left after propagation for each channel. In Fig. 4 , the channel transmission efficiencies are presented. The results are calculated by Eq. (23) and Eq. (25), and obtained from . As the mode order increases, less power can be maintained in the transmitted LG beams after propagation. The three curves matches quite well, demonstrating the accuracy of our formulas. The relative error between the data in  and our solutions is less than 8% for the “accurate” solution and less than 13% for the approximate solution. If our Mont Carlo simulation results with the enlarged simulation area are used for comparison, the relative error is less than 4% for the “accurate” solution and less than 16% for the approximate solution.
In Fig. 5 and Fig. 6 , the crosstalk coefficients between the channels, i. e. the off-diagonal terms in the transfer matrix M, are presented. Figure 5 gives the crosstalk from OAM state m to OAM state m-1, which is the crosstalk between neighboring channels, while Fig. 6 demonstrates the crosstalk from OAM state m to OAM state m-2 respectively. Again, the theoretical calculation, the data by our Mont Carlo simulations and the data form  match well.
To verify the robustness of the proposed CPE, different simulation parameters, such as the different atmosphere turbulence parameters, different beam initial waist, and different propagation distances are used.
Figure 7 and Fig. 8 show the channel efficiency under different atmosphere turbulence. The refractive-index structure parameter Cn 2 is 10−15 m-2/3 in Fig. 4 and 10-16 m-2/3 in Fig. 8 respectively. The beam initial waist is 1.6cm. It can be seen from the figures that the accuracy is quite promising both in the media with moderate (Cn 2 = 10−15 m-2/3) and weak (Cn 2 = 10−16 m-2/3) atmosphere turbulence. In case of moderate atmosphere turbulence (Fig. 7), the relative error is less than 5% for the “accurate” solution and less than 10% for the approximate solution. The relative error further reduces to 0.06% and 0.6% in the later case of weak atmosphere turbulence (Fig. 8).
Figure 9 demonstrates the channel efficiency when the beam initial waist is 3cm. Since the beam waist is larger in this case, the numerical simulation area is enlarged as 3m X 3m. The reason that the simulation area will impact the final results is that the Mont-Carlo simulation has used SSFM to simulate the beam propagation, which has assumed a periodical field distribution when applying the fast Fourier transform (FFT). It means that the plane waves keep interacting with each other within the simulation region. Such an interaction might introduce a “reverse coupling” effect, i.e. coupling of the power from the radiation modes back to the LG beams, which will not happen in real LG beam propagation, because the radiation modes will dissipate quickly and are impossible to be coupled back to the LG beams. The impact has been demonstrated in the pervious figures by comparing our simulation results in the enlarged area and the results in . It can also be seen from Fig. 9, which shows the results from the 3m X 3m simulation area are much closer to our theoretical predictions. Also, it can be seen from the figure that the CPE provides very good estimation for the channel efficiency. The relative error is less than 15% both for the “accurate” solution and the approximate solution. Further increase of simulation area can reduce the error, but computation cost becomes unbearable.
Figure 10 demonstrates the accuracy of the channel efficiency when the propagation distance varies for 100m to 2000m. For simplicity, only the efficiency of channel 0 is considered. The initial beam waist is 1.6cm and Cn 2 = 10−14 m-2/3. It can be calculated from Eq. (8) that the Rayleigh range is 946m. The Mont Carlo simulation area is still 70cmX70cm in order to be in accordance with . From the figure, it can be seen that as the propagation distance increases, the relative error between the theory and the Mont-Carlo simulation increases. But the relative error remains at a low level until the propagation distance exceeds the Rayleigh range both for the “accurate” solution and the approximate solution. When the propagation distance goes far beyond the Rayleigh range, the “accurate” solution is still accurate with the relative error of 10% in comparison with the Mont Carlo simulations. The increase of the error is mainly due to the computational error caused by beam broadening, which can be reduced by increasing the simulation area during the Mont Carlo simulation.
It can be seen from Figs. 4-10 that both the approximate solution to the CPE and the “accurate solution” to it have quite good accuracy with respect to the Mont Carlo simulation results within the Rayleigh range. This can be explained from the fact that coupling coefficients vary quite little during the propagation if the distance does not exceed the Rayleigh distance ZR. Therefore, from a practical point of view, the first order approximate solution is accurate enough to be used to predict the OAM state propagation characteristics in most of the applications. The discrepancy between the theory (the “accurate” solution) and the Mont Carlo simulations can be mostly attributed to the “reverse coupling effect” caused by the limited simulation area in the Mont Carlo simulations. Although the numerical parameters in  (which are also adopted in some of our Mont Carlo simulations) have already been optimized, the “reverse coupling” effect still plays a role especially when the random index perturbation is large and the propagation distance becomes long. It becomes even more significant when the propagation distance exceeds the Rayleigh range, because the beam will spread significantly and result in stronger “reverse coupling” in this case. Generally speaking, the “reverse coupling” effect is determined by the beam size, the strength of the turbulence, and the propagation distance. Larger simulation area will reduce the “reverse coupling” effect, but it will bring tremendous computational efforts.
In summary, we have presented an efficient theory to characterize the impact of the atmosphere turbulence during the OAM modes propagation in free space. Coupled average power equation is derived based on the coupled mode theory along with the analytical formulation of the coupling coefficients and modal losses. The CPE can be either solved numerically, which also saves computational cost by avoiding Mont Carlo simulations, or be approximated by its first order approximate solution. The results from the CPE are verified by the published data and the Mont Carlo simulation results. With the approximate analytical solution in combination with the coupled average power equation, it is possible to predict the evolution behavior of the OAM modes in a very efficient way.
This work is partially supported by the National science foundation of China (Grant No. 61201068) the Fundamental Research Funds for the Central Universities of China.
References and links
1. J. Wang, J. Yang, I. M. Fazal, N. Ahmed, Y. Yan, H. Huang, Y. Ren, Y. Yue, S. Dolinar, M. Tur, and A. E. Willner, “Terabit free-space data transmission employing orbital angular momentum multiplexing,” Nat. Photonics 6(7), 488–496 (2012). [CrossRef]
2. H. Huang, G. Xie, Y. Yan, N. Ahmed, Y. Ren, Y. Yue, D. Rogawski, M. Tur, B. Erkmen, K. Birnbaum, S. Dolinar, M. Lavery, M. Padgett, and A. E. Willner, “100 Tbit/s Free-Space Data Link using Orbital Angular Momentum Mode Division Multiplexing Combined with Wavelength Division Multiplexing,” in Optical Fiber Communication Conference/National Fiber Optic Engineers Conference 2013, OSA Technical Digest (online) (Optical Society of America, 2013), paper OTh4G.5. [CrossRef]
3. B. Guan, R. P. Scott, C. Qin, N. K. Fontaine, T. Su, C. Ferrari, M. Cappuzzo, F. Klemens, B. Keller, M. Earnshaw, and S. J. B. Yoo, “Free-space coherent optical communication with orbital angular, momentum multiplexing/demultiplexing using a hybrid 3D photonic integrated circuit,” Opt. Express 22(1), 145–156 (2014), http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-22-1-145. [CrossRef] [PubMed]
4. N. Bozinovic, P. Kristensen, and S. Ramachandran, “Long-range fiber-transmission of photons with orbital angular momentum,” in Proceedings of the Conference on Lasers and Electro-Optics (Optical Society of America, 2011), paper CTuB1. [CrossRef]
6. J. A. Anguita, M. A. Neifeld, and B. V. Vasic, “Turbulence-induced channel crosstalk in an orbital angular momentum-multiplexed free-space optical link,” Appl. Opt. 47(13), 2414–2429 (2008). [CrossRef] [PubMed]
7. G. A. Tyler and R. W. Boyd, “Influence of atmospheric turbulence on the propagation of quantum states of light carrying orbital angular momentum,” Opt. Lett. 34(2), 142–144 (2009). [CrossRef] [PubMed]
8. B. Rodenburg, M. P. J. Lavery, M. Malik, M. N. O’Sullivan, M. Mirhosseini, D. J. Robertson, M. Padgett, and R. W. Boyd, “Influence of atmospheric turbulence on states of light carrying orbital angular momentum,” Opt. Lett. 37(17), 3735–3737 (2012). [CrossRef] [PubMed]
9. Y. Zhang, M. Cheng, Y. Zhu, J. Gao, W. Dan, Z. Hu, and F. Zhao, “Influence of atmospheric turbulence on the transmission of orbital angular momentum for Whittaker-Gaussian laser beams,” Opt. Express 22(18), 22101–22110 (2014). [CrossRef] [PubMed]
10. D. Marcuse, “Mode conversion caused by surface imperfections of a dielectric slab waveguide,” Bell Syst. Tech. J. 48(10), 3187–3215 (1969). [CrossRef]
11. D. Marcuse, “Derivation of coupled power equations,” Bell Syst. Tech. J. 51(1), 229–237 (1972). [CrossRef]
12. D. Marcuse, “Power distribution and radiation losses in multimode dielectric waveguides,” Bell Syst. Tech. J. 51(2), 429–454 (1972). [CrossRef]
14. A. V. Oppenheim, A. S. Willsky, and S. Hamid, Signal and Systems, 2nd edition (Prentice Hall, 1997).
15. J. Zhou and P. Gallion, “A comprehensive analytical model to characterize randomness in optical waveguides,” J. Lightwave Technol. (submitted to).
16. D. Lenz, D. Erni, and W. Bächtold, “Modal power loss coefficients for highly overmoded rectangular dielectric waveguides based on free space modes,” Opt. Express 12(6), 1150–1156 (2004), https://www.osapublishing.org/oe/abstract.cfm?uri=oe-12-6-1150. [CrossRef] [PubMed]
17. V. P. Aksenov, “Scintillation index of a laser beam having the initial orbital angular momentum in the turbulent atmosphere,” in Imaging and Applied Optics 2015, OSA Technical Digest (online) (Optical Society of America, 2015), paper PM2C.5.
18. J. A. Anguita, M. A. Neifeld, and B. V. Vasic, “Spatial correlation and irradiance statistics in a multiple-beam terrestrial free-space optical communication link,” Appl. Opt. 46(26), 6561–6571 (2007). [CrossRef] [PubMed]
19. N. Baddour, “Two-Dimensional Fourier Transforms in Polar Coordinates,” Adv. Imaging Electron Phys. 165, 1–45 (2011). [CrossRef]