1. Introduction

Light that carries orbital angular momentum (OAM), i.e. optical vortices or twisted light, are a family of beams that have helical or twisted wavefronts. These beams carry OAM dependent upon the spatial field profile along with spin angular momentum. In contrast to spin, where the photon is limited to two states, OAM allows for, in principle, an infinite set of spatially orthogonal states [1]. This feature can be utilized for spatial division multiplexing (SDM) or quantum key distribution (QKD) to significantly increase the bandwidth of free space optical (FSO) communication systems [2–5] and remote sensing [6,7]. For such systems, propagation through atmospheric turbulence can significantly limit performance [8–10].

Atmospheric turbulence creates a temporally and spatially varying refractive index which can distort optical fields. In the context of OAM light, turbulence causes mixing between OAM states and thus can lead to crosstalk between SDM communication channels. The mixing process is illustrated Fig. 1 for a set of two channels (p₁,m₁) and (p₂,m₂), each in a different OAM mode, operating in parallel within the same beam profile. The total power in each mode (p,m) consists of the remaining pump power, denoted as P­­_{(p,m)→(p,m)}, plus the power transferred from all other modes (p’,m’) ≠ (p,m) due to turbulence, denoted as P­­_{(p’,m’)→(p,m)}. The ability to estimate the modal power spectrum of OAM beams after propagation through turbulence would aid in the design and characterization of free-space OAM systems.

 

Fig. 1. OAM mode mixing in turbulence. The net power in each mode is the remaining pump power plus the power injected by all other modes.

Download Full Size | PDF

In this work we consider Laguerre-Gauss (LG) modes, which are a class of OAM beams that satisfy the paraxial wave equation and are described in cylindrical coordinates (s, ϕ, z) by the normalized complex field envelope

The probability of a photon initially in the (p,m) mode evolving to carry m'ħ OAM after propagating through turbulence has been calculated [10]. In the context of LG beams and Fig. 1, this is equivalent to summing the power transfer across all radial modes, i.e. $\sum\nolimits_{p^{\prime} = 0}^\infty {{P_{(p,m) \to (p^{\prime},m^{\prime})}}}$. Theoretical [11–13], numerical [8,14], and experimental [15] studies have explored the spectral evolution of LG beams in the (0,m) basis subset under a variety of turbulent conditions. A more complete understanding of the full mode power spectrum, i.e. both radial p and azimuthal m, is desirable, and would allow for more accurate system performance predictions. In particular, LG mode sorters that utilize the radial mode basis as an additional degree of freedom for FSO communications systems have been developed [16–17].

We present a theoretical and numerical description of the full mode power spectrum of LG beams in Kolmogorov turbulence under the assumption of a collimated beam. The theoretical results are validated against split-step simulations and are shown to be accurate even in the strong turbulence regime. We also investigate the impact of the so-called quadratic approximation for the wave structure function on the accuracy of the final result. In addition, we consider the dependence of the input mode power on the RMS spot size of the beam.

2. Full modal power spectrum of LG beams in turbulence

Since LG modes constitute a complete and orthogonal set, one can expand a beam after propagating through turbulence as $U_{p,m}^{(t)}({\mathbf r}) = \sum\nolimits_{p^{\prime},m^{\prime}} {{c_{p^{\prime},m^{\prime}}}(z)\,\,U_{p^{\prime},m^{\prime}}^{(d)}({\mathbf r})}$, where the superscript (t) indicates the field after propagating through turbulence, (d) indicates diffraction-only propagation, and $|{c_{p^{\prime},m^{\prime}}}(z)\,{|^2}$ is the power in the (p’,m’) mode. As the LG modes in Eq. (1) are normalized, the fractional power transferred from an input LG mode (p,m) to mode (p’,m’) after propagating through turbulence is given by $\left\langle {|{c_{p^{\prime},m^{\prime}}}(z)\,{|^2}} \right\rangle$, i.e.

In the case of beams which are well-collimated, i.e. the total propagation distance is much less than the Rayleigh length z_R = kw²/2, one can discard the transverse Laplacian to yield $[{{\partial_z} - ik\,\delta n({\mathbf r})} ]U_{p,m}^{(t)}({\mathbf r})\, = 0\,.$ This equation can be expressed in the form

Defining x = 2s²/w², ε = w/2^1/2ρ₀ ≈$0.89w{({C_n^\textrm{2}{k^2}z} )^{3/5}}$, and ζ = ϕ’ - ϕ Eq. (7) can be written as

