1. Introduction

Atmospheric turbulence imposes a fundamental limit to the resolution of ground based telescopes [1, 2]. However it is possible to overcome this limitation thanks to Adaptive Optics (AO) [3], which allows compensation in real-time for this degradation by using a wavefront sensor linked to a deformable mirror.

As AO systems are very complex, they require extensive numerical simulation in order to design the architecture of the system [4]. The simulation of atmospherically distorted wavefronts is the starting point of all these models. Several methods currently exist to create atmospheric phase screens: modal techniques [5] or sample based. The latter can be subdivided into two subcategories: Fast Fourier Transform (FFT) based [6, 7] or covariance based [8].

For computing reasons, it is often preferred to use FFT based methods. These methods can also be used to study the temporal properties of atmospheric turbulence, thus allowing estimation of the temporal bandwidth of AO systems. This is accomplished by translating the generated phase screen over the telescope pupil. However, phase screens generated with these methods generally suffer from periodicity as well as poor representation of low spatial frequencies. A method to overcome this problem has been proposed [9], but is not straightforward to implement.

Moreover, the FFT methods suffer from two problems. First, because of the finite size of the computed phase screen, it is only possible to simulate finite exposure times, which become shorter as the simulated wind velocity increases. This is an acute problem when simulating very large apertures (more than 20 meters), for which the pupil has to be discretised over a large number of pixels (typically 512). As an example, if a 30 meter telescope is simulated over a grid of 1024×1024 pixels, and that we want to simulate a 1 mn exposure time, for a typical wind velocity of 10 m/s, the resulting phase screen computed with a FFT will have a size of 20480×20480 pixels, resulting in a memory usage of more than 3Gb for 64 bits float precision. Secondly, the FFT generated phase screens have static statistics, i.e. the Fried parameter r ₀ as well as the outer scale L ₀ are the same across the entire screen. In fact, for the studies of the AO systems for next generation of Extremely Large Telescopes (ELTs), applications such as exoplanet detection or long-exposure spectroscopy will require simulations of non-stationary turbulence, which are non possible with current methods.

To overcome these problems, we propose in this paper a new method to generate phase screens for very long exposures as well as for the simulation of non-stationary turbulence effects. But contrary to the existing methods, we do not need to store in memory the phase screen for the entire simulation. Indeed the goal is to use already computed phase values to generate new phase values consistent with the previous ones. Memory storage is required only for the phase screen needed at the specific simulation iteration.

2. Principle of the method

Let φ(u,v) be a square discrete phase screen computed via the FFT based method, with a size of d_p × d_p pixels. We have therefore 1 ≤ u,v ≤ d_p. Moreover we assume that these d_p pixels correspond to a physical size of D meters. The pixel scale is therefore p_S = D/d_p meters/pixel.

Now let φ the vector storing all the elements of the phase screen, and Z the vector storing the phase values of the last N_col columns. Therefore Z has d_p × N_col elements. Our goal here is to extend the phase screen, i.e. add a new column to it. Let us then call X the vector storing these new phase values. We assume the following relation between X and Z:

where A and B are two matrices, and β is a random vector with d_p independent random variables. Moreover, as the turbulent phase is zero mean and has Gaussian statistics [2], we assume that β is therefore a gaussian vector with zero mean and unit variance: 〈β〉 = 0, and 〈ββ^T〉 is equal to the identity matrix.

We now multiply Eq. 1 by Z ^T, and then take the average of the resulting equation. As β and Z are uncorrelated, we find the following relation:

thus allowing us to compute A:

We focus now on the matrices 〈ZZ ^T〉, 〈XZ ^T〉, and 〈XX ^T〉, which give the covariance of the phase for the X and Z vectors. We assume in the rest of this paper that atmospheric turbulence follows the Von-Kármán model, which introduces the wave-front outer scale L ₀, and which agrees with the latest turbulence measurements [10]. The phase structure function then has the following expression [11]:

