1.

Introduction

The propagation of optical waves in random media, including atmospheric turbulence, is governed by the stochastic Helmholtz equation. Applying weak perturbation theory, it is possible to achieve closed-form expressions for most field statistics as outlined by Tatarski.¹ Similarly, asymptotic theory^2,3 provides an excellent match to experiment when fluctuations are very large. In between these two regions, extended Rytov theory⁴ (ERT) provides a heuristic analytical model that has found a good match to some experimental data. Attempts at providing solutions in this region that do not rely on approximations involve complex multi-dimensional integrals and are difficult to generalize.^5–8

While not analytical, wave optics simulations (WOS) are known to be an accurate approximation to the stochastic Helmholtz equation. ⁹ This approximation is achieved by discretizing the propagation volume and collapsing the three-dimensional (3D) refractive index fluctuations in each volume segment into a thin phase screen. The optical field is then propagated between adjacent segments via a Fresnel propagation operator. The relative fluctuations at each screen are weak such that the approximation is equivalent to the small perturbation approximation used by Tatarski. The main strength of WOS is that it can provide insight into wave propagation statistics and the performance of systems that attempt to compensate for the effects of the media, and are relatively straightforward to execute.

Though the technique was previously described by others, relative to optical wave propagation, Martin and Flatté¹⁰ were the first to explore the use of WOS for optical waves. Their studies¹⁰ of plane and spherical waves¹¹ in random media described by power laws examines the intensity scintillation as a function of fluctuation strength in the media. Since that time, WOS has become the default tool used in modeling the performance of adaptive optics, beam projection, and free-space optical (FSO) communication systems.^12–16 Simulations involving phase screens are also commonly used to create synthetic image data of scenes blurred by turbulence.^17–19 Until relatively recently, it has been common to model atmospheric turbulence in WOS as a power-law medium described by a 3D energy spectral density defined only in terms of a spatial wavenumber with a $-11/3$ exponent. We would describe this as a purely Kolmogorov model in that it lacks an inner and outer scale.

In 1994, Dalaudier et al.²⁰ published evidence of turbulence in the upper atmosphere from balloon-borne measurements described as non-Kolmogorov in nature. Soon after, Beland²¹ examined its implications to existing turbulence theory. Stribling et al.²² then extended those results to weak fluctuation theory. In this context, non-Kolmogorov turbulence refers to deviation of the power law away from the 2/3 power-law structure function description of turbulence in the inertial subrange between the energy input (outer scale) and dissipation regions (inner scale). The equivalent 3D energy spectral density of turbulent index of refraction fluctuations has an $-11/3$ exponent as a function of spatial frequency. This latter quantity is more often found in modern analytical models that begin with the refractive index fluctuation energy spectrum starting with Toselli.^23–26 These works define non-Kolmogorov turbulence as a deviation in the power-law exponent, $\alpha$, away from the Kolmogorov value of $11/3=3.66$ over the range $3<\alpha <4$. In this way, a shift toward 3 is said to be smaller and to 4, larger. Another common thread of all the works cited below and the work that followed is the prominence of analytical methods as opposed to numerical methods such as WOS in exploring phenomenology.

While Martin and Flattè did use WOS to examine the effect of random media on intensity scintillation, it lacks a connection to modern atmospheric propagation work. It is also useful to understand the limitation of WOS and the effect of a change of power law on those limits. In earlier work,²⁷ we made an early attempt at making these connections and understanding the relationship between the power-law exponent of a random media and scintillation as a function of turbulence strength. Our conclusions were that an explicit inner and outer scale were necessary to make such a connection.

