Free-space optical communication^[1–4], quantum key distribution^[5–9], remote sensing^[10], and astronomical observation^[11,12] are important applications of atmospheric optics. However, the disturbance of atmospheric turbulence has become one of the main restrictions of these applications^[13–15]. The fluctuation of the refractive index leads to wavefront distortion, larger point spread function, random drift of the beam centroid, and intensity scintillation^[16]. To develop correction techniques for compensating its disturbance^[17–22], it is necessary to study the nature of atmospheric turbulence.

In theory, Kolmogorov’s law gives the power spectrum of atmospheric turbulence through basic assumptions^[23], and it has been widely used in turbulence modeling. Further, Taylor gives the frozen-flow hypothesis in the presence of crosswinds^[24] to describe the time evolution of atmospheric turbulence. To allow the research on atmospheric turbulence characteristics and the assessment of compensation technique performance in a laboratory setting^[25,26], the numerical simulation based on the phase screen (PS) transmission model is employed as an effective method for providing controlled and reproducible turbulence^[27]. The key to the simulation is to generate PSs that accurately represent the wavefront distortion caused by atmospheric turbulence.

Several methods have been proposed in the research of simulation over decades, including the spectral sampling method^[28] and the Zernike polynomial method^[29]. For time evolution, there are the long-frame PS translation method^[30], the Markov process-based frequency decay method^[31], and the coefficient evolution method based on the basis vector decomposition^[32]. However, there are some limitations to these methods, such as the error of spatial large-scale structures in the spectral method, high storage demand, and oversimplified time evolution in the translation method. Some of these shortcomings have been overcome by some of the latest methods, such as accurate correction of low-frequency components in the spectral method^[33] and the extrapolation algorithm to reduce memory requirements^[34]. Nevertheless, the statistical characteristics of time-evolving PSs remain unsatisfactory, necessitating heavy computational cost to enhance accuracy^[35], and these methods lack support in the underlying turbulence theory and are not fully validated against models^[36]. For more realistic and high-speed numerical simulation, more accurate and efficient time-evolving turbulence PSs are important and necessary.

We propose a time-evolving turbulence PS model by analyzing the intrinsic fractal nature of atmospheric turbulence. It has a good performance of accuracy, efficiency, and flexibility. First, we derive that scale transformation can adjust turbulence parameters flexibly. Second, we propose an interpolation algorithm for time-evolving PS and improve the spatial interpolation algorithm. Third, combining the covariance matrix^[37] and temporal spectra^[32], we propose an error-free time-evolving PS generation method on the grid. Then, using interpolation algorithms, the PS with low error and high spatio-temporal resolution can be obtained, whose computational overhead is proportional to the number of sampling points. Based on the time evolution under the ensemble average and efficient algorithm, our method can simulate the influence of atmospheric turbulence more realistically with a high frame rate, which is significant for atmospheric optics applications.

In the PS transmission model, atmospheric turbulence can be regarded as a random medium with an uneven refractive index. Based on the propagation equation of light, when the scale of the random medium is much larger than the wavelength, the effect of a section of the random medium can be approximated to phase modulation and propagation^[38].

We use the notation ϕ(r,t) to represent a two-dimensional PS perpendicular to the optical axis at time t, where r=(x,y) are space coordinates. For a turbulence layer of thickness L, refractive index structure constant CN2, and wavenumber k0=2π/λ, the spatial structure function (SF) of atmospheric turbulence phase is given^[39] by Dϕ(ρ)r:=⟨⟨⟨[ϕ(r,t)−ϕ(r+ρ,t)]2⟩⟩r⟩t=6.88(ρr0)5/3,(1)where ρ is the spatial displacement and its modulus is ρ, ⟨⟨·⟩⟩r represents averaging in the spatial domain, ⟨·⟩t represents averaging in the temporal domain, and r0=1.68(CN2Lk02)−3/5 is the Fried parameter^[39], which represents the coherent length of the wavefront. The spatial isotropy of atmospheric turbulence PS is shown by Eq. (1).

When there is crosswind V in the turbulent layer, the temporal SF can be derived based on the Taylor hypothesis^[24]: Dϕ(τ)t:=⟨⟨⟨[ϕ(r,t)−ϕ(r,t+τ)]2⟩t⟩⟩r=6.88(τt0)5/3,(2)where V=|V| and t0=r0/V represents the coherent time of the wavefront (see the Supplementary Material for details). It is worth mentioning that temporal statistical properties are independent of wind direction and location. The consistency of the SF under the spatial and temporal domains can be seen by the same form and power.

