Modelling and simulating complex media is a great challenge due to turbulent fluctuations in the temperature, density, salinity (water), and humidity (air) of the environment that occur in both time and space. Therefore, it is far easier to model turbulence by treating it stochastically, as a random process that is both homogenous and isotropic, following Kolmogorov turbulence theory [24]. This problem can then be further simplified by taking a continuous environment, and representing it in the form of discrete phase screens, which distort an optical wave only at that point in time and space. Between the phase screens, the optical wave propagates in a vacuum. This is the basis for the split-step optical wave simulation, which models the complex medium in the form of thin phase screens placed along the path of propagation, shown below.
Though the split-step method is useful for simulating accurate optical wave propagation with relative ease, it is not without its drawbacks. Primarily, the use of Kolmogorov theory of turbulence limits the strength of turbulence of each screen to weak fluctuations of the propagating beam [24]. In order to quantify the strength of the turbulence using the scintillation in the propagating beam, we make use of the Rytov method, assuming a line-of-sight system. For the Rytov method, the optical field can be written as shown in Eq. 4.1, where represents the propagation of the beam in a vacuum, and is the complex phase perturbation. The complex phase perturbation can be written as shown in Eq. 4.2 to calculate successive perturbations, such as the perturbation that is the result of each discrete screen. The complex phase perturbation can also be written in terms of the log-amplitude perturbation , and phase perturbation , given in Eq. 4.3 [28].
The Rytov method allows us to compute the total log-amplitude variance , or Rytov variance, of a propagating spherical wave using Eq. 4.4. The Rytov variance is a function of the wave number k, the total propagation distance , the propagation distance z, the number of screens n, and the refractive-index structure constant . The turbulence level must result in a total Rytov variance of less than 0.25, with each screen having a Rytov variance of less than 0.1, providing the first simulation constraint. The strength of turbulence in simulation is controlled by the rate of dissipation of kinetic energy per unit mass of fluid , the rate of dissipation of mean-squared temperature , which are used to estimate via Eq. 4.5, where the constant [29].
When designing the simulation, there are a number of other parameters that need to be controlled and changed depending on the environmental conditions, as well as other constraints that must be met to produce an accurate simulation. The primary parameters are the simulation resolution N, the side length of the source and observation planes D, grid spacing , the propagation distance between screens , the wavelength of the light λ, the beam waist radius , and the number of phase screens along the path of propagation n. Additionally, a value for c must be selected, which is a constant that represents the model’s sensitivity to turbulence. Typically, , where a value of 2 captures ~97% of the light and a value of 8 captures ~99% [28]. Selection of these parameters is informed by the three inequalities given in Eq. 4.6, 4.7. and 4.8.
To use these constraints, it is also necessary to know the spatial coherence diameter , commonly known as the Fried parameter, for the spherical wave propagating through the complex environment. The spatial coherence radius of a spherical wave propagating in oceanic turbulence, given in Eq. 4.9, is a function of the wave number k, propagation distance , rate of dissipation of kinetic energy , rate of dissipation of mean-squared temperature , and ratio of temperature and salinity contribution [30]. The spatial coherence radius of a spherical wave propagating in oceanic turbulence is related to the spatial coherence diameter by Eq. 4.10 [31].
The final constraint to meet is the minimum number of screens needed, related to the maximum distance between screens. This is given in Eq. 5.11, and is a function of the grid spacing , simulation resolution N, and wavelength λ. As a result of Eq. 5.32, the minimum number of screens is given in Eq. 5.12 [28].
Once all of these constraints are met, the simulation is accurately parameterized. However, to produce accurately estimated results, the thin phase screens placed along the path of propagation must be properly generated.
Thin phase screens are a discrete realization of a randomly-modeled turbulence. These means to generate a random thin phase screen, computer-generated random numbers must be transformed into a two-dimensional screen that has the same statistical properties as the modeled turbulence. This is most commonly done using a method based on the Fourier Transform (FT), where the numbers are drawn from a power spectrum representing the desired medium [28] [32]. For this work, the Nikishov spectrum is used because it models the oceanic environment, meaning that it accounts for changes to both the temperature and salinity of the water [33]. Using the method outlined by Korotkova in [29], the Nikishov power spectrum, given in Fig. 4.2, was individually plotted for just the temperature (Eq. 4.13) and just the salinity (Eq. 4.14). The full Nikishov power spectrum is given in Eq. 4.15. In these equations C and A represent constants, is the rate of dissipation of kinetic energy per unit mass of fluid, represents the rate of dissipation of mean-squared temperature and salinity, is the ratio of the temperature and salinity contribution (ranging -5 to 0, temperature to salinity), and η is the Kolmogorov inner-scale number.
When the spectra given by Eq. 4.13 and Eq. 4.14 are normalized and plotted, the figure given in Fig. 4.5 is produced.
Figure 4.5: Normalized temperature and salinity Nikishov spectra [29]
From this power spectrum, a phase screen such as the one given in Fig. 4.6 is generated. This can be done using the provided function “ft_nikishovPS.m” given in Appendix A, Sect. 10.
Figure 4.6: FT Nikishov Random Phase Screen
However, this screen alone does not produce results that match what is experimentally observed. This is because when the FT method is used to generate the phase screen, only the high spatial frequencies are captured, as the grid cannot be sampled with a high enough resolution to capture the low spatial frequencies. As can be seen in Fig. 4.5, there is significant power in the low spatial frequencies of the temperature-only Nikishov spectrum. To compensate, the subharmonic (SH) method is commonly used, combining random draws of Zernike polynomials with non-uniform sampling in the spatial-frequency domain. This is proposed in [34] and utilized by Schmidt [28]. The SH method produces screens such as the one given in Fig. 4.7, which stands in stark contrast to Fig. 4.6. A SH phase screen can be generated using “sh_nikishovPS.m” provided in Appendix A, Sect. 11.
Figure 4.7: SH Nikishov Phase Screen
In order to create an accurate phase screen, the FT and SH phase screens can simply be added together, producing the phase screen show in Fig. 4.8. This can be seen in “ft_sh_nikishovPS.m”, provided in Appendix A, Sect. 12.
Figure: 4.8: FT and SH Nikishov Phase Screen
When experimenting with phase screen generation, it was observed that the FT and SH phase screens have two different, but each well-documented, effect on the propagating beam. The FT screen, with power at the high spatial frequencies, causes the scintillation of the beam’s intensity, characterized by a shimmering of the irradiance of the light. Intuitively, it is clear how small-scale changes to the beam’s phase (Fig. 4.6) would create scintillation, as the shimmering is small-scale as well. The SH screen, with power at the low spatial frequencies, causes beam wander, pushing the beam around the receiver and causing the vortex of the LG beams to break, if strong enough. This too can be intuitively seen by the large-scale changes in phase. Both effects are typical to experimentally observe, and when the FT and SH screens are combined, both effects are observed in simulation. As an additional check for accuracy, 20 produced phase screens were correlated to check for proper randomization, and a low correlation of 0.411 was observed, shown in Fig. 4.9. The correlation analysis was completed in MATLAB using “Screen_Vali.m” provided in Appendix A, Sect. 13.
Figure 4.9: Correlation of 20 Nikishov Random Phase Screens
Before conducting a simulated evaluation of our beams, we studied the effect of changing simulation parameters to determine the optimal conditions and simulation parameters before the lengthy simulations were run. Depending on the resolution of the screens, and the number of realizations generated, the simulation can take 24-48 hours of continuous running to produce a full-sized, usable dataset. To examine the simulation, screens corresponding to three levels of turbulence were first generated, using . This is summarized in Table 4.2, with the screens in Fig. 4.10 corresponding to strong turbulence, Fig. 4.11 to moderate turbulence, and Fig. 4.12 to weak turbulence.
Table 4.2: Summary of Simulated Turbulence Levels and Parameters
Figure 4.10: Strong turbulence phase screens (A) FT phase screen, (B) SH phase screen, and (C) combined phase screen
Figure 4.11: Moderate turbulence phase screens (A) FT phase screen, (B) SH phase screen, and (C) combined phase screen=
Figure 4.12: Weak turbulence phase screens (A) FT phase screen, (B) SH phase screen, and (C) combined phase screen
As is clear from Figures 4.10, 4.11, and 4.12, there is a significant difference in the generated FT phase screens, with very little difference in the SH phase screens. This describes how the high-spatial frequencies dominate the distortion of the light in strong turbulence while the low-spatial frequencies dominate in weak turbulence.
A single LG beam, with order 0 and topological charge 3, was then propagated at varying distances in all three levels of turbulence, using N = 2048, D = 10 cm, n_scrns = 15 and lambda = 633 nm. The first propagation, through strong turbulence, is shown below in Fig. 4.13.
Figure 4.13: Beam in strong optical turbulence at distances Z = 5, 10, and 30 m
As the simulation is dominated by the scintillation, there is little utility in propagating a full alphabet through turbulence at this strength, as the beam quickly will become unrecognizable, to a human or the CNN.
The same LG beam was then propagated through moderate turbulence at varying distances, given in Fig. 4.14. This has more useful applications, because as can qualitatively be seen, the distortions are significant, but not overpowered by scintillation.
Figure 4.14: Beam in moderate optical turbulence at distances Z = 5, 10, and 30 m
Qualitatively, these results look very similar to what is experimentally observed in the lab, though in reality the propagation distance is much shorter. This is a byproduct of the weak fluctuation constraint of the split-step simulation.
Figure 4.15: Beam in weak optical turbulence at distances Z = 10 and 50 m
Figure 4.16: Beam in weak optical turbulence at distances Z = 100, 150, and 200 m
In figures 4.15 and 4.16, the same LG beam is propagated through weak optical turbulence over varying distances. For this case, it is particularly interesting to examine the effect of the turbulence on the phase of the beam, and the impact this has on the vortex of the beam. In the case of Z = 100 m, the beam starts to become noticiably stretched, but the integrity of the vortex remains intact. At Z = 150 m, the singularity in the phase that corresponds with the vortex of the beam begins to break, corresponding to an ovaling of the beam and thinning of the ring that characterizes the intensity pattern. At Z = 200 m, the singularity is fully broken, which corresponds to the broken vortex in the intensity distribution. It has been shown that the CNN relies on vortex structure to accurately classify the beam [25]. Though propagating over unrealistic distances, the simulation still provides us with valuable insight into what causes the vortex of a beam to break. In order to compensate for distortions, it is important to better understand their root cause.
Code and tutorials for the Split-Step Random Phase Screen Simulation can be found in the MATLAB section - linked here!
References
[24] L. C. Andrews and R. L. Phillips, Laser Beam Propagation through Random Media, 2nd Edition, SPIE, 2015.
[25] S. Avramov-Zamurovic, C. Nelson and J. M. Esposito, "Experimentally evaluating beam scintillation and vortex structure as a function of topological charge in underwater optical turbulence," Submitted to Optics Communication, 2021.
[26] Meadowlark Optics, "User Manual 1920x1152 XY Phase Series Spatial Light Modulator With PCIe Controller," Frederick, CO.
[27] Fastec Imaging, "TS3, TS4, and TS5 High Speed Cameras Operator's Manual," Fastec Imaging Corporation, San Diego, CA, 2015.
[28] J. D. Schmidt, Numerical Simulation of Optical Wave Propagation, Bellingham, WA: SPIE Press, 2010.
[29] O. Korotkova, "Light propagation in a turbulent ocean," Progress in Optics, vol. 64, pp. 1-43, 2019.