Randomized Spectral Sampling for Efficient Simulation of Laser Propagation through Optical Turbulence

TL;DR

This paper introduces a novel spectral sampling method for atmospheric turbulence phase screens that improves simulation efficiency and accuracy, especially for long-range optical propagation, by leveraging Fourier transform properties.

Contribution

The paper presents a new spectral sampling technique based on Fourier shift properties that enhances turbulence simulation accuracy without extra computational costs.

Findings

Method produces low spatial frequency distortions effectively.

Performance matches or exceeds existing methods for finite outer scale turbulence.

Advantages over standard methods for long-range atmospheric optical simulations.

Abstract

We present a new method for the generation of atmospheric turbulence phase screens based on the frequency shift property of the Fourier transform. This method produces low spatial frequency distortions without additional computation time penalties associated with methods using subharmonic subgrids. It is demonstrated that for simulations of atmospheric turbulence with finite outer scales, the performance of our method with respect to the statistical phase structure function of the screen meets or exceeds other methods with respect to agreement with theory. We outline small-scale accuracy issues associated with modelling non-Kolmogorov spectral power laws using existing techniques, and propose a solution. For simulations of long-range propagation through atmospheric optical turbulence, our method provides various advantages over standard methods.

Tables

Table 1. Table 1: Hybrid method percent RMS error compared to theory over region of interest computed using 5,000 phase screens’ x − limit-from 𝑥 x- and y − limit-from 𝑦 y- axes for 2048 × 2048 2048 2048 2048\times 2048 grid size, parameterized by outer scale, L 0 subscript 𝐿 0 L_{0} , as well as the number of subharmonic constellations, N p subscript 𝑁 𝑝 N_{p} . The L 0 subscript 𝐿 0 L_{0} ’s given in the first column have units of meters, and the simulated spatial domain was a one meter by one meter area.

-	Number of Subharmonic Constellations, $N_p$
${𝑳}_{\mathrm{𝟎}}$	1	2	3	4	5	6	7	8	9	10
${\mathrm{𝟏𝟎}}^0$	1.7	1.7	1.7	1.7	1.7	1.7	1.7	1.7	1.7	1.7
${\mathrm{𝟏𝟎}}^{.5}$	2.8	3.0	3.0	3.0	3.0	3.0	3.0	3.0	3.0	3.0
${\mathrm{𝟏𝟎}}^1$	2.4	2.7	2.7	2.6	2.6	2.6	2.6	2.6	2.6	2.6
${\mathrm{𝟏𝟎}}^{1.5}$	2.8	3.4	2.9	2.7	2.7	2.7	2.7	2.7	2.7	2.7
${\mathrm{𝟏𝟎}}^2$	2.5	1.1	1.4	1.2	1.3	1.3	1.3	1.3	1.3	1.3
${\mathrm{𝟏𝟎}}^{2.5}$	2.6	1.1	1.3	1.1	1.3	1.3	1.3	1.3	1.3	1.3
${\mathrm{𝟏𝟎}}^3$	7.5	5.3	1.4	1.6	1.5	1.5	1.5	1.5	1.5	1.5
${\mathrm{𝟏𝟎}}^{3.5}$	10	1.6	1.5	1.5	1.1	1.3	1.4	1.3	1.3	1.3
${\mathrm{𝟏𝟎}}^4$	14	1.8	1.8	0.8	0.9	0.9	1.0	1.2	1.2	1.1
${\mathrm{𝟏𝟎}}^{4.5}$	7.4	7.0	1.8	1.6	1.9	1.9	1.5	1.5	1.6	1.6
${\mathrm{𝟏𝟎}}^5$	14	6.4	3.0	1.4	1.4	2.0	2.4	2.5	2.4	2.6

Table 2. Table 2: Traditional method percent RMS error compared to theory over region of interest in similar configuration as described in Table 1 .

-	Number of Subharmonic Constellations, $N_p$
${𝑳}_{\mathrm{𝟎}}$	1	2	3	4	5	6	7	8	9	10
${\mathrm{𝟏𝟎}}^0$	2.1	2.2	2.2	2.2	2.2	2.2	2.2	2.2	2.2	2.2
${\mathrm{𝟏𝟎}}^{0.5}$	5.9	4.8	4.8	4.8	4.8	4.8	4.8	4.8	4.8	4.8
${\mathrm{𝟏𝟎}}^1$	5.4	4.5	4.9	4.8	4.8	4.8	4.8	4.8	4.8	4.8
${\mathrm{𝟏𝟎}}^{1.5}$	14	7.6	13	14	14	14	14	14	14	14
${\mathrm{𝟏𝟎}}^2$	24	3.3	11	15	16	16	16	16	16	16
${\mathrm{𝟏𝟎}}^{2.5}$	30	7.9	8.5	17	20	20	20	20	20	20
${\mathrm{𝟏𝟎}}^3$	31	9.7	6.1	16	22	24	24	24	24	24
${\mathrm{𝟏𝟎}}^{3.5}$	35	15	2.6	8.8	17	22	23	23	23	23
${\mathrm{𝟏𝟎}}^4$	36	15	2.2	8.2	15	20	23	24	24	24
${\mathrm{𝟏𝟎}}^{4.5}$	36	16	2.7	7.6	14	19	22	24	25	25
${\mathrm{𝟏𝟎}}^5$	37	18	5.0	5.3	12	16	19	22	23	24

Table 3. Table 3: Minimum percent RMS error ( ℰ ℰ \mathcal{E} ) and associated number of subharmonic constellations ( N p subscript 𝑁 𝑝 N_{p} ) for Randomized Hybrid Method (middle two rows) versus traditional method (bottom two rows) for N p ≤ 10 subscript 𝑁 𝑝 10 N_{p}\leq 10 . Data was collected over 50,000 phase screens for 1024 × 1024 1024 1024 1024\times 1024 grid, μ x = 1 subscript 𝜇 𝑥 1 \mu_{x}=1 and μ y = 2 subscript 𝜇 𝑦 2 \mu_{y}=2 .

$𝜶$	3.1	3.2	3.3	3.4	3.5	3.6	$\frac{\mathrm{𝟏𝟏}}{\mathrm{𝟑}}$	3.7	3.8	3.9
${𝑵}_{𝒑}$	2	3	4	3	7	10	9	5	7	10
$𝓔$	.65	.33	.82	.86	.26	.46	.61	.55	1.9	16
${𝑵}_{𝒑}$	1	2	2	2	2	3	3	4	5	9
$𝓔$	8.5	9.7	9.6	9.5	10	10	10	10	9.7	9.6