In this work, we describe the results of a comprehensive WOS campaign, started in Ref. 28 for plane waves propagating in uniform non-Kolmogorov turbulence volumes. Results are presented here for intensity scintillation as a function of both Rytov number and power-law exponent. Results are limited to a single propagation geometry featuring a finite inner and outer scale. Recent results have shown that failure to account for these features can result in inaccuracies when evaluating some quantities.^29,30 We find that peak scintillation or normalized intensity variance (NIV) increases with power law. Similarly, the Rytov number where peak NIV occurs also increases. This confirms earlier results²⁸ indicating a shift of the NIV curve up and to the right as power law increases. Across the campaign, WOS is run to their practical limit for this propagation geometry assuming an upper bound of $8192\times 8192$ phase screen samples. Thus, this work also describes the upper bound of WOS as a function of non-Kolmogorov power law using Martin and Flatté’s sampling constraints. Interestingly, here smaller power exponents are limited to Rytov numbers of 7 compared to 12 for Kolmogorov turbulence and larger power-law exponents.

The remainder of this paper is outlined as follows; following this introduction, we provide an overview of the WOS model and the sampling constraints used to define our simulation campaign. In Sec. 3, we describe the campaign in detail and describe our scoring quantities. Results are provided in Sec. 4 followed by conclusions and directions for future work in Sec. 5.

2.

Background

A goal of this work is to understand the practical limits of WOS via phase screens generated by filtering white Gaussian noise by the power spectral density (PSD) spectrum of the turbulence fluctuations and using the split-step propagation method. Excellent descriptions of the split-step propagation technique can be found elsewhere,^12,13 and are not included here.

In approaching this work, we use the sampling constraints outlined by Martin and Flatté¹⁰ for plane waves to ensure that the constraints account for non-Kolmogorov power-law media. This limit is based on the spatial bandwidth of the power-spectrum as sampled by the phase screen. As pointed out indirectly by Martin and Flatté and elsewhere, the power law affects the prominence of high versus low-frequency fluctuations in the media. For this reason, power-law random media with smaller power-law exponents are likely to have more high-frequency fluctuations and therefore require a higher sampling rate compared to the Kolmogorov default.

Before describing those constraints, let us first define the 3D PSD of the index of refraction fluctuations for a turbulence volume described by a power law with an inner and outer scale^31–33

It is common for WOS to use the Fried parameter, $r_0$, to set turbulence strength for phase screens. In previous work,³⁴ one of us compared two methods of setting the phase screen turbulence strength for studies of non-Kolmogorov turbulence. One method used the non-Kolmogorov equivalent Fried parameter so that each screen has the same effective spatial coherence properties. The other technique relied on first normalizing the spectral energy and then scaling the energy in the screen by the equivalent fluctuation energy that would be found in a Kolmogorov turbulence volume. A finding of this work³⁴ was that the latter method provides a better match to ERT when evaluating the scintillation index in terms of the plane-wave Rytov number ${\sigma }_R=(1.23C_n^2{k}^{7/6}{L}^{11/6}{)}^{1/2}$ where $k$ is the optical wavenumber and $L$ is the propagation distance.

In this work, phase screens are generated using the technique described originally in Ref. 33 and extended to WavePy.³⁵ The current version of WavePy uses the sub-harmonic method described by Johansson and Gavel.¹⁵ In previous work,³³ one of us noted a practical limit of around $\alpha =3.8$ for accurately modeling non-Kolmogorov phase screens with an infinite outer scale. Anecdotally, we find that this new method improves phase screen accuracy in terms of a structure-function match and allows us to explore out to $\alpha =3.9$. Though, as we will discuss later as $\alpha \to 4$, the turbulent disturbance becomes equivalent to a pure tilt.³⁶ It follows that improved subharmonic modeling would allow for larger power laws to be modeled more accurately.

In determining sampling rates for our WOS campaign, we use the spatial bandwidth definition described by Martin and Flatté to determine the required spatial bandwidth, $R_{\chi }$, in each phase screen as

Another requirement described by Martin and Flatté¹⁰ and others¹² limits the log-amplitude variance to not exceed 0.1% or 10% between propagation steps. This is effectively a forward-scattering requirement or the Rytov approximation and can be set in terms of $\alpha$ as

3.

Methods