where K _5/6 (r) is the modified Bessel function of the third kind, or McDonald function, Γ(x) is the gamma function, and r ₀ is the Fried parameter [12]. The phase structure function D_φ(r) is linked to the phase covariance function C _φ(r) = 〈φ(x)φ(x + r)) by the relation D_φ(r) = 2 [σ_φ ² - C_φ(r)], σ_φ ² being the phase variance. We have therefore:

Then, if (u_i,v_i) and (u_j,v_j) are the respective coordinates of the points i and j, and if we call r_ij = [(u_j - u_i)² + (v_j - v_i)²]^1/2 × p_s the physical distance between these two points, then the 〈Z_iZ^T_j〉, 〈X_iZ^T_j〉 and 〈X_iX^T_j〉 terms will be given by the value of the phase covariance function C_φ(r) at the distance r_ij.

Equation 1 gives also the following relation:

By taking the average of this equation, and as β and Z are uncorrelated, we have therefore:

and as 〈ZZ ^T〉 is a symmetric matrix, its inverse is also symmetric. Then, we can write

and we obtain:

We therefore have the expression of the BB ^T matrix. However we need the B matrix for our method. As BB ^T is symmetric, a singular value decomposition [13] allows to write it in the form BB ^T = UWU ^T, the columns of U being the eigenvectors of BB ^T and W being a diagonal matrix with the eigenvalues λ_i. Moreover, as W is also diagonal, it can be written as W = LL ^T , where L is a diagonal matrix, and whose diagonal values are the square roots of the λ_i. We can then write

thus B = UL.

We therefore know how to compute the A and B matrices, which can then be used to extend the phase screen. However the dimension of these matrices is directly linked to the number of previous columns N_col used to compute new phase values as well as the number of pixels used to sample the telescope pupil. We therefore focus in the next section on the number of columns needed in order to optimize the size of these matrices.

3. Maximum distance to take into account

The former section gave the mathematical formalism of our method of phase screen extension. We focus here on its practical implementation.

Equation 1 shows that our method uses two matrix products to compute the new phase values. The first matrix product is the multiplication of the A matrix by the Z vector. As Z has d_p × N_col elements, this implies that A is a d_p × (d_p × N_col) matrix. Therefore, the size of the A matrix evolves with N_col, implying an increase of the computation time when new phase values are computed.

To establish the constrainsts on the value of N_col, we performed several simulations, for different values of D, L ₀, and N_col. We used two methods to check that the phase screens computed have the correct statistics: the phase structure function as defined by Eq. 4, as well as the projection of the phase on the Zernike polynomials Z_i(2 ≤ i ≤ 66), in order to compute the variance of their coefficients a_i. The latter can be derived from the diagonal of the covariance matrix 〈a_ia_j〉, whose analytical expression is given in [14] for a turbulence with finite outer scale.

We started by simulating a D = 8 meter diameter pupil with a number of pixels d_p = 64, and we assumed an outer scale L ₀ of 16 meters (L ₀/D = 2). The ratio of (D/r ₀) was equal to 1, and was arbitrary as this ratio acts only as a scaling factor for the produced phase screen. The initial Von-Kármán phase screen was computed via the FFT based method. Then we used Eq. 1 to add one new column to the right of the phase screen and to remove the first column, so that the generated phase screen always has a size in memory of d_p × d_p pixels. This was repeated 15 000 times, so as to add 15 000 new columns and compute 15000 correlated phase screens with a size of d_p × d_p pixels. This gave us a sufficient number of samples allowing us to compute the phase structure function as well as the variance of Zernike polynomial coefficients. Moreover, we considered three different sizes for the vector Z in Eq. 1, corresponding to the case where this vector is storing the phase values of the last N_col columns, N_col being respectively equal to 1, 2 and 4. The result of these simulations is shown on the Fig. 1, where the theoretical and observed phase structure functions as well as Zernike polynomials variance are compared. We can see there that there is a very good agreement between theory and the numerical implementation, as soon as the last N_col = 2 columns of the phase screen are used to compute the new phase values.