Table 4. Table 4: Minimum percent RMS error ( ℰ ℰ \mathcal{E} ) and associated number of subharmonic constellations ( N p subscript 𝑁 𝑝 N_{p} ) for Randomized Hybrid Method (middle two rows) versus traditional method (bottom two rows). Data was collected over 5,000 phase screens for 2048 × 2048 2048 2048 2048\times 2048 grid, μ x = 1 subscript 𝜇 𝑥 1 \mu_{x}=1 and μ y = 2 subscript 𝜇 𝑦 2 \mu_{y}=2 .

$𝜶$	3.1	3.2	3.3	3.4	3.5	3.6	$\frac{\mathrm{𝟏𝟏}}{\mathrm{𝟑}}$	3.7	3.8	3.9
${𝑵}_{𝒑}$	4	3	4	4	5	7	12	8	17	17
$𝓔$	1.4	1.3	1.1	1.5	1.3	.95	1.5	1.2	1.3	4.0
${𝑵}_{𝒑}$	1	2	2	2	2	3	3	4	5	10
$𝓔$	11	11	11	8.2	9.5	13	12	11	13	9.1

Equations

$$D_n(\vec{r})=⟨{\left(f\left(\vec{r}+\overrightarrow{r_o}\right)-f\left(\overrightarrow{r_o}\right)\right)}^2⟩$$

$$D_{n}(\vec{r})=2\int_{-\infty}^{\infty}\Phi_{n}(\vec{\kappa})\big{[}1-\cos(\vec{r}\cdot\vec{\kappa})\big{]}\mathrm{d^{3}}\vec{\kappa}$$

$$\theta (j,l)=$$

$$=$$

$$\tilde{c}(\beta ,\gamma )=(a+ib)\cdot k\cdot \sqrt{2\pi \Delta z\Delta {\kappa }_x\Delta {\kappa }_y{\Phi }_n(\beta {\hat{e}}_x+\gamma {\hat{e}}_y)}$$

$${\theta }_R(j,l)=\sum _{n,m=-M/2}^{M/2-1}$$

$$\cdot exp\left[i\left(j\Delta x(n\Delta {\kappa }_x+\delta {\kappa }_x)+l\Delta y(m\Delta {\kappa }_y+\delta {\kappa }_y)\right)\right]$$

$$C(j,l)=M^2\cdot {\mathcal{F}}_2^{-1}[\tilde{C}(n,m)]$$

$${\theta }_R(j,l)=exp[i(j\Delta x\delta {\kappa }_x+l\Delta y\delta {\kappa }_y)]\cdot C(j,l)$$

$${\Phi }_n(\vec{\kappa })=0.033\cdot C_n^2\cdot |\vec{\kappa }{|}^{-\text{nicefrac}113}$$

$$\displaystyle\Phi_{n}(\vec{\kappa})=0.033\cdot C_{n}^{2}\cdot\textit{f}_{n}\bigg{(}\frac{\kappa}{\kappa_{l}}\bigg{)}\cdot\frac{\exp{\big{(}-\kappa^{2}/\kappa_{l}^{2}\big{)}}}{\big{(}\kappa^{2}+\kappa_{0}^{2}\big{)}^{\nicefrac{{11}}{{6}}}}$$

$${\text{f}}_n(x)=1+1.802x-0.254x^{\text{nicefrac}76}$$

$$D_{\theta }(\overrightarrow{r_{\perp }})=4\pi k^2\Delta z{\int }_{-\infty }^{\infty }{\Phi }_n\left(\overrightarrow{{\kappa }_{\perp }}\right)\left[1-\cos \left(\overrightarrow{r_{\perp }}\cdot \overrightarrow{{\kappa }_{\perp }}\right)\right]{\mathrm{d}}^2\overrightarrow{{\kappa }_{\perp }}$$

$$D_{\theta}(\rho)=8\pi^{2}k^{2}\Delta z\int_{0}^{\infty}\kappa_{\rho}\Phi_{n}(\kappa_{\rho})\big{[}1-J_{0}\big{(}\rho\kappa_{\rho})\big{]}\mathrm{d}\kappa_{\rho}$$

$${\mathcal{E}}_x=\sqrt{\frac{2}{M}\sum _{j=1}^{M/2}{\left(\frac{D_x(j\Delta x)-D_{\theta }(j\Delta x)}{D_{\theta }(j\Delta x)}\right)}^2}$$

$${\mathcal{E}}_y=\sqrt{\frac{2}{M}\sum _{l=1}^{M/2}{\left(\frac{D_y(l\Delta y)-D_{\theta }(l\Delta y)}{D_{\theta }(l\Delta y)}\right)}^2}$$

$$\mathcal{E}=100\%\times \frac{{\mathcal{E}}_x+{\mathcal{E}}_y}{2}$$

$$⟨\overline{{\theta }^2}⟩={\pi }^{2}L{\int }_0^{\infty }{\kappa }^3{\Phi }_n(\kappa )\left[1+\frac{2\pi }{(\kappa f{)}^2}sin\left(\frac{(\kappa f{)}^2}{2\pi }\right)\right]A\left(\frac{D\kappa }{2}\right)\mathrm{d}\kappa$$

$${\theta }_R(j,l)=\sum _{n,m=-M/2}^{M/2-1}$$

$$\cdot exp\left[i\left(j\Delta x\left(n\Delta {\kappa }_x+\delta {\kappa }_x\right)+l\Delta y\left(m\Delta {\kappa }_y+\delta {\kappa }_y\right)\right)\right]$$