The main result of this work is an extensive WOS campaign that evaluates the scintillation index as a function of the plane-wave Rytov number in non-Kolmogorov turbulence. Here, we intend to hold the path length fixed and vary only the turbulence strength of the volume in terms of the plane-wave Rytov number, ${\sigma }_R$, with the goal of performing simulations out to the limit allowed by the sampling constraints described in Sec. 2. For the purposes of this work, the maximum number of samples per screen will be $2^{13}=8192$ square. While $2^{14}=16384$ is feasible, for the scenario considered here that sampling is required only for conditions well within the saturation region and where wall-clock executions times are considerable.

The WOS campaign was run in two batches, the first was completed out to the maximum Rytov number allowed for one of three power laws, $\alpha$. Power laws values of 3.1 and 3.9 cover the lower and upper limits of practical simulation, while the Kolmogorov power value of 11/3 (3.66) was included as a baseline. A second, more limited, campaign was conducted for 25 values of $\alpha$ between 3 and 4. Values started at $\alpha =3.02$ and increased increments of 0.02 up to a value of $\alpha =3.2$. From there, $\alpha$ was incremented in steps of 0.1 up to 3.8 over the regime including the Kolmogorov value ($\alpha =3.66$) included as a special case. From $\alpha =3.8$ to 3.98 steps were reduced again to 0.02. The objective here being to identify values of power-law exponent where the simulation is no longer valid, or fidelity is reduced.

As mentioned in the previous section, we fixed the side length $D=1\mathrm{m}$ and varied the sampling rate $N$ as required. For the other WOS parameters, the wavelength is set to $\lambda =1\mu \mathrm{m}$, and the path length to $L=5\mathrm{km}$. For Kolmogorov turbulence, these parameters allow for a full exploration of the various propagation regimes described by theory without resorting to values of $C_n^2$ unlikely to be measured in the field, in the range from ${10}^{-17}$ to ${10}^{-12}{\mathrm{m}}^{-2/3}$.

Using Eqs. (3) and (5) in Sec. 2, we are able to evaluate the maximum Rytov number that can be evaluated for a given sampling rate, $N$. Simultaneously, for a given value of $\alpha$, we can find the minimum sampling rate needed to evaluate a specific Rytov number to ensure accurate results. As outlined in Ref. 27, it is possible to model atmospheres with power laws very near to the lower-bound of $\alpha =3$ via WOS. However, as the value of $\alpha$ approaches 3, the number of samples required to simulate even weak turbulence is very large. This both unnecessarily limits the range of volume turbulence strengths we can explore and increases wall-clock simulation time. Using a value of $\alpha =3.1$ as our lower-limit case provides some relief in these regards while allowing a detailed examination of the behavior at smaller power-law exponents.

In our previous related works,^29,30 we found that it is likely that WOS simulations without finite inner, $l_0$, and outer, $L_0$, scales are likely not valid. For this work, the outer scale size was set to $L_0=1\mathrm{m}$ and the inner scale to $l_0=0.005\mathrm{m}$. The former matches the fixed screen size, $D$, while the latter is slightly larger than the largest spatial sampling rate of $D/N$ for $N=256$ of $\mathrm{\Delta }x=0.0039\mathrm{m}$.

3.1.

Limits on WOS in Non-Kolmogorov Turbulence

In Fig. 1, we evaluate Eqs. (3) and (5) to visualize the limits on WOS in non-Kolmogorov turbulence as described in this paper. In Fig. 1(a), the maximum Rytov number allowable for each screen size of $N=2^n$ where $n=\mathrm{1,2},\dots ,14$ is presented. Of note here is that the maximum Rytov number that can be simulated with a screen size of $N=2^{14}$ is close to ${\sigma }_R=24$. Also, that for all values of sampling rate, the power law allowing the highest turbulence strength is close to $\alpha =3.8$ it is also at this point where additional sampling most increases the range of Rytov numbers covered. While as $\alpha \to 3$, the same benefit is not conveyed and the maximum turbulence strength is ${\sigma }_R=6$ for $N=2^{13}$ and ${\sigma }_R=7$ for $N=2^{14}$. Though not captured in our campaign, we also note that as $\alpha \to 4$ the max Rytov number drops to zero. Again, this is because media with a power law of exactly 4 are pure tilts and cannot be described via a structure function. This also explains the difficulty observed in Ref. 18 of generating accurate phase screen statistics at larger values of $\alpha$ and further underscores the practical limit of $\alpha \approx 3.9$ for WOS modeling.