The integral in Eq. (8) is oscillatory and difficult to evaluate, especially as the azimuthal mode separation increases. Making the quadratic approximation, i.e. ${\varepsilon ^{5/3}}{(x + x^{\prime} - 2{(xx^{\prime})^{1/2}}\cos \zeta )^{5/6}} \approx {\varepsilon ^2}(x + x^{\prime} - 2{(xx^{\prime})^{1/2}}\cos \zeta )$, the integral over ζ can be carried out using $\exp [\alpha \cos \zeta ] = \sum\nolimits_{\nu ={-} \infty }^\infty {{I _\nu }[\alpha ]\exp [i\nu \zeta ]}$ where I_ν is the modified Bessel function of the first kind. Equation (8) becomes

3. Validation via split-step simulation

Simulations based on the Fourier split-step method [21–23] were performed to validate Eqs. (8) and (9). Our simulations do not employ the quadratic approximation. As an example, we consider a collimated beam initially in the (p = 2, m = 2) LG mode with wavelength λ = 875 nm and zero-order spot size w = 33.25 mm (RMS spot size 88 mm) propagating 500 m in air. With these parameters, the total distance traveled corresponds to 0.13 z_R, and thus the beam remains well collimated. The structure constant is taken to be constant throughout the medium, and multiple simulations are performed with $C_n^2$ values ranging from 5.1 × 10⁻¹⁵ to 6.4 × 10⁻¹³ m^-2/3, i.e. weak to strong turbulence. Since an individual phase screen can only accurately represent weak turbulence, 60 evenly spaced screens were used, each constructed using the randomized spectral sampling technique described in Ref. [24]. For all simulations a 2048 × 2048 grid is used with a grid spacing of 0.35 mm. For each simulation, the beam is decomposed into LG modes as shown in Eq. (2) every 50 m. The mode basis included p’ = 0 to 6 and m’ = -2 to 6 for a total of 63 modes. The number of iterations used for the ensembles ranged from 800 (for weak turbulence) to 8400 (for strong turbulence).

The results of the simulations are compared to the theoretical relationships in Eqs. (8) and (9) in Figs. 2–4, where ${P_{(2,2) \to (p^{\prime},m^{\prime})}}(z)/{P_{(2,2)}}(0)$ is plotted as a function of ε. The results are separated out for individual values of p’ and for m’ ≤ m and m’ ≥ m. The filled circles are the simulation results, the dashed lines Eq. (8) (without the quadratic approximation), and the solid lines Eq. (9) (with quadratic approximation). Note that each figure does not necessarily have the same y-axis bounds. The parameter ε does not have a one-to-one mapping to the Rytov variance $\sigma _R^2$, but for our parameter range ε = 0.1 to 10 corresponds to $\sigma _R^2$ = 3 × 10⁻⁴ to 7 and thus represents very weak to strong turbulence.

 

Fig. 2. The fractional power transferred from LG mode (2,2) to modes (p’,m’) as a function of ε = w/2^1/2ρ₀, for p’ = 0 to 2 and m’ = -2 to 6. The filled circles represent the split-step simulation results, the dashed lines Eq. (8) (without quadratic approximation), and the solid lines Eq. (9) (with quadratic approximation). The left figures (a,c,e) correspond to m’ ≤ m, and the right (b,d,f) to m’ ≥ m.

Download Full Size | PDF

 

Fig. 3. The fractional power transferred from LG mode (2,2) to modes (p’,m’) as a function of ε = w/2^1/2ρ₀, for p’ = 3 to 4 and m’ = -2 to 6. The filled circles represent the split-step simulation results, the dashed lines Eq. (8) (without quadratic approximation), and the solid lines Eq. (9) (with quadratic approximation). The left figures (a,c) correspond to m’ ≤ m, and the right (b,d) to m’ ≥ m.

Download Full Size | PDF

 

Fig. 4. The fractional power transferred from LG mode (2,2) to modes (p’,m’) as a function of ε = w/2^1/2ρ₀, for p’ = 5 to 6 and m’ = -2 to 6. The filled circles represent the split-step simulation results, the dashed lines Eq. (8) (without quadratic approximation), and the solid lines Eq. (9) (with quadratic approximation). The left figures (a,c) correspond to m’ ≤ m, and the right (b,d) to m’ ≥ m.

