An Effective Way for Simulating Oceanic Turbulence Channel on the Beam Carrying Orbital Angular Momentum

In this paper, we present an effective way for simulating oceanic turbulence channel on the beam carrying orbital angular momentum (OAM). The influence caused by oceanic turbulence channel on the phase and intensity of the propagation beam is equivalent to that the beam passing through several individual phase screens generated by power spectrum inversion method at regular intervals. A modified subharmonic compensation method is then further balance the phase screen for the losses of lower frequency components in the power spectrum inversion method. The feasibility is verified by the theoretical phase structure function and the propagation characteristics of an OAM beam in underwater environment. The results show that the phase structure function and the propagation characteristics of the OAM beam evaluated by the phase screen model all coincide with those theoretical results at high spatial frequency. Simultaneously, the low frequency components could be effectively compensated by the modified subharmonic method. With the increase of the subharmonic order and sample level, the performance evaluated by the phase screen model are closer to the theoretical ones. It has provided an effective way for simulating oceanic turbulence channel for the underwater optical communications.

Since a relatively low attenuation in blue-green range (450-550 nm) was found for optical communications in an underwater environment 1 , underwater wireless optical communication (UWOC) has attracted much attentions 2 . UWOC lends an advantage of broader bandwidths and a potential for higher levels of security compared to its acoustic counterpart 3 . Besides, UWOC can include multiplexing schemes that further increase the capacity, where the known multiplexing scheme is space division multiplexing (SDM), involving spatially modifying different data-encoded light sources for co-aligned propagation and subsequent receiver discrimination and retrieval.
One form of SDM in free space optical communications that has recently received much attention is orbital angular momentum (OAM) [4][5][6][7][8] . A light beam with a helical wavefront carries an OAM value corresponding to   per photon, where  is the reduced Planck's constant and  is an unbounded integer that represents the number of 2π phase changes in the azimuthal direction. Due to different OAM modes are inherently orthogonal 9 , the OAM-multiplexed optical communications can handle a huge amount of information and tremendously increase the capacity of communication link [4][5][6][7][8] .
Very recently, OAM multiplexing technique has been introduced to UWOC 10-14 . For example, J. Baghdady et al. reported 3-Gbit/s UWOC employing 2 OAM modes multiplexing 10 . Y. Ren et al. further increased the UWOC transmission capacity to 4 Gbit/s (directly modulated LD) and 40 Gbit/s (external modulation & frequency doubling) by multiplexing 4 green OAM modes 11 . Also 2-meter underwater transmission of 4-fold green light (520 nm) OAM modes was demonstrated for a UWOC multicasting link in 12 .
At the same time, many theoretic research works on the evolution of optical beam carrying OAM in the underwater environment have been reported [15][16][17][18][19] . For instance, Xu et al. 15 deduced the expression of the cross-spectral density matrix for the propagation of random electromagnetic vortex beams in oceanic turbulence, while Cheng et al. 16 used Rytov approximation theory to calculate the influence of ocean turbulence on the transmission of Laguerre-Gaussian (LG) beam.
www.nature.com/scientificreports www.nature.com/scientificreports/ For all the experimental results [10][11][12][13][14] , the underwater conditions are emulated by using a fixed length tank filled with tap water, and only limited underwater conditions, such as, scattering/turbidity, current, and thermal gradients could be setup; For the theoretical results [15][16][17][18][19] , the analytical solutions are always derived by the complex mathematical calculation on beam expression with Rytov approximation method. The derivations are generally complicated and sometimes there are no analytical solutions. It is desirable to have an easy way to simulate the optical propagating in underwater channel for exploring the performance of a UWOC system or a UWOC system carrying OAM.
In the paper, we present an effective way for simulating oceanic turbulence channel on the propagation beam carrying orbital angular momentum (OAM) using power spectrum inversion method. The influence caused by underwater channel on the intensity and phase of the propagation beam is equivalent to the propagation results when the beam passes through several individual phase screens at regular intervals. At the same time, a modified subharmonic compensation method is used to solve the low frequency loss problem in a large-scale phase fluctuations. The oceanic fluctuation spectrum in the index of refraction proposed by Nikishov et al. 20 is adopted in the model. The theoretic phase structure function and the theoretic propagation characteristics of an OAM beam in oceanic turbulence are utilized to verify the proposed method.