Fig. 1

Limits on the maximum plane-wave Rytov number, ${\sigma }_R$, that can be simulated for a specific power-law exponent, (a) $\alpha$ or (b) screen size, $N$. (b) The three values of $\alpha$ indicated are those chosen for detailed evaluation.

4.

Results

Figure 2 shows the results of the first of our two WOS campaigns. Each specified value of $\alpha$ scintillation as measured in our WOS is plotted out to the maximum Rytov number allowed according to Eq. (2) for a screen with $N=8192$ samples. As described in Table 1, $\alpha =3.1$ is limited to a maximum of 7 while $\alpha =3.66$ and 3.9 go to ${\sigma }_R=12$. The results here are similar to those presented elsewhere and previously,^23,27,28 but in this instance include, an explicit inner and outer scale in the turbulence spectrum sized based on the simulation geometry.

Fig. 2

Plane-wave simulation of scintillation index, ${\sigma }_I^2$, as a function of Rytov number, ${\sigma }_R$, for $\alpha =3.1,3.66,3.9$. For each power law, simulations were conducted out to the maximum Rytov number specified for a sampling rate of $N=8192$. In this scenario, the path length was $5000\hspace{0.17em}\mathrm{m}$ and the wavelength was $\mu \mathrm{m}$. Error bars indicate the variance about the mean at each data point. For all values of $\alpha$, the outer-scale size was $L_0=1\mathrm{m}$ and inner-scale was $l_0=5\mathrm{mm}$.

Examining Fig. 2, we see that the curves shift up and to the right relative to power-law exponent. Earlier onset of saturation for the $\alpha =3.1$ case is also observed though without the sharp peak in the focusing region^9,37 observed in Ref. 28. We attribute this difference to the explicit, finite, inner-scale in the turbulence spectrum in Eq. (1). This spectral feature filters, or limits, the degree of small-scale fluctuations that drive scintillation and loss of spatial coherence resulting in saturation. It is worthy note that even if an inner scale is not an explicit feature of the turbulence power spectrum one is included, implicitly, by the simulation sampling rate. In WOS campaigns, like ours, where the sampling rate increases with turbulence strength, the inner scale would decrease with each change in sampling rate. Analytic theory indicates³⁷ that this change should hasten onset of saturation and reduce peak scintillation. Here, our chosen value of $l_0=5\mathrm{mm}$ is larger than the largest sampling rate, $\mathrm{\Delta }x$, used in the WOS to avoid this complication.

In some previous works,^23,28 it was noted that scintillation may increase without bound as $\alpha \to 4$. We observe that this behavior is missing in Fig. 2. Here, also, the inclusion of an explicit finite outer scale likely plays a role by limiting the relative energy in large-scale fluctuations.²⁵ However, the observation holds that, all things being equal, peak scintillation increases with power law and occurs at higher equivalent turbulence strength. So that, relative to Kolmogorov turbulence, an increase the value of the power-law exponent increases peak scintillation and results in relatively higher scintillation in the saturation region. Conversely, if the value of the power-law exponent of the medium is decreased peak scintillation in the focusing and saturation regions are smaller.

On the other hand, for plane waves propagating in weak-to-moderate strength turbulence volumes (${\sigma }_R<1$) the situation is reversed. In this regime, the relationship between turbulence strength and intensity scintillation is driven by the power law of the medium. Consequently, we observe that mediums with smaller power laws experience higher fluctuations in intensity. Likewise, mediums with smaller power laws experience less focusing and enter saturation sooner as a function of turbulence strength.

The results of the second batch of our WOS campaign are shown in Fig. 3 and confirm these results further. All parameters are the same as those in Fig. 2, but the value of ${\sigma }_R$ is limited to a maximum of 4 while $\alpha$ was varied as described in Sec. 2. Figure 3(a) can be compared directly to Fig. 2 and confirms that the behavior observed there is approximately continuous and increases monotonically as a function of $\alpha$.