We conclude therefore that using the last 2 columns is sufficient to compute new phase values consistent with the already computed ones, allowing us to extend the phase screen as long as required. To examine any dependence on L ₀, we repeated the simulation, this time assuming an outer scale L ₀ of 64 meters (L ₀/D = 8), and considering again two values for N_col: 2 and 4. The results of this study are shown in Fig. 2. Here again, we can see there that there is a perfect agreement between the observed and theoretical phase structure functions and Zernike polynomial variance for N_col = 2.

 

Fig. 1. Phase structure function (left) and variance of the Zernike polynomials coefficients (right) computed from 15000 generated phase screens using Eq. 1, for several values of N_col: 1 (top), 2 (middle), and 4 (bottom). A D = 8 meter pupil was simulated on a square grid of 64×64 pixels, and the outer scale was equal to L ₀ = 16 meters. The broken line on the figures on the left corresponds to the theoretical phase structure function, whereas the full line corresponds to the computed one. The full line on the figures on the right shows the theoretical variance of the coefficients of the Zernike polynomials 2 to 66, and the squares the computed values.

Download Full Size | PDF

 

Fig. 2. Phase structure function (left) and variance of the Zernike polynomials coefficients (right) computed from 15000 generated phase screens using Eq. 1, for two values of N_col: 2 (top) and 4 (bottom). A D = 8 meter pupil was simulated on a square grid of 64×64 pixels, and the outer scale was equal to L ₀ = 64 meters (L ₀/D = 8).

Download Full Size | PDF

This study shows therefore that our method is insensitive to the (D/L ₀) ratio: as long as two columns of the previous phase screen are used to compute a new column in order to extend the phase screen, meaning that we use phase values up to a distance of 0.25 meters, the newly produced phase screen has the correct statistical properties, which is confirmed by the comparison of the phase structure function as well as the Zernike polynomial variance.

Until now, we considered a D = 8 meter pupil simulated on a square grid of d_p × d_p pixels (d_p = 64), leading to a spatial sampling of p_s = 0.125 meters per pixel, and for this spatial sampling, two columns are sufficient to create a new column, independent of the outer scale L ₀ . To explore the requirements for finer spatial sampling, we performed new simulations, but this time considered a D = 4 meters aperture simulated on a grid of 64×64 pixels, leading to a spatial sampling p_s = 0.0625 meters per pixel.

We first considered the same outer scale value as in our previous first simulation case, i.e. L ₀ = 16 meters (L ₀/D = 4), and then used either N_col = 2 or N_col = 4 columns to create the new column of the new phase screen. The comparison of the observed and theoretical phase structure function as well as Zernike polynomial variance is shown on the Fig. 3. This time we see that the observed and theoretical phase structure functions for both values of N_col do not match perfectly, especially for distances greater than D/4, but however are in good agreement with each other. In the case where N_col = 2, the observed phase structure function is greater than the theoretical one for high separations, whereas we observe the opposite effect when N_col = 4. This is in fact confirmed by the Zernike polynomial variance for both cases, and especially for the tip and tilt modes, where the variances are respectively over or underestimated. But the relative error remains very small (less than 10%), so that the phase screens have the correct statistics. This was confirmed when we did the same study for the second value of the outer scale studied previously, i.e. L ₀ = 64 meters. The results are shown in Fig. 4 and we can see there that, again, the use of the last N_col = 2 columns allows us to produce a new column in the phase screen with the correct statistics, as shown by the phase structure function as well as the Zernike polynomial variance.