Random Phase Screens Model for Simulating Oceanic Turbulence Channel
In the section, we present the random phase screens model to simulate the beam propagating through the underwater environment. Several individual phase screens generated by power spectrum inversion method at regular intervals is utilized to simulate the influence of ocean turbulence on the propagation beam.
The random phase screens model. Oceanic turbulence environment is a complex system of physical, chemical and biological combinations. The effect of these factors on beam propagation can ultimately be attributed to kinetic energy and temperature gradients, which will caused different refraction index in different regions. Due to the fluctuations in density, salinity and temperature of the oceanic turbulence environment, the variation of the refraction index along the propagation path would lead to the beam intensity fluctuations in a random manner. These intensity fluctuations would then severely influenced beam propagation properties, like intensity distribution, beam wander and scintillation index, and finally degrade the performance of an UWOC system. Hence, an useful optical beam propagation model for oceanic turbulence is necessary for predicting the performance of the UWOC system.
Here, we use successive phase screens at regular intervals to represent the oceanic turbulence, named random phase screens model [21][22][23] , which is numerical simulation of the light field where the continuous random media is represented as a series of random phase screens transverse to the propagation direction. In the model, multiple individual phase screens are introduced to the beam with a regular interval. Each phase screen is a thin sheet that adjusted the phase of the beam. The perturbation effects, such as those induced by oceanic turbulence, could be incorporated by a split-step method which treating propagation and phase perturbations separately in discrete steps along the propagation axis. In each step, the phase perturbations caused by the refractive index fluctuation arising from turbulence are implemented by multiplying the input field with a phase exponential function (a phase screen), and then the beam undergoes a free-space propagation between the adjacent two phase screens. Additionally, the phase screens are statistically independent of the spatial statistics that matched the index of refraction fluctuations. Figure 1 is the schematic diagram of the oceanic turbulence phase screens model, the phase screen is in the x − y plane, and the direction of beam propagation is in z-axis.
Consider an OAM beam in the plane transmitter (z = 0), the initial beam in the spatial domain (the beam field) is A 0 (x, y, z = 0), where x, y are spatial domain coordinates at z = 0. Before the first phase screen, the beam field should be ′ ′ − A x y z ( , , ) 1 according to the Fresnel diffraction formula. www.nature.com/scientificreports www.nature.com/scientificreports/ is the imaginary unit, x′, y′ are the spatial domain coordinates at z = Δz (Δz is the transmission distance between the transmitter to the first phase screen), λ is the wavelength, k = 2π/λ is the wave number (wave vector). After the beam passes through the first phase screen ϕ 1 (x, y), the random phase would impress on the beam's field, and the beam field turns to The phase screen ϕ 1 (x, y) in the simulation may be expressed as a N x × N y array of random complex numbers, whose statistics match those fluctuations of oceanic turbulence.
This procedure is then repeated for the second phase screen by replacing A 0 (x, y, z) with in Eq. (1), and replacing . The process should be repeated until the last phase screen is reached. The above phase screens ϕ 1 (x, y), ϕ 2 (x, y), …, ϕ M (x, y) can be produced by filtering a white Gaussian noise with a phase screen spectrum, followed by Fourier transform. Generally, the relation between the phase screen spectrum F ϕ (k x , k y ) and the spectrum of refractive-index fluctuations spectrum of the refractive index variations Φ(k x , k y ) can be described as x y x y 2 where k = 2π/λ is the wave vector, Δz is the spacing distance between the subsequent phase screens.
Therefore, a phase screen in simulation x y with size L x = N x × Δx, L y = N y × Δy (N x and N y are the array sizes of phase screen, Δx and Δy are the spatial-sampling intervals at x, y directions), is then described as where i is also the imaginary unit, a(m, n), b(m, n) are zero mean Gaussian uncorrelated random numbers, which are described as x y x y 2 2 Here, 〈·〉 denotes ensemble average, Δk x , Δk y are the spectrum-sampling intervals, Δk x = 2π/L x , Δk y = 2π/L y , respectively. a(0, 0) = b(0, 0) = 0, since the zero frequency component does not change the spatial statistics of the fields.
The modified subharmonic compensation method. The split-step method has the advantages of much faster speed and smaller memory requirements, and is suitable for generating large-scale phase screen. However, it suffers from a drawback of the under-sampling for the low spatial frequency components. Near the origin in Fourier space, the spectrum is changing so rapidly that the sampling near the origin is not enough to simulate all low frequency components, especially those frequencies with periods greater than the size of phase screen, such as, the spectral region (−Δk x /2, Δk x /2) and (−Δk y /2, Δk y /2), which would result in an inadequate simulation of the ocean turbulence random phase screen in the large-scale phase fluctuations.
By re-sampling the spectrum near the origin (the central point of phase screen), the subharmonic compensation method can balance the loss of low frequency components inherent in power spectrum inversion method. Here, a modified subharmonic compensation method is presented. Figure 2 shows the modified subharmonic compensation method, denoted by the subharmonic order p and the sample level q, and the spectral-sampling interval would now turn to Δk x(y) = 2π/(q p L x(y) ). Figure 2(a) denotes the subharmonic compensation method for 3 × 3 sample level. Figure 2(b) shows the first order subharmonic compensation in this sample level. The mid area in Fig. 2(a) has now been divided into nine sub areas, and there is 8 new interpolated spectral samples there. Figure 2(c) shows the second order subharmonic compensation in this sample level, the mid area in Fig. 2(b) has further been divided into nine smaller areas.
For the 3 × 3 sample level, there are 8 new interpolated values for the first subharmonic order compensation, and the second subharmonic order compensation is operated on the first compensation result, therefore, there exists 16 new interpolated low-frequency samples in the second subharmonic order compensation. The www.nature.com/scientificreports www.nature.com/scientificreports/ interpolated low-frequency samples are increased when the subharmonic order increase. The more the compensation order is, the more samples at the central point of the power spectrum has, and the more fruitful the low frequency components the phase screen has.
On the other hand, the mid area is divided into 25 smaller areas for the first subharmonic order compensation at 5 × 5 sample level, and there exists 24 new interpolated samples for the first subharmonic order compensation. Continually, the mid area is then divided into 49 smaller areas for the first subharmonic order at 7 × 7 sample level, and there exists 48 new interpolated samples for the first subharmonic order. When the sample level increase, the new interpolated subharmonic samples values grow up quickly.
With the modified subharmonic compensation, the phase screen with low frequency components can be subsequently expressed as 24  Here q × q is the sample level of the modified subharmonic compensation. When q = 3 the Eq. 7 falls back to the techniques proposed by Lane et al. 25 based on the addition of subharmonics with random complex amplitudes to improve the large-scale statistics. a(s, t, p), b(s, t, p) are also zero mean Gaussian uncorrelated random numbers, they are described as Here, the spectrum of the phase screen is In the simulation, the compensated sample values are first subjected to discrete Fourier transform to obtain the added phase distribution generated by the modified subharmonic compensation, and then they are added to the random phase screen created by using power spectrum inversion method. The final phase screen can be represented as the superposition of the initial phase screen and the compensated phase screen, θ = ϕ + θ sh .