As Lane et al.^[28] figured, the fractal nature of a Kolmogorov PS is a direct consequence of the SF being described by a power law. It means that, regardless of the spatial or temporal scale of viewing, all Kolmogorov PSs look similar and have SFs of the same shape. That is say, the fractal nature of atmospheric turbulent PS can be seen from SF after time or space scaling. {Dϕ(αr)r=α5/3Dϕ(r)rDϕ(βt)t=β5/3Dϕ(t)t,(3)where α and β are the scalar factors of the spatial and temporal domains. Additionally, when turbulence parameters r0 and t0 change, it can be equivalent to the scale change of the spatial and temporal domains. Next, we can use this formula to achieve arbitrary scale transformation.

Under spatio-temporal sampling x=mΔr, y=nΔr, and t=lΔt, there is a set of sampled Kolmogorov PSs Φ(m,n,l):=ϕ(mΔr,nΔr,lΔt) that satisfies Eqs. (1) and (2), in which turbulence parameters are r0 and t0. Then, PSs Φ′ with new parameters r0′ and sample intervals Δr′ still should satisfy Eq. (1). According to Eq. (3), the equation is ⟨⟨[Φ′(m0,n0,l)−Φ′(m0+m,n0+n,l)]2⟩⟩m0,n0=Fr5/3⟨⟨[Φ(m0,n0,l)−Φ(m0+m,n0+n,l)]2⟩⟩m0,n0, where Fr=r0Δr′r0′Δr, and we call it the scalar factor. Same to the spatial domain, for PSs Φ′′ with new parameters t0′ and Δt′, we can get the equation according to Eq. (3): ⟨[Φ′′(m,n,l0)−Φ′′(m,n,l0+l)]2⟩l0=Ft5/3⟨[Φ(m,n,l0)−Φ(m,n,l0+l)]2⟩l0, where Ft=t0Δt′t0′Δt. By multiplying the scalar factor before averaging, we can obtain the solution of Φ′,Φ′′: {Φ′(m,n,l)=Fr5/6Φ(m,n,l)Φ′′(m,n,l)=Ft5/6Φ(m,n,l).(4)

Now we can obtain new PSs with new parameters just by scaling without recalculating, which increases the flexibility of existing PS generation methods. In simulation scenarios where turbulence parameters need to be swept, using Eq. (4) can reduce computational overhead.

Since Kolmogorov PSs are fractal and can be actually described by the fractional Brownian motion (fBm)^[40], an efficient algorithm named random midpoint interpolation^[41] can be used to improve the resolution of the PS at the expense of small error. A one-dimensional fBm BH(t) is a continuous-time Gaussian process characterized by an incremental variance^[42]: ⟨[BH(T)−BH(T+t)]2⟩∝t2H,(5)where H is the Hurst parameter and is in the range of 0–1. The fractal dimension of the fBm is F=E+1−H^[43], where E is the dimension of parameter t. According to Eq. (2), it is clear that the Kolmogorov PS at any r in the temporal domain is an fBm with H=5/6 and F=7/6.

The principle of the random midpoint interpolation is new=avg+scaleH×random, where avg is the average value of adjacent points and random is the Gaussian random number. Thus, we can propose the interpolation algorithm to enhance temporal resolution: For a given time-evolving turbulence PS sequence Φt,in(r,l) with length L and interval Δt, the interpolated PS sequence Φt,out(r,l′) has 2L−1 frames with interval Δt/2. The process is divided into two steps.

First, 2L−1 frames of PS sequence Φt,out are created for assigning each frame of Φt,in(r,l) to the odd l′. Then, the even l′ frame of the PS is computed by $$\begin{gathered} Φ(r,l′)=12[Φ(r,l′−1)+Φ(r,l′+1)]+Ct(Δt)5/6σ(r) \\ =ItΦ+Ct(Δt)5/6σ(r), \end{gathered}$$(6)where It represents the calculation coefficient vector, Ct is a constant, σ(r) is the normalized Gauss random field with zero mean, and its power spectrum is the same as Kolmogorov turbulence’s k−11/3 power spectrum^[44].

For the Kolmogorov PS in the spatial domain, the interpolation in the spatial domain has been proposed by Lane et al.^[37], and it is actually a two-dimensional fBm with H=5/6 and F=13/6according to Eq. (1). Thus, we proposed the modified interpolation algorithm to enhance spatial resolution: the random displacement is proportional to (Δr)5/6, and the modified calculation coefficient matrix is Ir4×4=[−0.0017−0.0341−0.0341−0.0017−0.03410.31990.3199−0.0341−0.03410.31990.3199−0.0341−0.0017−0.0341−0.0341−0.0017].(7)

Here the sum of the calculation coefficient of neighbors is a unit. Otherwise, a proportional cumulative error cannot be negligible after a multi-round iteration of interpolation, especially if the absolute value of the initial phase is large. Figure 1 shows the difference in the PS before and after the modification of the interpolation algorithm. Small-scale structures obtained by interpolation deviate from theory because (Δr)5/3 is smaller than (Δr)5/6 when Δr is less than 1. The chessboard-like error occurs in the interpolated PS after adding the 64π offset phase in the red box area because of the non-normalization of coefficients, and it remains the same under modulo 2π after correction.

Here is the modified interpolation algorithm. Assuming Φin(m,n) is the initial PS with a grid size N×N and Δr, and the interpolated PS Φout(m′,n′) has (2N−3)×(2N−3) grid size with spacing Δr/2, as shown in Fig. 2, the process can be divided into three steps.

First, (2N−1)×(2N−1) PS Φ′(m′,n′) is created for assigning each value of Φ(m,n) to the odd coordinates, which are illustrated by a solid circle. Second, the shaded circles are computed in parallel by the following equation using their nearest neighbors: $$\begin{gathered} Φ′(m′,n′)=14[Φ′(m′−1,n′−1)+Φ′(m′+1,n′−1) \\ +Φ′(m′−1,n′+1)+Φ′(m′+1,n′+1)] \\ +Cr(Δr)5/6σ \\ =IrΦ′+Cr(Δr)5/6σ. \end{gathered}$$(8)

Next, the dotted circles are computed in parallel by the following equation: $$\begin{gathered} Φ′(m′,n′)=14[Φ′(m′−1,n′)+Φ′(m′,n′−1) \\ +Φ′(m′+1,n′)+Φ′(m′,n′+1)] \\ +Cr(Δr/2)5/6σ \\ =I˜rΦ′+Cr(Δr/2)5/6σ, \end{gathered}$$(9)where Ir represents the calculation coefficient matrix. After 45° rotation, it becomes I˜r. Cr is a constant, and σ is a Gaussian random variable with zero mean and unit variance. Last, all edge points of the Φ′(m′,n′) are dropped out to obtain interpolated Φout. This is because, for the finite-size grid, edge points do not have enough neighbors to interpolate. A detailed description of the algorithm can be found in the Supplementary Material.

However, an accurate initial PS is required for applying interpolation. Therefore, based on the covariance matrix^[37] and temporal spectra^[32], we propose an error-free time-evolving turbulence PS generation method on the grid. In short, we can derive time-evolving PS Φ(t)=Uλ(t)=∑iMNλi(t)Ui by decomposition of the phase fluctuation covariance for a given sampling grid M×N, and the coefficient λi(t) can be described by temporal power spectrum density (PSD) wλ,i(ν). Please see the Supplementary Material for details. Although our method has a high computational load and limited time-evolving sequence length, the interpolation algorithm can compensate for these deficiencies effectively. The grid-based decomposition makes our method more accurate than the Zernike polynomial^[32] under the same computing resources.

Through numerical simulation, we confirmed that the Kolmogorov PSs obtained by the generation method and interpolation algorithm are consistent with the theory.

First, we calculated the complete orthogonal bases Ui of the PS for the 8×8 grid, shown in Fig. 3(a). To validate the decomposition and superposition, a set of 104 static PSs are generated, and its SF is calculated by Eq. (1). Unlike time-evolving PSs, these PSs are independent of each other. We use R(ρ=Δrm+n2):=⟨D(m,n)r⟩/Dϕ(ρ)r, the ratio between the calculated value and theoretical value of the SF to show results. For the accurate PS, R(ρ)=1 at arbitrary ρ. As shown in Fig. 4(b), we can see that both small and large scales are in good agreement with the theory for different grids, in which the maximum deviation is less than 2%. The standard deviation of R(ρ) decreases with the root of sampling number N shown in Fig. 4(d), which accords with the central limit theorem. Thus, it is shown that our generation method is theoretically accurate with only statistical error, which is a significant advantage.

Then, by applying spatial interpolation, we can iteratively enhance the spatial resolution of PSs, which is shown in Fig. 4(a). As shown in Fig. 4(c), the deviation of R(ρ) of interpolated PSs remains less than 2%. However, R(ρ) at the nearest neighbor distance is higher than 1, which is caused by the independently distributed random variable σ. More importantly, Fig. 4(e) shows that the ⟨R(ρ)⟩ρ of interpolated PSs are always within an acceptable range, and the fluctuation does not increase with iterations of interpolation. In our practice, the interpolation process can be iterated at least 15 times without appreciable deviation.

For a sufficiently long initial time-evolving PS sequence, we choose the 32×32 grid to get 32 initial frames of time-evolving PSs. Assuming V=5m/s, r0=0.1m, and Δr=0.01m, then t0=r0/V=20ms, Δt=Δr/V=2ms, and we get the temporal PSD and time sequence of the coefficient shown in Figs. 3(b) and 3(c). It can be seen that coefficients of low-order Karhunen–Loeve (K–L) functions only have low-frequency components but a large variance. Conversely, coefficients of high-order K–L functions have full frequency components but a small variance.

By performing temporal interpolation, we can get more frames in the PS sequence. As shown in Figs. 5(a) and 5(b), when 32 initial frames of time-evolving PSs undergo two iterations of temporal interpolation, we can obtain a sequence of 125 frames with Δt=0.5ms. To validate the temporal characteristics, 103 sets of time-evolving PS sequences are generated, and the temporal PSD of phase fluctuation is calculated by wΦ(ν)=⟨⟨|Φ^(m,n,ν)|2⟩⟩m,n, where Φ^(m,n,ν) is the Fourier transform of Φ(m,n,t). Figures 5(c)–5(e) show the comparison between the average PSD ⟨wΦ(ν)⟩ and the theoretical PSD wϕ(ν)∝ν−8/3^[44]. First, it is verified that spatial interpolation has no effect on the temporal characteristics. Second, due to temporal interpolation, the spectrum extends to high frequency, and it is in good agreement with the theoretical curve. All misalignment introduced by interpolation occurs at the Nyquist frequency of the initial sequence and the interpolated sequence. We speculate that the reason is that random midpoint interpolation breaks the periodicity generated by inverse discrete Fourier transform (IDFT). Last, the interpolation simultaneously in the temporal and spatial domains has been verified.

In summary, the directly generated PS has only statistical fluctuations. The spatial interpolation does not introduce perceptible error, and there are no cumulative errors in iterations. The temporal interpolation only introduces a negligible number of high-frequency errors.

The methods for generating turbulence PS are divided into static PSs^[28,37] and time-evolving PSs^[31,32]. Our method can generate both and has advantages in accuracy and efficiency. Next, we compare it with existing methods to verify the advantages of the method proposed in this work. Please see the Supplementary Information for details.

We use the IDFT method based on the spatial power spectrum^[28] and the sub-harmonic method corrected by Frehlich^[27] to generate 1000 static PSs whose size is 512×512, Δr=2mm, and r0=0.1316m. As shown in Fig. 6(a), our method is consistent with theory at both the small and large space scales compared with others. Then, we perform spatial shifting and frequency decay on static PSs to obtain 1000 frames of the time-evolving PS^[31] in which Δt=1ms. Here, the crosswind is set to 5 m/s, and the coherence time t0 is 26.32 ms. The results and theoretical curve are shown in Fig. 6(b). It is clear that our method keeps errors with the theory within the smallest range compared with others, and the spatial shifting reveals the periodic error caused by fast Fourier transform (FFT)-based methods. Besides, related works^[45,46] demonstrated that the Zernike polynomial method had no advantage over FFT-based methods in terms of accuracy and efficiency. Compared with the temporal spectra method^[32], which is one of the recent representative works, our method also shows advantages in accuracy both in the spatial and temporal domains.

Then, we compared our simulation results with experimental data^[47] in the 460 m field channel. The angle of arrival (AoA) temporal power spectrum of field-test results is shown in Fig. 7(a), in which λ=1550nm, the beam radius w=30mm, and the receiver diameter D=80mm. The measured turbulence parameter is r0=3.7cm. Our simulation results are shown in Fig. 7(b) under the same conditions, and the theoretical −2/3 power asymptote^[44] is also plotted. The results demonstrate that the simulation is well consistent with the real turbulence in the real-world, which verifies the effectiveness of our method. Compared to turbulence in the field channel, the PS model can provide reproducible parameter-controlled turbulence, which is useful for developing compensation techniques and evaluating system performance^[47]. In future work, it is worth enhancing the performance and expanding the application of our model through extensive experimental analysis.

In this work, we first achieve scale transformation based on a fractal model for time-evolving turbulence PSs, which facilitates parameter scanning in simulations for various scenarios. Second, the random midpoint interpolation algorithm is proposed to enhance the spatio-temporal resolution of PS efficiently. Further, we propose an error-free time-evolving PS generation method on the grid, and then the errors and time cost of the generated PS and interpolation algorithm are analyzed quantitatively, which proves that our method has remarkable accuracy and efficiency.

Additionally, the evolution based on temporal spectra is more realistic compared with the long-frame translation because it means that the time interval Δt is equivalent to translating the average distance of VΔt in the ensemble, and each sampled frame actually has fluctuations in series. Furthermore, the fractal nature of the atmospheric turbulence PS and random midpoint interpolation algorithm have wide applicability, which could greatly reduce the workload in PS simulation.