Download Full Size | PDF

We find that the theoretical result without the quadratic approximation (Eq. (8)) is accurate for all turbulence strengths and decomposition modes tested. The theoretical result with the quadratic approximation (Eq. (9)) is accurate for all decomposition modes tested and ε ≳ 1 but can lead to an underprediction of mode power when ε < 1. For our parameter sets, ε ≈ 1 corresponds to $\sigma _R^2$ ≈ 0.02 to 0.1. Some key exceptions are the input mode (2,2), the adjacent azimuthal modes (2,1) and (2,3), and (1,3) and (3,1), for which the quadratic approximation appears accurate throughout the tested turbulence range. Note that these modes have the most power at ε = 0.1 out of the modes tested: all are above 0.01×P_(2,2)(0). The accuracy of the quadratic approximation in the moderate and strong turbulence regime is valuable, as Eq. (9) is less oscillatory and thus more readily evaluated numerically than Eq. (8).

There are simulation results that appear to produce small discontinuities near ε = 4 to 5 and are somewhat divergent from the theory curves beyond ε = 5, in particular (3,6), (4,6), and (5,6). The discontinuity can be attributed to the merging of the different ensemble datasets. In addition, the mode power being a function of ε assumes a perfectly collimated beam, which is not strictly true throughout the 500 m of propagation for the simulations. Indeed, the data point at ε = 4 corresponds to the last sample of the $C_n^2$ = 1.3 × 10⁻¹³ m^-2/3 dataset, while ε = 5 corresponds to the first plotted sample of the 6.4 × 10⁻¹³ m^-2/3 dataset.

4. Impact of RMS spot size on power remaining in the input mode

In Ref. [13], we found that the fractional power remaining in an input (0,m) mode is largely dependent on the ratio of the RMS spot size to the coherence diameter, regardless of turbulence strength. The same feature was also found in Ref. [10], but for the likelihood that a photon initially in the (p,m) mode maintains mħ OAM (regardless of the final radial mode). To determine whether this behavior extends to arbitrary LG mode (p,m), we performed additional simulations for beams initially in the (0,0), (3,-3), and (4,4) states with RMS spot size ${(2p + |m|+ 1)^{1/2}}w$ = 88 mm. Remaining simulation parameters are the same as described in Section 3. The results are shown in Fig. 5 alongside the predictions based on Eq. (9). We find that the power remaining in the input (p,m) mode for collimated beams does indeed largely depend on the ratio of the RMS spot size to the coherence diameter, even into the strong turbulence regime. Thus, the topological charge only impacts the power remaining in the input mode insofar as it relates to the RMS spot size. Note that the curves based on Eq. (9) for the (2,2), (3,-3), and (4,4) are very similar and are difficult to individually distinguish in the figure.

 

Fig. 5. The fractional power remaining in various input LG mode as a function of ${({2p + |m|+ 1} )^{1/2}}\varepsilon $, where ε = w/2^1/2ρ₀ and ${({2p + |m|+ 1} )^{1/2}}w$ is the RMS spot size of the (p,m) LG mode. The different tick marks represent the split-step simulation results, and the solid lines Eq. (9).

Download Full Size | PDF

We can also expect that beams with longer wavelengths will be more robust against turbulent mixing since ε ∝ k^6/5. Though we note that longer wavelengths also lead to greater diffractive spreading angles (i.e. shorter Rayleigh lengths) and thus, for a given input RMS spot size, the collimated beam requirement of our analysis leads to a shorter range of validity for Eqs. (8) through (10).

5. Conclusions

In summary, we presented a theoretical model that quantifies the full mode power spectrum of collimated LG beams in Kolmogorov turbulence. The model is validated against split-step simulations and found to be accurate into the deep turbulence regime. We also find that use of the quadratic approximation introduces minimal error when ε $\mathrm{\ \mathbin{\lower.3ex\hbox{$\buildrel> \over {\smash{\scriptstyle\sim}\vphantom{_x}}$}}\ }$1, i.e. above $\sigma _R^2$ ≈ 0.05 for our parameter sets, but can lead to underprediction of mode power in weaker turbulence. We also found that the power remaining in the input LG mode largely depends on the ratio of the RMS spot size to the coherence diameter, regardless of turbulence strength. Our results should prove useful in the design and characterization of free-space OAM systems.