$$\displaystyle\theta_{out}(j,l,p)=3^{-2p}\sum\limits_{n,m=-1}^{1}\big{(}1-$$

$$\displaystyle\cdot\tilde{c}\bigg{(}\frac{n\Delta\kappa_{x}+\delta\kappa_{x}}{3^{p}},\frac{m\Delta\kappa_{y}+\delta\kappa_{y}}{3^{p}}\bigg{)}$$

$$\displaystyle\cdot exp\Bigg{[}i\bigg{(}j\Delta x\frac{n\Delta\kappa_{x}+\delta\kappa_{x}}{3^{p}}$$

$$\displaystyle\hskip 36.98857pt+l\Delta y\frac{m\Delta\kappa_{y}+\delta\kappa_{y}}{3^{p}}\bigg{)}\Bigg{]}$$

$${\theta }_{in}(j,l)=3^{-2(N_p+1)}$$

$$\displaystyle\cdot exp\Bigg{[}i\bigg{(}j\Delta x\frac{n\Delta\kappa_{x}+\delta\kappa_{x}}{3^{N_{p}+1}}$$

$$\displaystyle\hskip 36.98857pt+l\Delta y\frac{m\Delta\kappa_{y}+\delta\kappa_{y}}{3^{N_{p}+1}}\bigg{)}\Bigg{]}$$

$${\theta }_H(j,l)=$$

$$\displaystyle D_{n}(x,y,z)=\tilde{C}_{n}^{2}\ \bigg{(}\frac{x^{2}}{\mu_{x}^{2}}+\frac{y^{2}}{\mu_{y}^{2}}+z^{2}\bigg{)}^{\frac{\alpha-3}{2}}$$

$$\displaystyle\Phi_{n}(\kappa_{x},\kappa_{y},\kappa_{z})=\frac{A(\alpha)\tilde{C}_{n}^{2}\mu_{x}\mu_{y}}{\big{(}\mu_{x}^{2}\kappa_{x}^{2}+\mu_{y}^{2}\kappa_{y}^{2}+\kappa_{z}^{2}\big{)}^{\nicefrac{{\alpha}}{{2}}}}$$

$$\displaystyle A(\alpha)=\frac{\cos{\big{(}\frac{\pi\alpha}{2}\big{)}}\Gamma(\alpha-1)}{4\pi^{2}}$$

$$\sigma_{x}^{2}=2\pi\Delta zk^{2}\iint_{K_{out}}\Phi_{n}(\kappa_{x},\kappa_{y},0)\big{[}1-\cos(\Delta x\cdot\kappa_{x})\big{]}\mathrm{d}\kappa_{x}\mathrm{d}\kappa_{y}$$

$$\sigma_{y}^{2}=2\pi\Delta zk^{2}\iint_{K_{out}}\Phi_{n}(\kappa_{x},\kappa_{y},0)\big{[}1-\cos(\Delta y\cdot\kappa_{y})\big{]}\mathrm{d}\kappa_{x}\mathrm{d}\kappa_{y}$$

$$\displaystyle K=\Big{\{}(\kappa_{x},\kappa_{y}):$$

Full text

1 Introduction

The split-step propagation method for modeling optical propagation through atmospheric turbulence has been widely used in statistical analysis of beam propagation since its introduction by Fleck, et al. [1, 2]. This method remains popular in simulations of long-range, linear optical propagation due to computation time advantages associated with using the Fast Fourier Transform (FFT) algorithm to compute the Discrete Fourier Transform (DFT), which are used in both the optical propagation and atmospheric distortion algorithms of the cited method. In simulations of nonlinear optical propagation phase screens are also widely used [3, 4, 5, 6], including when supporting studies centering on filamentation [7, 8, 9, 10, 11, 12]. However, due to circular shift symmetry and aliasing affects associated with the DFT, significant effort has gone into the development of computational methods which add subharmonic, low-spatial-frequency components to the atmospheric screens [13, 14]. Additionally, Zernike-polynomial-based methods [15] and other creative methods [16] have been pioneered partly to address this issue, include methods using randomized sampling of the turbulence energy spectrum introduced by Charnotskii [17]. However, due to the computational efficiency of leveraging the FFT algorithm, DFT-based methods for phase screen generation remain popular [18, 19, 20, 21, 22, 23].

We present a modified method which exploits the Fourier transform shift theorem [24], which also extends to the DFT [25], in order to include low frequency components in an FFT-centric method in a straightforward manor. For many applications, this method may provide sufficient phase screen accuracy relative to theory without additional computational penalty associated with subharmonics and other methods. Additionally the method can be combined with the subharmonic method of Lane, et. al [13] in order to give very accurate results across a range of spectral models of practical and theoretical importance. Section 2 discusses the basic method and results of use for bounded spectral models where it is well suited. Section 3 discusses the use of the method in concert with subharmonics and other improvements, and includes a detailed analysis of computation time considerations. Finally, Section 4 summarizes results and discusses possible applications.

2 Randomized FFT-based Sampling

2.1 Algorithm

Central to the study of optical turbulence is the index of refraction structure function, $D_{n}(\vec{r})$, defined as [26]:

[TABLE]

As for any real function, the structure function is related to the field’s three dimensional energy spectrum, $\Phi_{n}(\vec{\kappa})$, by the following relationship [26]:

[TABLE]

Traditional phase screen simulations using square grids [1, 2, 13, 14, 19] approximate the continuous energy spectrum as discrete, and generate complex screens as per:

[TABLE]

where $i=\sqrt{-1}$, $\Delta x$ and $\Delta y$ are the grid spacings in the $x-$ and $y-$directions, $M$ the number of grid points along each axis, and $\Delta\kappa_{x}$ and $\Delta\kappa_{y}$ are the spatial wavenumber grid spacings in the $x-$ and $y-$directions. We note the relation, $\Delta x\Delta\kappa_{x}=\Delta y\Delta\kappa_{y}=2\pi/M$. $\tilde{c}(n\Delta\kappa_{x},m\Delta\kappa_{y})$ is a random function defined by:

[TABLE]

where $\beta$ and $\gamma$ are dummy-variables, $k=2\pi/\lambda$ is the optical wavenumber (with $\lambda$ the wavelength in vacuum), $a$ and $b$ are Gaussian random variables with variances of one, and $\hat{e}_{x}$ and $\hat{e}_{y}$ are unit vectors in the $x-$ and $y-$directions.

Eq. 3 represents a Fourier series using elements which are all harmonic across the spatial domain, creating screens which are periodic [13]. Investigating the effect the $\tilde{c}(0,0)$ term has on the summation in Eq. 3, we note that it results in only the addition of a constant phase, piston term across all of $\theta(j\Delta x,l\Delta y)$. This piston term term does not contribute tip, tilt, focus or defocus effects at any scale, or otherwise contribute to the behavior of the propagating field. Should $\Phi_{n}(0,0)$ have a large enough value quantization error [25] will result. $\tilde{c}(0,0)$ is commonly set to zero in practice [14, 19], which avoids these issues.

In defining a new type of complex phase screen, $\theta_{R}$, we propose a more meaningful use of the point closest to the $\kappa$-space origin by virtue of:

[TABLE]

where $\delta\kappa_{x}$ and $\delta\kappa_{y}$ are random variables described by a uniform distribution bound by $\pm\Delta\kappa_{x}/2$ and $\pm\Delta\kappa_{y}/2$, respectively. This offsets the lowest wavenumber grid point away from the origin, along with also translating the rest of the sampling grid in the frequency domain. By allowing $\tilde{c}(n\Delta\kappa_{x}+\delta\kappa_{x},m\Delta\kappa_{y}+\delta\kappa_{y})$ to define the elements of a matrix, $\tilde{C}(n,m)$, we find that Eq. 5 is implementable via inverse FFT as per:

[TABLE]

[TABLE]

where $\mathcal{F}_{2}^{-1}[\ \cdot\ ]$ in Eq. 6 denotes the two dimensional inverse FFT operator. Alternatively, the FFT can be used directly with appropriate conditioning of the $C$ matrix.

Fig. 1 juxtaposes the sampling methods discussed, where the convention of setting $\tilde{c}(0,0)$ to zero is reflected by the lack of a traditional spectral sampling grid point at the origin. $\theta_{R}(j,l)$ represents a single complex-number valued phase screen, with the real and imaginary parts therein defining a pair of real-number valued phase screens. Simulated atmospheric turbulence distortion is applied via multiplication of our complex propagating beam or wave by $exp\left(i\cdot Re\left[\theta_{R}(j,l)\right]\right)$ or $exp\left(i\cdot Im\left[\theta_{R}(j,l)\right]\right)$, where $Re$ and $Im$ functions represent taking the real and imaginary parts of an array, respectively.

$\theta_{R}(j,l)$ and the real valued phase screens it produces no longer exhibit periodicity, and will have domain-wide low spatial frequency distortions. Additionally, we find the increase in computational delays associated with the use of the algorithm given in Eq. 5 relative to Eq. 3 the context of split-step wave optics simulations to be only 25.4% – 28.5% for the grid sizes displayed in this section. The phase screens shown in Fig. 2 and Fig. 3 were generated using the popular approximation to the Hill Spectrum [27] developed by Andrews [28, 29], known commonly as the modified atmospheric spectrum. This spectral model is discussed in detail in subsequent sections. Note that tip and tilt components can be seen across the $x=0$ and $y=0$ axes, respectively, in the screen shown in Fig. 3a. Additionally, Fig. 2a displays periodicity [13], in that should one circularly shift [25] the phase screen in either or both directions no sharp discontinuities would be apparent within the boundaries of the screen. To help visualize this, we have added Fig. 2b and 3b to illustrate the presence and lack of periodicity resulting from the relevant algorithms.

2.2 Results for Bounded Spectral Models

The most widely used three dimension spectral model of atmospheric turbulence is derived from A. Kolmogorov’s famous $\nicefrac{{2}}{{3}}$’s law as [26, 29]:

[TABLE]

This model is popular due to its simple formulation, and approximate accuracy when the beam statistics of interest are within the inertial subrange of turbulence. However, this spectrum diverges at the $\kappa$-space origin leading to unphysical properties such as containing infinite energy, divergent covariances, lack of a viscosity driven minimum feature size, and lack of a maximum feature size [29]. For these reasons, we will refer to this type of spectral model as unbounded. As we will demonstrate in the next section, additional modifications to the FFT-based algorithm may be required to accurately model unbounded spectral models using phase screens, as is the case for Kolmogorov turbulence.

We turn our attention to a practical atmospheric turbulence spectral model which accounts for inner scale, $l_{0}$, and outer scale, $L_{0}$, bounds on the inertial subrange, as well as intricacies of the experimentally observed energy spectra at higher spatial frequencies [30, 31, 27]. The modified atmospheric spectrum is given by Andrews [28, 29]:

[TABLE]

where $\kappa_{l}=3.3/l_{0}$, $\kappa_{0}=2\pi/L_{0}$, $\kappa=|\vec{\kappa}|$, and we define the function $f_{n}$ as:

[TABLE]

As this spectral model does not present the same complications as that of Eq. 8, we refer to this as a bounded spectral model.

To assess the accuracy of the revised method, we must designate our metrics of interest. We had previously defined the refractive index three dimensional structure function, $D_{n}(\vec{r})$, via the spectral model of interest in Eq. 1. The structure function we are interested in, however, is that of an atmospheric phase screen which approximates the cumulative effects of optical propagation through a finite propagation distance, $\Delta z$. We denote this function as $D_{\theta}(\vec{r_{\bot}})$, where $\vec{r_{\bot}}=x\hat{e_{x}}+y\hat{e_{y}}$. $D_{\theta}(\vec{r_{\bot}})$ is defined by the two dimensional integral over all $\vec{\kappa_{\bot}}=\kappa_{x}\hat{e_{x}}+\kappa_{y}\hat{e_{y}}$ as per [14]:

[TABLE]

For the modified atmospheric spectrum, though we are aware of closed-form approximations of turbulent structure functions for plane waves applicable to our analysis [32], we have instead developed our theoretical structure function via numerical integration of the equivalent form for isotropic turbulence:

[TABLE]

where $\rho=|\vec{r_{\bot}}|$, $\kappa_{\rho}=|\vec{\kappa_{\bot}}|$, and $J_{0}$ denotes the zeroth-order Bessel function of the first kind.

It is well documented that aliasing effects associated with the FFT-based propagation step of the split-step algorithm make parts of the simulation domain unusable [1, 19, 33, 34]. For this reason, a region of interest must be defined, which drives properties of the simulation. Number of grid points, simulated resolution, as well as the propagation distance between screens must be chosen carefully [19]. This requires consideration of many factors, including wavelength, coherence lengths, aperture sizes, etc. As a practical matter, many studies explicitly dedicate half of the $x-$ and $y-$domain of simulation as guard bands to protect against edge effect aliasing [35, 36]. Additionally, the requirement of grid sizes greater than or equal to twice the size of the limiting apertures (or regions of interest) is explicit in some analysis of simulated propagation using changes of scale between the source and observation planes [33, 34]. In our own simulations of Gaussian beam propagation [20, 37, 38], we typically constrain the beam diameter to half the domain of simulation in each $x-$ and $y-$direction in order to avoid edge aliasing effects. In order to present our results in a simple fashion, we assume that most users would have a region of interest defined by approximately this inner portion of the simulation domain.

Defining the measured $x-$direction structure function along the $M/2^{th}$ row from the $M/4^{th}$ point to the $M/4+j^{th}$ point as $D_{x}(j\Delta x)$, and $y-$direction structure function along the $M/2^{th}$ column for its corresponding points as $D_{y}(l\Delta y)$, we can define our percent root mean square (RMS) error metric, $\mathcal{E}$, in terms of the $D_{\theta}$ defined by Eq. 12 via the equations:

[TABLE]

[TABLE]

[TABLE]

It should also be noted that as part of this study, the diagonal direction structure function was also assessed, with similar results. However, because the grid diagonals are not orthogonal to the $x-$ and $y-$directions those metrics are not included in our overall statistics.

We have found that for the range of outer scale values from one to one thousand times the domain of simulation, the RMS error as a percent assessed over half the simulation domain is constrained to less than 4%. More precisely, errors observed over the range $1\leq L_{0}/(M\Delta x)\leq 10^{3}$ ranged from 0.34% - 3.79%. For the unrandomized grid, errors ranged from 2.57% - 61.51% over the same region. Fig. 4 displays the Monte-Carlo simulation results over 25,000 complex phase screens. All data in this study was collected using MATLAB. Because each complex screen contains a real and imaginary component, and structure function is computed over orthogonal $x-$ and $y-$directions, this simulation set contains 100,000 independent samples per point. Results are not shown for $512\times 512$ or $1024\times 1024$ traditional grids due to overlap of the plotted results, i.e. the results were largely indiscernible from the $2048\times 2048$ traditional grid results. Additionally, in order to impress a sense of proportionality upon the reader we have included Fig. 5 which is parametrized by the grid size, outer scale, inner scale, and effective coherence length, $\rho_{o}$, given for isotropic turbulence and Kolmogorov 11/3rds spectral power laws as [19] $\rho_{o}=\left(1.46k^{2}\Delta zC_{n}^{2}\right)^{-3/5}$. On the logarithmic scale, the randomized method follows the theoretical structure function very closely relative to the traditional method.

Although the structure function is typically the key metric of concern when evaluating phase screen algorithms, the overall probability density function of the resultant screens may be of interest to parties concerned about the overall accuracy of wave optics simulations using the proposed method. In order to evaluate this, we have run Monte Carlo trials of our proposed method while logging several key parameters:

Samples of $\theta_{R}(\frac{M}{2},\frac{M}{2})$, the absolute phase of the resultant screens at a central point of the screen. 2. 2.

Samples of $\theta_{R}(\frac{3M}{4},\frac{M}{2})$ - $\theta_{R}(\frac{M}{4},\frac{M}{2})$, the phase difference across the region of interest in the $x-$direction. 3. 3.

Samples of $\theta_{R}(\frac{M}{2},\frac{3M}{4})$ - $\theta_{R}(\frac{M}{2},\frac{M}{4})$, the phase difference across the region of interest in the $y-$direction.

Histogram analysis of the results show some interesting characteristics of the collected samples. For small $L_{0}$’s, the probability density of the absolute phase was clearly Gaussian in nature, however as $L_{0}$ was increased this property quickly faded. Analysis of the natural logarithm of $|\theta_{R}(\frac{M}{2},\frac{M}{2})|$ showed a clear log-normal characteristic as $L_{0}$ was increased beyond $M\Delta x=10$ m. Although a log-normal characteristic of absolute phase may seem undesirable, we do not believe this to be a significant issue with the randomized method as the phase difference histograms maintain approximately zero mean Gaussian distributed characteristics. Results from the trials for several $L_{0}$’s are shown in Fig. 6. We note that the relative phase statistics in the $x-$ and $y-$directions were combined to create the plots on the right hand side of Fig. 6. Repeating this exercise over several different grid sizes, as well as different separations over which the phase differences were measured, yielded qualitatively similar results. Since the statistical fluctuations of interest (e.g. scintillation, beam wander, etc.) the of light undergoing atmospheric turbulence distortion are induced by relative phase, we do not believe that the log-normal characteristic of the absolute phase measured will harm the overall statistical properties of the simulated light provided the relative phase maintains zero mean Gaussian characteristics.

In order to demonstrate that the refinement in structure function accuracy using the randomized method reliably improves statistics of propagating light, we have performed wave optics simulations quantifying the angle of arrival (AoA) fluctuations of plane waves propagating through optical turbulence using both our randomized method and the traditional FFT-based method. Closely following the methodology used by Voelz [23] we have simulated 500 nm plane waves propagating through multiple phase screens, using a number of Rytov variances, and collected the aperture averaged AoA’s for each trial by focusing the light collected over varying aperture sizes and examining the location of the centroid in the focal plane. We have deviated from Voelz’s method, however, in that instead of fixing the relationship between the propagating plane wave’s spatial domain size and aperture diameter, $D$, such that $D/(M\Delta x)=.4$, we have always chosen a spatial domain of $M\Delta x=1$ m prior to aperturing. Once the propagating light reaches the aperture plane, a focusing transmittance function [39] is applied followed by circular aperture functions of several diameters ranging from .5 to 50 cm. The resulting waves from each aperture function are separately propagated to the focal plane using the angular spectrum propagation method discussed by Schmidt [19] in order to allow for differing spatial domain sizes between the aperture and focal planes. For all simulation results and theoretical curves shown, the modified spectrum with inner scale $l_{0}=1$ cm and outer scale $L_{0}=100$ m was used.

A theoretical expression for the variance of AoA for apertured plane waves is given by Cheon [40] as:

[TABLE]

where $f=\sqrt{\lambda L}$ is the Fresnel length, and $A(x)=\left(\frac{2J_{1}(x)}{x}\right)^{2}$ is the Airy function with $J_{1}$ denoting the first-order Bessel function of the first kind. We note that we believe there may have been a typo in [40], and have added a factor of $\pi$ to Eq. 5 in the reference in order to make the equation compatible with later derivations therein. This theoretical metric has been compared to the results of the wave optics simulations in Fig. 7. In general, we find a marked improvement of results relative to theory when the randomized method is in use for all Rytov variances and aperture sizes.

3 Hybrid method for use with bounded and unbounded spectral models

3.1 Core Algorithm

The modified algorithm discussed in the previous section was first investigated with regards to unbounded, anisotropic, non-Kolmogorov spectral models [38]. Therein it was discovered that for structure function power laws greater than the $\nicefrac{{2}}{{3}}$’s of Kolmogorov the randomized algorithm alone was not sufficient to ensure accurate statistics of observed simulated structure function. For this reason, we have developed an algorithm utilizing both FFT-based frequency sampling randomization and subharmonic frequency sampling randomization. We define the following:

[TABLE]

[TABLE]

[TABLE]

In Eq. 19 $N_{p}$ is the number of subharmonic constellations of sampled frequencies (groups of 8 subharmonics chosen from common subgrid boundaries) and $\boldsymbol{\delta}[n,m]$ is the two dimensional discrete Dirac delta function ($\boldsymbol{\delta}[n,m]=1$ for $n=m=0$, otherwise $\boldsymbol{\delta}[n,m]=0$) which we use to ignore the DFT frequency domain origin and the central point of each constellation. It is very important to note that in Eqs. 18 and 19, we choose a different $\delta\kappa_{x}$, $\delta\kappa_{y}$ for each element of the summation (see Fig. 8). That is, for any index ($n$, $m$) change in in Eq. 18 or 19 we choose a new $\delta\kappa_{x}$, $\delta\kappa_{y}$ according to a uniform distribution.

The final hybrid phase screen, $\theta_{H}$, is given by:

[TABLE]

In Eq. 20, the summation over $\theta_{out}$ represents the contributions of each of the $N_{p}$ subharmonic constellations to the phase screen, and $\theta_{in}$ provides the final low-frequency contribution to the phase screen from a spectral sample closest to the origin in $\kappa-$space. The sampling approach described by Eq. 20 can be visualized by Fig. 8 for $N_{p}=1$. In Fig. 8, the red dot closest to the origin represents the $\theta_{in}$ spectral sample, and the red dots in the surrounding eight grid partitions represent the $\theta_{out}$ spectral samples associated with the subharmonic constellation. We have found by choosing the correct number of subharmonic constellations, $N_{p}$, Eq. 20 yields very accurate results for any reasonable spectral model. We shall demonstrate results for both bounded and unbounded spectral models later in this section.

3.2 Addition of White Noise to Phase Screens to Support Subresolution Inner Scales

Since the introduction of generalized spectral models by Kon [41], much theoretical work has gone into the study of wave propagation through anisotropic, non-Kolmogorov turbulence defined by unbounded refractive index spectra [42, 43, 44, 45, 22, 46, 47]. These spectral models are derived from the structure function of refractive index of the form:

[TABLE]

where $\mu_{x}$, $\mu_{y}$ are the anisotropy parameters in the $x-$ and $y-$directions, respectively, $\alpha$ is the three dimensional spectral power law, and $\tilde{C}_{n}^{2}$ is the generalized refractive index structure function constant with units $m^{3-\alpha}$. Although occasionally studies state this structure function model is only valid for $l_{0}\ll\sqrt{x^{2}+y^{2}+z^{2}}\ll L_{0}$ [44], as applied to integrals for calculating second and fourth order beam statistics the inner and outer scales appear as zero and infinity, respectively. In order for the phase structure function integral definitions of Eqs. 11 and 12 to converge $\alpha$ is typically limited to the range $3<\alpha<4$. Also note that for $\alpha=\nicefrac{{11}}{{3}}$, Eq. 21 simplifies to the $\nicefrac{{2}}{{3}}$’s law of Kolmogorov. It can be shown [42] that Eq. 21 corresponds to a three dimensional energy spectrum:

[TABLE]

[TABLE]

where $\Gamma$ denotes the gamma function.

These spectral models do not address practical matters of maximum feature sizes (outer scales) or Kolmogorov microscales (inner scales), where the internal subrange ends and dissipation is the primary form of energy transfer [48, 49]. However, these models are useful for studies of non-classical turbulence when the inertial subrange can be approximated as infinite. For these unbounded cases, great attention has so far been devoted to modeling low spatial frequency components. In order to explain why this is necessary, we note that insertion of the energy spectrum $\Phi_{n}(\kappa_{x},\kappa_{y},0)$ given by Eq. 22 into the structure function identity given in Eq. 11 results in an integrand which diverges as $|\vec{\kappa_{\bot}}|$ approaches zero. For a more thorough explanation of issues involving use of spectral models which diverge at the zero-frequency point, please see [38].

Very little emphasis, however, has been placed on high frequency components outside of the simulated $\kappa$-space. As demonstrated in Figs. 9 and 10, this results in a sag of the phase screen structure function relative to theory over small distances. In order to explain why this sag occurs, we note that the DFT formation of the phase screen algorithm given by Eqs. 3, 5, and 17 include spectral energy contributions only within its frequency domain sample space, which we later define explicitly as $K_{in}$. Spectral energy outside this sample space, $K_{out}$, is not typically included in the discretized simulations, although the structure function integral definition given by Eq. 11 includes spectral energy across all $\kappa$-space. We attempt to resolve this problem via addition of white noise to the screen, in order to simulate spectral energy not included in the $\kappa$-space sampling grid or subharmonic subgrids. Recalling the formation of the structure function in Eq. 11, we calculate the variance of the white noise to be added to the screen as per the set of area integrals:

[TABLE]

[TABLE]

where $K_{out}$ represents the region spanning all of the $\kappa_{z}=0$ plane, which we define unambiguously via:

[TABLE]

In practice, the variances of Eqs. 24 and 25 can be evaluated numerically as the sum of several integrals. For the data sets in this article, four integrals per parameter set were used spanning from the each corner of $K_{in}$ to a $|\kappa_{x}|,|\kappa_{y}|=\infty$ point in an adjacent quadrant of $\kappa-$space. Finally, the variances of two white noise processes are calculated as:

[TABLE]

In order to whiten our phase screens, a $M\times M$ matrix of white noise generated variance $\sigma_{1}^{2}$ is added to the screen, and is followed by addition of a random number with variance $\sigma_{1}^{2}$ across each column (if $\sigma_{x}^{2}>\sigma_{y}^{2}$) or row (if $\sigma_{y}^{2}>\sigma_{x}^{2}$) of the grid. For isotropic turbulence $\sigma_{x}^{2}=\sigma_{y}^{2}$, $\sigma_{2}^{2}=0$, and the second step can be negated. This method ensures the small scale structure function across the $x-$, $y-$, and diagonal directions is improved relative to theory. All random elements pertaining to white noise are generated using a zero mean Gaussian distribution.

Figures 9 and 10 show qualitative results of using this method. As lower power law $\alpha$ values place a higher portion of their spectral energy at high frequencies, we have chosen to display a power law of $\alpha=3.1$. In order to reduce the number of independent variables specifying each plot, Figure 9 has been given in terms of effective coherence length, $\rho_{o}$, which for non-Kolmogorov, anisotropic turbulence is given by [45]:

[TABLE]

Additionally, we have included a curve in Fig. 10 showing the phase screen structure function, normalized by theory, produced using the hybrid method but without the addition of white noise. This is intended to demonstrate the positive effects of our white noise algorithm on the relative error at small separations (i.e. near 1% of the spatial domain and below).

3.3 Results for Bounded Spectral Models

Returning to the modified atmospheric spectrum discussed in Section 22.2, we observe a marked difference in accuracy of the hybrid method vs the subharmonic method of Frehlich [14], which we refer to interchangeably as the traditional subharmonic method. We have chosen to compare with this specific subharmonic method, as opposed to other candidates [13, 50] due to its improved convergence with theory [14] by virtue of weighting the subharmonic amplitude variances using area integrals of spectral models of interest, as opposed to (non-randomized) spectral samplings. Results for several values of the outer scale, $L_{0}$, are shown in Fig. 11, with full results given in Tables 1 and 2. For each case, the size of the outer scale has been set as a factor of the total simulated $x-$, $y-$domain, which was always one meter for this simulation set (i.e. $M\Delta x=1$ m for all bounded spectrum data sets). Because it is impractical to present data for all grid sizes, outer scales, and number of subharmonic constellations, we have focused on the $2048\times 2048$ grid case. The inner scale for the simulations of this subsection was fixed at $l_{0}=M\Delta x/100$ = 1 cm. Because the inner scale of interest was many times larger than the resolution of the grid we elected not to include the white noise algorithm in the results of this subsection.

The data in this subsection represents statistics taken from a large sampling of phase screens, along the $\frac{M}{2}+1$ ordered row and column of each screen. Because the sampling directions are orthogonal the sample set sizes are, essentially, 10,000 trials. We find that for each case the hybrid method outperforms the traditional subharmonic method, which can be verified by close inspection of Tables 1 and 2. In general, the RMS error over the region of interest in the phase screens can be driven to 3% or below with the addition of enough subharmonics. Again comparing Tables 1 and 2, we find that the error ratio of the hybrid method with the subharmonic method can be as low as 13.3%.

3.4 Results for Unbounded Spectral Models

We also wish to assess the accuracy of our revised algorithm for the cases of generalized anisotropic, non-Kolmogorov turbulence spectra discussed in Section 33.2. Although in the previous Section the theoretical structure function metric was computed using numeric integration, in this case simple closed-form solutions exist. By applying a change of variable to Eq. 11 and we obtain:

[TABLE]

$B(\alpha)$ is defined in terms of a Bessel function integral identity [51]:

[TABLE]

where sec denotes the secant function. Substituting $D_{\theta}(j\Delta x,0)$ and $D_{\theta}(0,l\Delta y)$ for the theoretical error expressions of Eqs. 13 and 14, respectively, allows us to again use Eq. 15 as our error metric.

Due to the combination of our assessment of error as a ratio relative to a theory, as well as the scale invariance [13] of this section’s turbulence models, the specific $\tilde{C}_{n}^{2}$ and domain lengths, $M\Delta x$ and $M\Delta y$, do not affect results. The results are, however, sensitive to the $\alpha$ in use and number of subharmonic constellations. Figure 12 displays a comparison of the hybrid method, including the addition of white noise, versus the Frehlich subharmonic method for several $\alpha$’s over 50,000 independent, $1024\times 1024$ resolution phase screens using numbers of constellations between zero and ten. We note that data was also collected for $\alpha$’s of 3.2, 3.4, 3.6, 3.7, and 3.8 up to $N_{p}=10$ however that data is not plotted. Table 3 summarizes the minimum errors observed in testing for both schemes, as well as associated $N_{p}$. We have observed that our hybrid method outperforms the Frehlich subharmonic method on the $1024\times 1024$ grid for any number of subharmonic constellations, except for $\alpha=3.9$. Noting that our revised method falls short for the $\alpha=3.9$ case (at least for $N_{p}\leq 10$), comparing the minimum $\mathcal{E}$ observed using each method for $\alpha=3.1$ to $3.8$ we note that the average ratio of our hybrid method’s $\mathcal{E}$ to that of the Frehlich method is 7.4%.

We have also included results for a $2048\times 2048$ grid case, as shown in Fig. 13 and in Table 4. These metrics were gathered using a smaller number of total phase screens (5,000) due to the longer computation times associated with the larger grid sizes. For $\alpha\geq\nicefrac{{11}}{{3}}$ we have included results up to $N_{p}=20$ in order to demonstrate that for higher $\alpha$’s performance improvements appear to continue as subharmonic constellations above 10 are added. Results for select $\alpha$’s are shown in Figure 13, with We should note that due to the smaller number of total phase screens assessed, the statistical trends appear noisier than the $1024\times 1024$ grid case. Additionally, due to the overall lower error numbers for $\alpha<3.8$ collected using the $1024\times 1024$ grid over a greater number of sampled, we believe it is good assumption the higher resolution $2048\times 2048$ grid results would be improved with a larger sampling phase screens under test. For the $2048\times 2048$ grid case, the average ratio of the randomized hybrid method’s $\mathcal{E}$ to that of the traditional subharmonic method is 15.1% across all $\alpha$’s under test.

3.5 Effect of Direct Fourier Series on Computation Time

In the previous two subsections we have demonstrated that the addition of subharmonics to the FFT-based screens is necessary to achieve accurate RMS error statistics for bounded spectral models with large outer scales, or unbounded spectral models with high power laws. Although we have parameterized the simulation results so far using the number of subharmonic constellations in use, $N_{p}$, we note that the subharmonic modification can also be thought of as a non-uniform, discrete Fourier series with $N_{f}=8\cdot N_{p}$ frequency elements added to each FFT-based screen. We have chosen to introduce $N_{f}$ in order to facilitate a comparison with methods using direct non-uniform Fourier series later in this subsection.

When using subharmonic algorithms, it is quite noticeable to users that computation times are increased relative to the use of a purely FFT-based phase screen algorithm. In order to quantify the computation time effects, we have designed a simulation routine where in-the-loop computational delays are measured in a cumulative fashion over all trials. By in-the-loop delays, we are referring to computational steps which have to, by necessity, be performed each trial of the simulation, such as the creation of the random phase screens from trial to trial. Out-of-the-loop delays, on the other hand, do not need to be performed from trial to trial and in most cases do not have a considerable effect on computation time, and as such are not measured. For phase screen simulations which begin with a propagation of an undistorted Gaussian beam, as an example, we typically perform the propagation to the first screen just once outside of the Monte Carlo loop, meaning it is an out-of-the-loop delay.

Although wave optics simulation parameters will vary depending on application, we have chosen nominal settings of ten screens and ten propagation stages for these measurements as these settings are somewhat typical for our purposes. Simulations were performed in MATLAB using a Windows 10 computer equipped with a 3.6 GHz, 64-bit processor and 32 GB of RAM. These tests utilized the unbounded Kolmogorov spectral model using twenty five subharmonic constellations ($N_{p}$). One hundred simulation trials were run using grid sizes of 256 $\times$ 256, 512 $\times$ 512, 1024 $\times$ 1024, 2048 $\times$ 2048, and 4096 $\times$ 4096, and the various steps of the wave optics simulations were cumulatively timed for later comparison. Each trial consisted of the creation of ten phase screens, phase modulation of the propagating wave using each phase screen, and ten propagations using the FFT-based angular spectrum propagation method [19]. For these computations, one real-valued phase screen derived from each complex screen was used in the simulation. As the angular spectrum propagation method is designed for use with expanding and contracting spatial domains, we find that this is a common practice because two real-valued screens derived from one complex screen cannot be applied over differing spatial domains without additional interpolation steps.

Fig. 14 documents the results of the the simulations. Subplot 14a plots the growth in computation time necessary to produce each phase screen as a function of $N_{f}$ normalized by the time to produce a purely FFT-based screen, and averaged over one thousand screens. Comparing the first and last points of the curves indicates rates of computation time increase ranging from 8.84% per direct Fourier series component for the 256 $\times$ 256 grid to 5.94% for the 4096 $\times$ 4096 grid. Subplot 14b places the results in the context of the overall computation time associated with a full wave optics simulation. This subplot plots the growth in computation time associated with additional direct Fourier series components, as in subplot 14a, but is normalized by the total in-the-loop computation time of the wave optics simulation, which includes phase modulation of the simulated field with the phase screen distortion and propagation using the angular spectrum propagation algorithm. The increases in total computation time ranged from 4.47% per direct Fourier series component for the 256 $\times$ 256 grid to 2.70% for the 4096 $\times$ 4096 grid. Finally, Subplot 14c shows the absolute times associated with each ten screen and ten propagation trial of the wave optics simulation, averaged over 100 trials. We note that because the addition of white noise to the screens discussed in Subsection 33.2 may not be preferable for all spectral models we have not included this factor in the results shown, however we add that including this step increases the time associated with creation of the FFT-based screens by approximately 20% for all grid sizes.

As previously mentioned, the metrics presented in this subsection have been parameterized by $N_{f}$, the number of direct Fourier series components added to the screen using the subharmonic method, for the purpose of comparison with other methods. Charnotskii [17] has previously designed methods for creation of atmospheric turbulence phase screens using randomized spectral sampling, which have been shown to produce phase screens with accurate statistics provided a sufficient number of non-uniform Fourier series elements are used [52]. We note, however, that analysis for this method was performed using different spectral models than those used in this work, and also that RMS averaged or integrated statistics were not provided as functions of power law or inner-to-outer scale ratio. Additionally, the results of simulations in [17] and [52] are compared to asymptotic approximations of the structure function at large and small scales, as opposed to the numerically computed (for bounded spectral models) or closed form solutions (for unbounded spectral models) as we have performed herein. These factors make a direct comparison between these methods rather difficult. We note in general, however, that in order to accurately approximate the structure function metrics of interest such that there is near overlap when plotted on the logarithmic scale a $N_{f}$ of 1,000 to 5,000 must be used.

For the grid sizes which we have analyzed in this subsection, this number of discrete Fourier elements implies a significant increase in computation time relative relative to our method. According to our analysis, the computation time difference associated with $N_{f}=1,000$ ranges from 27.0 times larger than that of the purely FFT-based simulation for the 4096 $\times$ 4096 grid, to 44.7 to times larger than that of the purely FFT-based simulation for the 256 $\times$ 256 grid. Although our randomized algorithm without the use of subharmonics may produce acceptable results for many purposes, in Subsections 32.2 and 33.4 we have demonstrated that our hybrid algorithm produces very accurate results for $N_{p}\leq 10$, $N_{f}\leq 80$ for all but one case, which was for non-Kolmogorov turbulence with $\alpha=3.9$. Comparing the difference between the total delay associated with a wave optics simulation using the $N_{f}=80$ hybrid screens (including FFT-based component) and the extrapolated delay associated with a non-FFT $N_{f}=1,000$ screen predicts that the latter method would result in delays 7.84 times higher for the 4096 $\times$ 4096 grid, and 9.04 times higher for the 256 $\times$ 256 grid. Obviously, we assume that the delays associated with $N_{f}=5,000$ would be five times higher than for the $N_{f}=1,000$ cases, which may lead to unacceptable computation time delays for larger grid sizes.

4 Conclusion and Discussion

In this work we have demonstrated a straightforward method to correct low spatial frequency and periodicity issues associated with FFT-based phase screen methods. We have further demonstrated that low and high spatial frequency contributions can be further improved by combining our core technique with supplemental methods. In previous work [38] we have estimated that the addition of every subharmonic constellations adds approximately the same computation time of generating the purely FFT-based phase screen. As in many cases the core algorithm sans subharmonics outperforms other methods utilizing multiple subharmonics, significant computation time associated with creating the screens may be conserved while garnering improved structure function accuracy. In cases where accuracy is paramount, our hybrid method is able to produce very low error percentages across the region of interest.

Aside from the stated application of simulation of atmospheric optical turbulence, utilizing randomized spectral sampling in concert with the FFT may have additional applications. When modelling processes containing divergences or nulls in their spectral representations we believe a modification of our algorithm giving in Section 22.1 should be considered. Simulations of processes approximated by fractional differencing [53], such as radio frequency oscillator phase noise [54], are a potential candidate. The $1/f^{\ \alpha}$ spectral model (with $f$ the time frequency) of the stochastic process is similar to our own. Additionally, straightforward applications would be to apply the techniques outlined in this study to simulations of partially coherent sources [55, 56], simulations using three-dimensionally correlated phase screens [57, 58], and simulations of optical propagation in underwater turbulence, which has its own unique spectral representations [59, 60, 61]. The heightened low spatial frequency accuracy our techniques will have a direct impact on beam wander statistics, which have been shown as an important factor in free space optical (FSO) communications [62, 63, 64]. Additionally, temporal statistics are often at the center of studies regarding FSO performance metrics [65, 66, 67, 64], which has lead to research on long phase screens to simulate aperiodic turbulence moving at the wind speed [62, 68, 69]. We note that through a combination of circularly shifting the $C$ matrix of Eq. 6 and applying multiplying the exponential term in Eq. 7 by a factor, we have successfully demonstrated the creation of aperiodic moving phase screens derived from the original components of $\theta_{R}$, but not requiring the use of additional FFT’s after the first complex screen realization. Extending this capability to include subharmonic constellations is not challenging. This feature can be directly applied to simulations of time-domain turbulence affects, or combined with other techniques [69, 70] to improve statistics.

5 Funding Information

This work was supported by the Office of Naval Research (ONR) Atmospheric Propagation Studies for High Energy Lasers (APSHEL) program under grant N000141812008.