Computer simulations of optical turbulence in the weak- and strong-scattering regime: angle-of-arrival fluctuations obtained from ray optics and wave optics

Abstract. It is known that certain geometrical-optics predictions often agree well with optical turbulence field observations even though theoretical constraints for ignoring diffraction may be violated. Geometrical optics assumptions can simplify analyses, and ray optics can significantly reduce simulation computation time. Here, an investigation into angle-of-arrival fluctuations is presented involving wave optics and geometrical (ray) optics computer simulations of a plane wave of visible light propagating through a turbulent refractive-index field. The simulation and Rytov-based theory results for the variances of aperture-filtered angle-of-arrival fluctuations generally agree well for weak scattering (Rytov variance, σR2≲0.2), but for increasing Rytov variance, the simulation results demonstrate a positive slope that can be significantly shallower than that predicted by the theory. For weak-to-moderate scattering regimes (σR2≲2.67), a comparison of the ray and wave results show they match for aperture diameters greater than about two Fresnel lengths. This result is consistent with a previous theoretical analysis by Cheon and Muschinski. For the strongest scattering case studied (σR2=26.7), the wave and ray simulations match for aperture diameters greater than about 10 Fresnel lengths. For smaller apertures, we attribute the disparity between the wave and ray simulation results to a Fresnel filtering effect.


Introduction
Understanding the physics of optical wave front distortions caused by turbulence in the atmospheric refractive-index field is essential for the design and operation of optical systems used in free-space optical communication, surveillance, navigation, remote sensing, astronomy, and directed-energy technologies. [1][2][3][4][5][6] If particle scattering, absorption, coupling, and depolarization effects are negligible then there are two phenomena that determine the optical field: refraction and diffraction. The combined effects of refraction and diffraction lead to phase-gradient fluctuations (angle-ofarrival fluctuations) and intensity fluctuations (scintillation), both of which limit the performance of a variety of optical systems.
The analysis of propagation scenarios through both deterministic and stochastic refractive-index fields may be substantially simplified if diffraction effects can be neglected, and considerable theoretical effort has been devoted to establishing criteria for the negligibility of diffraction effects. [1][2][3]7,8 These theoretical constraints, however, tend to be excessively conservative, as pointed out by Strohbehn 50 years ago. 9 In other words, certain geometrical-optics predictions often agree well with field observations, but it remains unclear, from the theoretician's perspective, why that is so.
In addition to simplifying theoretical analyses, geometrical optics and ray tracing can also substantially reduce computation time for numerical simulations of propagation through turbulence. The utility of geometrical and ray optics is exemplified in a variety of publications from the last several decades. For example, Churnside and Lataitis 10 presented a simple geometrical optics equation for the variance of beam displacements caused by propagation through weak turbulence, and geometrical optics simulation models have been developed by Yuksel et al. 11 to investigate effects of aperture averaging on intensity fluctuations. More recently, Lachinova et al. 12 have compared the performance of a brightness function-based numerical simulation method with a wave optics-based method for incoherent imaging through turbulence. The brightness function at the receiving pupil is estimated using ray-like trajectories that are perturbed by the turbulent refractive index gradients. These authors found that the brightness approach can reduce computation times by several orders of magnitude. Other recent examples include ray-tracing methods applied for simulating laser beam-steering effects in turbulent media, 13 the derivation of turbulence strength profiles from radiometer temperature measurements with a geometrical ray-tracing analysis, 14 and the application of ray tracing through phase screens for the simulation of images from optical survey telescopes. 15 Since the 1980s, powerful phase-screen techniques have been developed [16][17][18][19] that allow optical propagation through atmospheric turbulence to be accurately simulated well into the "strong-scattering" regime, that is, for scenarios where the negligibility of diffraction effects and the "weak-scattering" assumption become invalid, or at least questionable. Such simulations can be used to test and refine the traditional, often excessively conservative, criteria for the negligibility of diffraction effects.
Here, we present and discuss computer simulations of plane waves of visible light (λ ¼ 500 nm) propagating through fully turbulent, homogeneous, and isotropic refractive-index fields. We study the optical angle-of-arrival (AOA) fluctuations by means of both a wave-optics propagation simulation and a ray-tracing simulation through a turbulent path that is characterized by a sequence of phase screens with Kolmogorov statistics. The path length assumed is L ¼ 2 km and the refractive-index structure parameter C 2 n ranges from 10 −16 m −2∕3 to 10 −13 m −2∕3 , which corresponds to Rytov variances σ 2 R that span from 0.0267 (weak scattering) to 26.7 (strong scattering), where σ 2 R ¼ 1.23C 2 n k 7∕6 L 11∕6 and k ¼ 2π∕λ is the optical wave number. For the two simulation methods, we examine the variance and histograms of the aperture-filtered (aperture-averaged) AOA fluctuations for aperture diameters ranging from 0.6 to 13 Fresnel lengths. We compare the simulation results with theoretical predictions 20 based on the Rytov theory, which accounts for both refraction and diffraction effects but does not capture strong-scattering effects. Our specific interests are: (1) to better understand under what conditions the Rytov theory for AOA fluctuations becomes invalid, (2) to reveal the behavior of the AOA fluctuations in strong scattering conditions, and (3) to investigate the consequences of ignoring diffraction and the applicability of ray tracing for modeling AOA fluctuations.