Fig. 3

(a) Result of the second WOS campaign for 25 values of $\alpha$ in the range $3<\alpha <4$ between $0<{\sigma }_R<4$. (b) Surface plot of (a).

In Fig. 4, the peak value of scintillation index [Fig. 4(b)] and the plane-wave Rytov number where the peak scintillations were observed [Fig. 4(a)] are plotted as power law, $\alpha$, is varied. In both figures, the least-squares linear fit to the data is also plotted. The relationship between $\alpha$ and maximum scintillation is observed to have a slope of 0.9 while the slope of the Rytov number where peak scintillation occurs and $\alpha$ is 1.8. In this latter case, values necessarily are restricted to the specific values of Rytov number simulated during the campaign resulting in binning or quantization errors. In both cases, as we noted earlier, WOS accuracy is not certain as we exceed the bounds of $3.1<\alpha <3.9$ also contributing some uncertainty to these results. Regardless, the evidence of a positive linear relationship between power law and both peak intensity scintillation and the volume turbulence strength at which that peak scintillation occurs is clear. For this propagation scenario, it is also clear that the turbulence value of peak scintillation moves to the right at about double the rate of the peak scintillation increases.

Fig. 4

(a) Value of ${\sigma }_R$ where maximum scintillation occurs and (b) peak value of ${\sigma }_I^2$ as a function of power-law index $\alpha$.

To our knowledge, this quasi-linear, positive, monotonic relationship between power law of the medium and peak scintillation in terms of Rytov number has not been reported elsewhere. Toselli et al.²³ used ERT to generate analytical models of plane-wave scintillation as a function of power law and Rytov number. However, this model does not include finite inner and outer scales. Though, that work is consistent with these findings and our previous work.^27,28,34 In subsequent works by Toselli et al.³⁸ and others^{24–26,39,40} on non-Kolmogorov turbulence findings are often reported for a fixed value of ${\tilde{C}}_n^2$ in the weak to moderate regime while varying propagation distance instead of normalized volume turbulence strength as is done here. As pointed out by Charnotskii⁴¹ because the units of ${\tilde{C}}_n^2$, vary with $\alpha$ it is not possible to use a common length scale for comparing results in any meaningful way. As far as we can ascertain, all previous works report a relationship similar to Ref. 38 where peak scintillation occurs at a $\alpha <11/3$ in the region of $3.2<\alpha <3.3$ and is lower on either side of the peak going to zero as $\alpha \to 3$ and to a small value as $\alpha \to 4$.

To place our work in this context, in Fig. 5, we have plotted the normalized scintillation index, ${\sigma }_I^2/\mathrm{Max}({\sigma }_{I,\alpha }{)}^2$, as a function of power-law index, $\alpha$. As pointed out earlier, we see that in the weak turbulence regime, for a fixed volume turbulence strength, scintillation decreases with $\alpha$. In contrast, scintillation increases if the volume is characterized as strong or deep (${\sigma }_R\gg 1$). However, in this figure, when ${\sigma }_R=1.5$ scintillation peaks near the center of the range at $\alpha =3.5$ and decreases if $\alpha$ is increased or decreased.

Fig. 5

Scintillation index for three fixed values of integrated turbulence strength in terms of ${\sigma }_R$ as a function of $\alpha$. Each trace is normalized to the maximum value of scintillation index over all values of $\alpha$.

Thus, relative to those other works the relationship between $\alpha$ and ${\sigma }_I^2$ varies with integrated turbulence strength in the volume. Also, because of the ambiguity introduced via a fixed $C_n^2$, this behavior may change depending on the power law so that any of the behaviors observed in Fig. 5 may apply. For example, if $\lambda =1.55\mu \mathrm{m}$, $L=1000\mathrm{m}$, and ${\tilde{C}}_n^2=7\times {10}^{-14}{\mathrm{m}}^{3-\alpha }$, as in Ref. 38, the Rytov number varies from 3.8 to 1.18 to 0.73 for $\alpha =3.1,3.66$, and 3.9. Consider, then, that for this one example all three of the behaviors in Fig. 5 may be observed. This observation further emphasizes that some measure of equivalent volume turbulence strength must be used when trying to understand the impact of turbulence power-law exponent on beam propagation and imaging.