We have focused in this section on the practical implementation of our new phase screen generation which consists of extending an already computed phase screen by using known values to produce new values consistent with the previous ones. Here we have studied the quality of the resulting phase screen for different values of the outer scale L ₀, the telescope diameter D as well as the number of pixels d_p used to simulate the pupil, as a function of N_col, the number of previous columns of the phase screen used to compute new phase values. We found that the use of the last 2 columns was sufficient to produce an infinitely long phase screen with the correct statistics, independent of all the other parameters.

It should be emphasized here that the phase screens always used d_p × d_p pixels, by removing the first column before adding a new column at each iteration. However, the method can also be used to create rectangular phase screens. This requires initially that the memory space needed by the rectangular array is allocated, before it is filled by extending the initial square phase screen. In this case the array generated can be stored and subsequently used later and repeatedly.

4. Simulation of non-stationary phase screens

We have shown in the previous section that the last 2 columns of a phase screen can be used in order to extend it as much as is required, via Eq. 1. We show in this section how the same equation can be used to extend the phase screen, and at the same time to simulate the evolution of the atmospheric parameters.

We focus in particular on the simulation of a non-uniform Fried parameter r ₀. Equation 3 gives us the definition of the A matrix: A = 〈XZ ^T〉 〈ZZ ^T〉^-1. As discussed in section 2, the generic term of each matrix 〈XZ ^T〉 or 〈ZZ ^T〉^-1 is given by the phase covariance function C_φ(r), whose expression is:

 

Fig. 3. Phase structure function (left) and variance of the Zernike polynomials coefficients (right) computed from 15000 generated phase screens using Eq. 1, for two values of N_col: 2 (top) and 4 (bottom). A D = 4 meter pupil was simulated on a square grid of 64×64 pixels, and the outer scale was equal to L ₀ = 16 meters (L ₀/D = 4)

Download Full Size | PDF

 

Fig. 4. Phase structure function (left) and variance of the Zernike polynomials coefficients (right) computed from 15000 generated phase screens using Eq. 1, for two values of N_col: 2 (top) and 4 (bottom). A D = 4 meter pupil was simulated on a square grid of 64×64 pixels, and the outer scale was this time equal to L ₀ = 64 meters (L ₀/D = 4)

Download Full Size | PDF

Therefore, when the matrix product between 〈XZ ^T〉 and 〈ZZ ^T〉^-1 is performed, the terms (L ₀/r ₀)^5/3 cancel. This means that if we assume a constant outer scale L ₀, A is insensitive to r ₀.

Let us then focus on the B matrix. From Eq. 9, we have:

Therefore, the matrices 〈XX ^T〉 and A〈ZX ^T〉 make use of the phase covariance function C_φ(r), and thus the (L ₀/r ₀)^5/3 term which can be introduced as a factor in the above equation. Therefore, if we compute the BB _T matrix for (L ₀/r ₀) = 1, multiply the resulting B matrix by the appropriate factor (L ₀/r ₀)^5/6 , and then use this matrix in Eq. 1 to create a new column in order to extend the phase screen, the turbulence in this column will be characterized by this r ₀. If we repeat this operation and change the r ₀ at each iteration, then this allows us to create a phase screen with a different r ₀ for each column, and thus a non-stationary phase screen.

A numerical simulation was performed to check the validity of this method, in which we assumed a pupil of 40 cm simulated over 32 pixels, an outer scale L ₀ equal to 4 meters (L ₀/D = 10), and a Fried parameter r ₀ = 15 cm. We first used the Fourier method to create an initial phase screen of 800 × 800 pixels, which was then cropped to 800 × 32 pixels, corresponding to a physical length of 10 meters. We then updated the last column of this phase screen 37200 times, equivalent to a phase screen 480 meters long, assuming the same r ₀ value of 15 cm. Next we repeated the operation, but now we assumed a sinusoidal r ₀ oscillating between 10 and 20 cm, with a spatial period of 240 meters. This was then equivalent to a phase screen with a virtual length of 960 meters (76000 × 32 pixels). Moreover, at each iteration where the phase screen was extended, we projected the phase on the Zernike polynomials 2 to 36, thus providing us a set of 76000 × 35 coefficients.