Approach
A plane wave is propagated through a set of split-step phase screens that model atmospheric turbulence [ Fig. 1(a)] and similarly, a bundle of parallel rays is traced through the same screens [ Fig. 1(b)]. The wave optics simulation provides a coherent wave result that accounts for refraction, diffraction, and interference. On the other hand, the ray-tracing simulations provide the trajectories of a large number of individual rays, where the phases and amplitudes are left undefined.
The simulation parameters of interest are defined in Table 1 and the values used for the simulations are listed. The turbulence along the path is assumed to be statistically homogeneous and isotropic. The path is divided into equal length segments and a phase screen is positioned at the center of each segment. The screens are created with the Fourier filtering method and have Kolmogorov statistics. To model an infinite outer scale, a random tilt component is added to each screen that compensates for turbulent wave-front tilts with length scales greater than the grid width (see Ref. 18 and references therein). Zero inner scale is assumed for Fig. 1 Illustration of the split-step phase screens that model turbulence in (a) wave and (b) ray simulations. L is the propagation path length and D is the diameter of the aperture that is applied in the observation plane for the aperture-averaged AOA. the turbulence spectrum; however, the period of the highest spatial frequency of a screen is twice the grid spacing, so the modeled inner scale is no smaller than 2 mm.
The number of turbulence screens N required along the path is an integer and should be greater than ð10σ 2 R Þ 6∕11 , which ensures each screen corresponds to weak scattering, and therefore phase-only screens are adequate. 19 The minimum value of N for the strongest scattering case is 22 screens. Fewer screens can be used for the weaker scattering cases, which reduces computational time. However, we implemented more than the required minimum number of screens for the weak scattering cases to provide a better representation of the continuous volume turbulence that is assumed along the path.
The grid parameters are chosen to provide a grid spacing of Δx ¼ λL∕S, where S is the grid side length. This sampling choice has been shown to make artifacts that arise from the discrete Fourier transform implementation of the Fresnel diffraction integral to be sufficiently small to be negligible. 21 The resulting grid spacing Δx of 1 mm allows about 3× sampling of the minimum Fried parameter (r 0 ) value encountered in our study. The propagation length L, wavelength λ, and four choices of turbulence structure parameter C 2 n (10 −16 , 10 −15 , 10 −14 , and 10 −13 m −2∕3 ) correspond to Rytov variances σ 2 R of 0.0267, 0.267, 2.67, and 26.7, which range in scattering strength from weak (σ 2 values for the four cases are calculated to be 21.3, 5.3, 1.33, and 0.33 cm.
Ultimately, our interest is to examine the variance and histogram of the fluctuations of the aperture-averaged AOA. This requires running the simulation many times for independent turbulence realizations, obtaining the apertureaveraged AOA for each realization, and then computing the variance.

Wave Optics Simulation
The initial plane wave at the entrance plane, z ¼ 0, is a unit-amplitude, zero-phase field that is propagated in steps between the entrance plane and the observation plane using a discretized form of the Fresnel diffraction integral. 21 At the i'th intermediate plane, the phase screen ψ i is applied in the usual way as u 0 i ¼ u i expðjψ i Þ, where u i is the incident field and u 0 i is the field that exits the plane. The discrete Fourier transform-based simulation produces a periodicity of the electric field in the x-y plane across the grid boundaries. In addition, the compensating random tilt component in the turbulence phase screens creates some "ringing" artifacts in the receiving plane near the grid boundaries. To minimize unwanted periodicity and edge effects, the aperture is only applied to the center 0.4 × 0.4 fraction of the receiving plane grid.
One approach to find the aperture-averaged AOA for a given turbulence realization is to numerically compute the field phase gradients at the observation plane grid points and average the results within the aperture region. However, direct calculation of the phase gradients at the grid points is problematic for strong scattering because of branch points (implying discontinuous phase changes from 2π to 0, or 0 to 2π) that appear in the wave front and correspond to theoretically infinite phase gradients.
For the results presented here, we use an approach where the wave front within the aperture for a given turbulence realization is numerically focused with a focus transmittance function. 21 The xand y-direction centroid positions of the intensity spot at the focal plane are found and the apertureaveraged AOA is the angle defined by the ratio of the centroid displacement (from the optical axis) and the focal length. This AOA is computed for many independent turbulence realizations and the variance is calculated. This approach has the advantage that it is representative of a practical AOA measurement that can be made with a lens and camera system. 22 It also effectively weighs the angle contributions by the spatial field amplitude in the pupil, where, for example, the amplitude is vanishing near a branch point so the sometimes extremely large gradients associated with these features are not significant for the AOA variance calculation.

Ray Tracing Simulation
For the ray simulation, as shown in Fig. 1(b), the source is comprised of a set of rays where each ray is positioned at a grid point corresponding to the wave simulation and is initially directed parallel to the optical axis. For a given turbulence realization, the rays are traced from one plane to the next where at the i'th plane, the phase screen changes the x-direction ray angle as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 1 ; 3 2 6 ; 1 8 3 where θ i and θ 0 i are the incoming and outgoing ray angles and dψ i dx is the phase-screen gradient in the x-direction. Ray angle changes in the y-direction are computed similarly using the y-direction gradients. A three-point calculation is used to compute the screen gradients, and the gradients are interpolated at the position where the ray strikes the plane. The xand y-direction angles for the rays that eventually fall  within the aperture region at the observation plane are averaged to obtain the aperture-averaged AOA. Our approach for tracing rays through the split-step phase screens is similar to the method presented by Lachinova et al., 12 however, we do not compute the brightness function that is part of their incoherent image simulation.

Theory
The variance of the aperture-filtered AOA fluctuations for a circular aperture and for a plane wave propagated through homogeneous and locally isotropic inertial subrange turbulence is given by Cheon and Muschinski 20 as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 2 ; 3 2 6 ; 6 9 4 where D is the aperture diameter, q ¼ D∕ ffiffiffiffiffi ffi λL p is the aperture diameter normalized by the Fresnel length, β ¼ 0.5216 is a constant, and Γ is the gamma function. Equations (2) and (3) were developed assuming the Rytov approximation for the propagating field. This result is applicable for all values of q and assumes infinite outer scale and zero inner scale.

Observation Plane Illustration
Before presenting the AOA results, we first examine the character of the rays and waves arriving at the observation plane for a single turbulence realization. (1) ray-density fields obtained from the ray simulation and (2) intensity and phase fields obtained from the wave optics simulation. The ray density fields were obtained by counting the rays within areas defined by the binning values. The binning is relative to the grid spacing Δx, for example, 10 × 10 binning indicates an area of 10 mm × 10 mm. The bin size was selected to emphasize the most apparent characteristic feature sizes within the frame. A contour plot algorithm was also applied to the ray density plot to reduce the blockiness of the display.
Common spatial features are apparent in the intensity and ray-density fields as shown in Fig. 2. The corresponding features are related to the effects of refraction as ray tracing does not account for the interference and diffractive effects provided in the wave simulation. The phase cross-sections, displayed modulo-2π, illustrate the increasing wave front complexity as a function of the Rytov variance. As noted previously, after the ray and wave fields are generated at the observation plane as shown in Fig. 2, a circular aperture of diameter D is imposed and the methods described in Sec. 2 are applied to find the aperture-averaged AOA field for the particular turbulence realization.

Variances of Aperture-Filtered Angle-of-Arrival Fluctuations
The variances of aperture-averaged AOA fluctuations are computed from many independent turbulence realizations. Figure 3 shows results for the theoretical [Eq. (2)] and simulated aperture-averaged AOA variances as a function of the normalized aperture diameter q for the four values of C n 2 listed in Table 1. Each data point was estimated from 1000 turbulence realizations. The general behavior is that the AOA variance decreases as a function of increasing aperture diameter. Both simulation results generally agree well with the AOA theory 20 for weak scattering (e.g., σ 2 R ≲ 0.2) but decrease significantly below the theory curve as the Rytov variance increases. The results also show the agreement between theory and simulation is better in stronger scattering for larger apertures (e.g., σ 2 R ¼ 2.67 and q ≳ 10). The over-prediction of the AOA variances by the Rytov theory is likely a result of the well-known failure of the Rytov approximation in the strong-moderate-to strong-scattering regime. We note that the Rytov theory is also known to over-predict intensity scintillation as a function of the Rytov variance in the moderate-and strong-scattering regimes.
In general, the ray-based and wave-based simulation variances compare favorably in Fig. 3. The largest discrepancy is a 15% higher variance for the rays in strong scattering with a small aperture (σ 2 R ¼ 26.7 and q ¼ 0.63). For the weak-tomoderate scattering cases [ Figs. 3(a)-3(c)], the ray and wave results match well for q ≳ 2. That is, diffraction effects are negligible in the weak-and moderate-scattering regime if the aperture diameter is larger than about twice the Fresnel length. This result is consistent with the theoretical analysis by Cheon and Muschinski 20 who found that for weak scattering, the AOA variance predicted by geometrical optics deviates from that predicted by the Rytov theory by <1% if q > 1.65. For smaller q, the wave simulations predict a lower AOA variance than the ray simulations. For the strongest scattering case [ Fig. 3(d)], the wave and ray variance values converge for q ≳ 10. Figure 4 shows an alternative presentation of the results where aperture-averaged AOA variances for the small aperture (q ¼ 0.63; D ¼ 0.02 m) and the large aperture (q ¼ 12.6; D ¼ 0.4 m) are plotted as a function of the Rytov variance, σ 2 R . The Rytov theory curves demonstrate the linear dependence described by Eq. (2). For the small aperture [ Fig. 4(a)], the simulations show a positive slope as a function of scattering strength that is significantly shallower than the theory line. For large aperture diameters, the Rytov-based theory is expected to collapse to the geometrical optics solution. 20 This behavior is suggested in the large aperture case [ Fig. 4(b)], where the simulation results follow the theory line much more closely. Figure 4 also shows a small difference between the ray-based and wave-based results for the small aperture (q ¼ 0.63), but essentially no difference for the large aperture (q ¼ 12.6). As noted previously, geometrical optics is known to be questionable for q < 1, which is consistent with the fact that our wave-optics simulations produce smaller AOA variances than the AOA variances produced by our ray-tracing simulations. We attribute this to "Fresnel filtering." The Fresnel-filter kernel in the wave-number representation of the AOA variance for plane waves predicted by the Rytov theory has the form K F ðxÞ ¼ 1 þ ð2π∕x 2 Þ sinðx 2 ∕2πÞ, 20 where x ¼ κðλLÞ 1∕2 is the wave number normalized by the Fresnel wave number 1∕ðλLÞ 1∕2 . The case x ≪ 1 (Fresnel length is small compared with κ −1 , the wave-front tilt length scale under consideration) is the geometrical-optics limit where the Fresnel length has no effect on the AOA variance. For x comparable to and >1, the Fresnel length (and therefore diffraction) does play a role: K F ðxÞ becomes smaller than K F ðxÞ in the geometrical-optics limit. Therefore, diffraction causes the AOA variance predicted by wave optics to be smaller than the AOA variance predicted by geometrical optics.
That is, a finite Fresnel length acts as a characteristic filter scale similar to the filtering effects due to a finite inner turbulence length scale or a finite aperture diameter. Essentially, diffraction causes wave-front tilt components with length scales comparable to and smaller than the Fresnel length scales to be smoothed out. The Fresnel-filtering effect is significant when the aperture size and the inner scale of turbulence are both comparable to or smaller than the Fresnel length, and it is negligible otherwise. The Fresnel filtering is accounted for in the wave-optics simulations but not in the ray-tracing simulations.
To examine the AOA distribution results in detail, we generated histogram estimates of the probability density   Gaussian distributions. The two different slopes in the q-q curves in Fig. 5(d) correspond to a difference in the AOA variance for the wave and ray results, which was discussed previously for small q. Histogram and q-q plots for strong scattering are shown in Fig. 6. These distributions are also Gaussian and q-q plot slope differences in Figs. 6(d)-6(f) are consistent with the variance results in Fig. 3(d). We expect a Gaussian distribution for the aperture-averaged AOA if the propagation through turbulence and aperture averaging represents a large sum of independent tilt components such that the central limit theorem applies.

Conclusions
The simulations illustrate that after propagation through homogeneous turbulence, wave and ray results share obvious intensity/ray-density and AOA characteristics related to the effects of refraction. We find that both the wave and ray simulation results generally predict variances of the apertureaveraged AOA fluctuations that agree well with the Rytov theory for weak scattering (e.g., σ 2 R ≲ 0.2), but as scattering strength (i.e., the Rytov variance) increases, the simulation results demonstrate a positive slope that can be significantly shallower than that predicted by the weak-scattering (Rytov) theory. This over-prediction by the Rytov theory is likely a result of the failure of the Rytov approximation in the strongscattering regime.
Some insights related to the negligibility of diffraction effects are gained from a comparison of the wave and ray simulation results. For weak-to-moderate scattering regimes (σ 2 R ≲ 2.67), the ray and wave results match well for apertures where q ≳ 2. This result is consistent with the theoretical analysis by Cheon and Muschinski 20 who found that for weak scattering, the AOA variance predicted by geometrical optics deviates from that predicted by the Rytov theory by <1% if q > 1.65. For the strongest scattering case studied here (σ 2 R ¼ 26.7), the wave and ray variance values agree well for q ≳ 10. For smaller apertures, we find the wave simulation variances are smaller than those produced by the ray simulations, and we attribute this to a "Fresnel filtering" effect associated with diffraction. Both the wave and ray AOA distributions are found to be Gaussian. David Voelz is a professor in the Klipsch School of Electrical and Computer Engineering at New Mexico State University. His research interests include spectral and polarization sensing, laser beam propagation through atmospheric turbulence, laser communications, imaging theory, and astronomical instrumentation development. He earned his PhD EE degree from the University of Illinois in 1987. He is a fellow of SPIE.
Erandi Wijerathna is a PhD student in the Klipsch School of Electrical and Computer Engineering at New Mexico State University. Her research interests include laser beam propagation through atmospheric turbulence, imaging theory, and polarimetric lidar. She earned her BS (honors) in physics from the University of Colombo, Sri Lanka, in 2011. She received her MS degree in physics and in electrical engineering from New Mexico State University in 2016. in meteorology, U Hannover, Germany, 1998) is an atmospheric physicist. He was a CIRES research scientist at the University of Colorado at Boulder, Colorado (1998Colorado ( -2004) and a professor of electrical engineering at UMass Amherst (2004-2011). Since 2011, he has been a senior research scientist at NorthWest Research Associates, Boulder, Colorado. During the last 28 years, he has conducted research on atmospheric turbulence and atmospheric wave propagation.
Xifeng Xiao received her BS and MS degrees in physics from Xiamen University Fujian, China, in 1998 and 2001, respectively. She earned her MS and PhD degrees in electrical engineering from New Mexico State University in 2004 and 2008, respectively. Currently, she is a software engineer at Ball Aerospace, Albuquerque, NM. Her research interests include simulation and modeling of free-space laser communication, liquid-crystal polarization, and demonstration and implementation of AOTF technologies.