5.

Conclusions and Future Work

In this work, we explored the limits of WOS plane waves propagating in non-Kolmogorov turbulence and with power-law exponents in the range $3<\alpha <4$. At the upper bound, the medium becomes a pure tilt, and therefore, WOS is limited by the size of the screen or the fidelity of the low-spatial frequency compensation mechanism; sub-harmonics for example. Conversely, as the power law approaches the lower-bound the medium become spatially uncorrelated, and the spatial bandwidth required to properly simulate the medium becomes very large. Accordingly, the number of samples required by the simulation is also large. Both WOS requirements can be ameliorated by the inclusion of a finite inner and outer scale. Indeed, as we have shown elsewhere,³⁰ if these values are not defined explicitly in the power-spectrum model they are implicit in the simulation model. Finally, if we take the practical limit of WOS to limited to $N=16384$ samples, the maximum plane-wave Rytov number that can be simulated in any power-law medium is $\sigma =24$ when $\alpha =3.8$.

The results of our analysis on the limits of WOS inform two simulation campaigns exploring the interplay between Rytov number, $\alpha$, and scintillation index. The first campaign aimed to evaluate scintillation as a function of Rytov number for values of power law $\alpha =3.1,3.66,3.9$. For each value of $\alpha$ simulations were undertaken out the maximum ${\sigma }_R$ allowed for a screen size of $N=8192$. Consistent with other results, the maximum scintillation observed and the Rytov number where the peak occurs increase with power law. At small power-law exponents, scintillation quickly saturates as the volume turbulence strength increases. In this work, we did not observe a strong focusing peak seen in our previous works as $\alpha \to 3$. We attribute this to the use of a finite inner scale larger than the sampling rate. Similarly, in previous works scintillation appeared to increase without bound at power laws near $\alpha =4$. However, if we include an explicit outer scale on the order of the screen size peak scintillation increases but eventually rolls off into saturation. These findings are consistent with theory that attributes small-scale fluctuations to scintillation strength and large-scale fluctuations to the peak of intensity fluctuations in the focusing region.

Our second simulation campaign aimed to empirically evaluate the effect of medium power-law index on scintillation. We found here that, in contrast to some previous works, intensity scintillation increases monotonically with the power law. Likewise, the ${\sigma }_R$ where peak scintillation occurs moves to the right at twice the rate peak scintillation increases. We assert that the differences between these results and previous works are due to the proper scaling of turbulence strength with the power law used here. The aforementioned inner and outer scale may also be contributing factors.

Over the course of this work, the contradictory nature of our findings and previous works led us to revisit our results to ensure our WOS parameters were correct. As a result of these explorations,²⁷ we were led to conclude that, at least for the simulation scenario explored here, these parameters are a mostly a second-order effect. That is to say, the overall trends remain as long even if the sampling rate or number of screens changes. This finding has been recently confirmed by Wijerathna.⁴² In one of his last works, Flatté and Gerber⁹ indicated that WOS is an exact solution to the stochastic Helmholtz equation over all bounds when properly configured. Much has been made of the qualification, but based on this work and others, we feel WOS will usually provide accurate results even if they are not precise.²⁷

There is more work to be done in this area. For example, it may be worthwhile to further explore Fresnel zone effects by varying wavelength and path length. This work looked only at plane-wave propagation. Therefore, it may be interesting to see if further accommodations¹² are needed to account for beam expansion and beam-wander for diverging sources. This work also explicitly does not include a Hill bump or other features in the dissipation range. This exclusion is purposeful as it is not at all clear the nature of this feature when the medium is non-Kolmogorov. Finally, a straight-forward extension would be to explore the effect of larger values of $l_0$ and smaller values of $L_0$ for different power laws as $\alpha$ is varied.