Model Verification and Validation
In this section, we verify the validity of the model from two aspects. One is the phase structure function, the other is the propagation property of beam carrying OAM modes through the underwater environment.
Phase structure function. It is shown the mean squared difference for any two points in the index of refraction only depends on the distance between the two points, with no relation to the direction in Kolmogorov model 26 . Hence, the statistics characteristics of the underwater turbulence can be described by a structure function.
According to the results in 27 , the theoretical phase structure function of a plane wave was described as  www.nature.com/scientificreports www.nature.com/scientificreports/ where r 0 was a plane-wave spatial coherence radius, defined as the separation distance at which the modulus of the complex degree of coherence fell to 1/e. Here, the plane-wave spatial coherence radius is  where ε, χ T , ω are defined as Eq. 6. Equation 10 indicates that under Rytov approximation, the Kolmogorov five-thirds power law of wave structure function is also valid for oceanic turbulence in the inertial range if the power spectrum of oceanic turbulence proposed by Nikishov is adopted.
On the other hand, we also can evaluate the phase structure function with the beam field simulated under the random phase screens model. Therefore, the phase structure function can be evaluated as 2 where 〈·〉 is an ensemble average, ϕ(r) is the field amplitude at r point, r, r′ are the two points in the phase screen.
To compare the evaluated phase structure function with the theoretical one, we used 1000 random phase screens to calculate the assemble average in the numerical simulations. Simultaneously, the oceanic turbulence simulation parameters were set as following. ε = 0.01 m 2 /s 3 , χ T = 10 −10 K 2 /s, ω = −4, η = 0.0001, z = 50 m. The wavelength of the beam was 488 nm. The size of the phase screen was 0.8 m × 0.8 m and the sampling points were 81 × 81. Figure 3 shows the phase structure function evaluated under the random phase screens model, the theoretical phase structure function with Eq. 10, and those results with the modified subharmonic compensation method. The oceanic turbulence simulation parameters were set as ε = 0.01 m 2 /s 3 , χ T = 10 −10 K 2 /s, ω = −4, η = 0.0001, z = 50 m. The modified subharmonic compensation method had 3 × 3 sample level, and the subharmonic order was the first order, the second order and the third order. The results showed that the phase structure function of phase screens was very close to the theoretical one (r ≫ η), when the distance between two points was smaller, say, less than 10% of the size of phase screen; With the distance increases, the difference between the theory one and the evaluated one increased. When the distance between two points on the phase screen was close to the size of the phase screen, the evaluated structure function (red curve) dropped very quickly, was far from the theoretical value (black curve).
On the other hand, the evaluated phase structure function from the improved phase screens model had better performance than that one without compensation. As the more and more subharmonic components were added to phase screens model, the evaluated phase structure function approached the theoretical one, which indicated that the subharmonic compensation method could effectively overcome the lack of low frequency components problem in power spectrum inversion method.
We further discuss the effect of the sample level on the compensation. In order to analyze the compensation effect, we define a parameter CE to qualify the proximity between two curves.  www.nature.com/scientificreports www.nature.com/scientificreports/ Figure 4 shows the compensation effect with different subharmonic order (p) and sample level (q), where Fig. 4(a) illustrates the phase structure function with different sample level, Fig. 4(b) shows the compensation effect against the sample level. Simultaneously, the oceanic turbulence simulation parameters were set as same as that in Fig. 3, ε = 0.01 m 2 /s 3 , χ T = 10 −10 K 2 /s, ω = −4, η = 0.0001, z = 50 m. The results showed that the evaluated phase structure function was approaching to the theoretical one when the sample level (q) was increased. For the lower subharmonic order, the compensation effect CE increased quickly with the sample level. However, the compensation effect became constants when the subharmonic order was greater than 5. The higher sample level was, the better the compensation effect was. Figure 5 further discusses the influence of turbulence parameters in phase screens model on the compensation effect. Figure 5(a-c) describe the compensation effect against the rate of dissipation of kinetic energy ε, the rate of dissipation of mean-squared temperature χ T , and the ratio of temperature and salinity contributions to the refractive index spectrum ω, respectively. The results showed that the compensation effect varied little with these parameters (ε, χ T and ω), which indicated that the random phase screens model had a good stability and was available for simulating the oceanic turbulence. where m is the number of radial nodes in the intensity distribution, l is topological charge, giving an OAM of   per photon, πω λ = z / When a LG beam carrying topological charge  0 passes through the oceanic turbulence, not only the intensity of beam, but also its phase distribution will change. OAM is no longer pure and its energy will spread to the adjacent OAM states. The crosstalk would occur between different OAM modes, and the resultant beam can be considered as a superposition of all OAM modes. Since OAM modes with different topological charges are mutually orthogonal, at the receiver, the intensity of  ( ) 0 OAM state can be obtained by using an inverse spiral phase mask ϕ −  i exp( ) 0 since the resultant beam has turned back to Gaussian mode. Moreover, the intensity for the detection probability of each OAM state should be normalized. Hence, the probability of mode  ( ) 0 in the resultant beam in theory can be detected as  www.nature.com/scientificreports www.nature.com/scientificreports/