We then used these Zernike coefficients to compute their variance, in order to measure the Fried parameter r ₀. Indeed, for a given (D/L ₀), the Zernike polynomial variance is only proportional to (D/r ₀)^5/3 [14]. Therefore, we first computed the theoretical variance of the Zernike for (D/r ₀) = 1. Then we performed a least-square fit of the Zernike coefficient variance to the theoretical one so as to compute the (D/r ₀)_5/3 scaling factor, and hence r ₀ at each iteration. The evolution of r ₀ as a function of the iteration number is shown in the top panels of Fig. 5, where we can see that there is an excellent agreement between the theoretical and measured values. In particular the sinusoidal evolution of r ₀ is clearly visible on these figures, showing that we effectively created a non static phase screen. Moreover the error between observation and theory remains very small. Indeed, the histogram of the r ₀ error shows that it can be fitted with a Gaussian distribution, with an average error of 0.34 cm and a standard deviation σ = 0.81 cm.

We argue that this error is mainly due to the finite number of samples used to compute r ₀ from the Zernike polynomials variance. Indeed, the bottom panels on the Fig. 5 show the result of the measure of r ₀ from Zernike polynomials variance, but this time computed on a phase screen with the same virtual length than before (960 meters) and with an uniform r ₀ = 15 cm. In particular, if we look again at the histogram of the error, we find that it can also be fitted by a Gaussian distribution, with properties very similar to the case where r ₀ is not static: we find that the r ₀ error has an average value equal to 0.40 cm, with a standard deviation σ = 0.75 cm.

We have therefore shown in this section that our method allows simulation of non-stationary turbulence effects, contrary to traditional methods which assume the same turbulence conditions all over the computed phase screen. In particular we have shown here that the simulation of a phase screen with a spatially non uniform Fried parameter r ₀ is straightforward to implement. The simulation of a non-static outer scale L ₀ can also be implemented via a similar approach. However it will be more difficult as it requires computation of the A and B matrices at each simulation iteration, as the phase covariance function C_φ(r) does not scale linearly with the outer scale.

 

Fig. 5. Comparison (left), error (middle), and histogram of the error (right) between the the theoretical and observed Fried parameter r ₀. The top panels correspond to a phase screen where r ₀ is equal to 15 cm in the first half of the screen and oscillating between 10 and 20 cm in the second half. The bottom panels correspond to a phase screen with constant r ₀ = 15 cm. The black line on the left figure shows the expected r ₀, whereas the red line shows the r ₀ value deduced from Zernike polynomials variance, both as a function of the iteration number. The dashed blue line on the r ₀ error histogram shows the result of a fit of the error by a Gaussian distribution

Download Full Size | PDF

5. Conclusion

We have demonstrated in this paper a new method to simulate atmospheric turbulence. The method permits simulation of arbitrarily long exposure times by extending a given phase screen by as much as is needed. We have shown that only the last two columns of the phase screen are required to compute a new column, independent of the outer scale or the spatial sampling. This implies that the phase screen for the whole simulation does not need to be stored in memory, but only the phase values required at each simulation iteration for a given array size. This makes the method particularly suitable for the simulation of Adaptive Optics systems for Extremely Large Telescopes, where memory management will be a critical issue because of the large number of points required to simulate the telescope pupil.

Another key improvement provided by this method is the ability to simulate realistic turbulence conditions, i.e. a turbulent phase with properties evolving as a function of time. In particular, It is straightforward to simulate with this method a turbulent phase with a non static Fried parameter r ₀, which is again well suited to studies of AO systems for Extremely Large Telescopes. Indeed very long exposures (more than one hour) are very likely to be performed on such telescopes for applications such as exoplanet detection or the observation of very distant galaxies. In this case a realistic simulation will be required in order to examine in detail the sensitivity of the AO system to the variation of atmospheric parameters.