# Method for simulating infinitely long and non stationary phase screens with optimized memory storage

François Assémat, Richard Wilson, Éric Gendron

► To cite this version:

François Assémat, Richard Wilson, Éric Gendron. Method for simulating infinitely long and non stationary phase screens with optimized memory storage. Optics Express, 2006, 14, pp.988-999. 10.1364/OE.14.000988 . hal-03801614

HAL Id: hal-03801614

https://hal.science/hal-03801614

Submitted on 30 Dec 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Method for simulating infinitely long
and non stationary phase screens with
optimized memory storage

François Assémat and Richard W. Wilson

Centre for Advanced Instrumentation, Department of Physics, Durham University,
South Road, Durham DH1 3LE, United Kingdom
francois.assemat@durham.ac.uk, r.w.wilson@durham.ac.uk

Eric Gendron

Laboratoire d'Etudes Spatiales et d'Instrumentation en Astrophysique, Observatoire de Paris,
5 place Jules Janssen, 92195 Meudon, France
eric.gendron@obspm.fr

**Abstract:** A new method to simulate turbulent phase is investigated in this paper. The goal of this method is to be able to simulate very long exposure times as well as time evolving turbulence conditions. But contrary to existing methods, our method allows to simulate such effects without using the whole memory space required by the simulated exposure time, making it particularly suited to the simulation of adaptive optics systems for very large apertures telescopes.

© 2006 Optical Society of America

OCIS codes: (010.1330) Atmospheric turbulence; (350.5030) Phase.

**References and links**

1. D. L. Fried, "Optical Resolution Through a Randomly Inhomogeneous Medium for Very Long and Very Short Exposures," J. Opt. Soc. Am. **56**, 1372–1379 (1966).

2. F. Roddier, "The Effects of Atmospheric Turbulence in Optical Astronomy," Prog. Optics **19**, 281–376 (1981).

3. F. Roddier, *Adaptive optics in astronomy* (Cambridge University Press, 1999).

4. M. Carbillet, C. Vérinaud, B. Femenía, A. Riccardi, A. and L. Fini, "Modelling astronomical adaptive optics - I. The software package CAOS," Mon. Not. R. Astron. Soc. **356**, 1263–1275 (2005).

5. N. Roddier, "Atmospheric wavefront simulation using Zernike polynomials," Opt. Eng. **29**, 1174–1180 (1990).

6. B. L. McGlamery, "Computer simulation studies of compensation of turbulence degraded images," in *Image processing*, Proc. SPIE, 225–233 (1976).

7. G. Sedmak, "Implementation of Fast-Fourier-Transform-Based Simulations of Extra-Large Atmospheric Phase and Scintillation Screens," Appl. Opt. **43**, 4527–4538 (2004).

8. C. M. Harding, R. A. Johnston and R. G. Lane, "Fast Simulation of a Kolmogorov Phase Screen," Appl. Opt. **38**, 2161–2170 (1999).

9. R. G. Lane, A. Glindemann and J. C. Dainty, "Simulation of a Kolmogorov phase screen," Waves Random Media **2**, 209–224 (1992).

10. A. Ziad, J. Borgnino, F. Martin, J. Maire and D. Mourard, "Towards the monitoring of atmospheric turbulence model," Astron. Astrophys. **414**, 33–36 (2004).

11. A. Tokovinin, "From Differential Image Motion to Seeing," Publ. Astron. Soc. Pac. **114**, 1156–1166 (2002).

12. D. L. Fried, "Statistics of a Geometric Representation of Wavefront Distortion," J. Opt. Soc. Am. **55**, 1427–1435 (1965).

13. W. H. Press, S. A. Teukolsky, W. T. Vetterling and B. P. Flannery, *Numerical recipes in C. The art of scientific computing. 2nd ed.*, (Cambridge University Press, 1992).

14. N. Takato N. and I. Yamaguchi, "Spatial correlation of Zernike phase-expansion coefficients for atmospheric turbulence with finite outer scale," J. Opt. Soc. Am. A **12**, 958–963 (1995).
# 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_0$ as well as the outer scale $L_0$ 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 $\varphi(u, v)$ be a square discrete phase screen computed via the FFT based method, with a size of $d_p \times d_p$ pixels. We have therefore $1 \le u, v \le 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 $\varphi$ 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 \times 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$:

$$X = AZ + B\beta \quad (1)$$

where $A$ and $B$ are two matrices, and $\beta$ 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 $\beta$ is therefore a gaussian vector with zero mean and unit variance: $\langle \beta \rangle = 0$, and $\langle \beta \beta^T \rangle$ is equal to the identity matrix.

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

$$ \langle XZ^T \rangle = A \langle ZZ^T \rangle \qquad (2) $$

thus allowing us to compute A:

$$ A = \langle XZ^T \rangle \langle ZZ^T \rangle^{-1} \qquad (3) $$

We focus now on the matrices $\langle ZZ^T \rangle$, $\langle XZ^T \rangle$, and $\langle XX^T \rangle$, 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_0$, and which agrees with the latest turbulence measurements [10]. The phase structure function then has the following expression [11]:

$$ D_{\phi}(r) = \left(\frac{L_0}{r_0}\right)^{5/3} \times \frac{2^{1/6}\Gamma(11/6)}{\pi^{8/3}} \left[\frac{24}{5}\Gamma\left(\frac{6}{5}\right)\right]^{5/6} \times \left[\frac{\Gamma(5/6)}{2^{1/6}} - \left(\frac{2\pi r}{L_0}\right)^{5/6} K_{5/6}\left(\frac{2\pi r}{L_0}\right)\right] \qquad (4) $$

where $K_{5/6}(r)$ is the modified Bessel function of the third kind, or McDonald function, $\Gamma(x)$ is the gamma function, and $r_0$ is the Fried parameter [12]. The phase structure function $D_{\phi}(r)$ is linked to the phase covariance function $C_{\phi}(r) = \langle \phi(x)\phi(x+r) \rangle$ by the relation $D_{\phi}(r) = 2[\sigma_{\phi}^2 - C_{\phi}(r)]$, $\sigma_{\phi}^2$ being the phase variance. We have therefore:

$$ C_{\phi}(r) = \left(\frac{L_0}{r_0}\right)^{5/3} \times \frac{\Gamma(11/6)}{2^{5/6}\pi^{8/3}} \left[\frac{24}{5}\Gamma\left(\frac{6}{5}\right)\right]^{5/6} \times \left(\frac{2\pi r}{L_0}\right)^{5/6} K_{5/6}\left(\frac{2\pi r}{L_0}\right) \qquad (5) $$

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)^2 + (v_j - v_i)^2]^{1/2} \times p_s$ the physical distance between these two points, then the $\langle Z_i Z_j^T \rangle$, $\langle X_i Z_j^T \rangle$ and $\langle X_i X_j^T \rangle$ terms will be given by the value of the phase covariance function $C_{\phi}(r)$ at the distance $r_{ij}$.

Equation 1 gives also the following relation:

$$ XX^T = AZZ^T A^T + AZ\beta^T B^T + B\beta ZZ^T A^T + B\beta\beta^T B^T \qquad (6) $$

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

$$ \langle XX^T \rangle = A \langle ZZ^T \rangle A^T + BB^T \Leftrightarrow BB^T = \langle XX^T \rangle - A \langle ZZ^T \rangle A^T \qquad (7) $$

and as $\langle ZZ^T \rangle$ is a symmetric matrix, its inverse is also symmetric. Then, we can write

$$ A^T = \langle ZZ^T \rangle^{-1} \langle ZX^T \rangle \qquad (8) $$

and we obtain:

$$ BB^T = \langle XX^T \rangle - A \langle ZZ^T \rangle \langle ZZ^T \rangle^{-1} \langle ZX^T \rangle = \langle XX^T \rangle - A \langle ZX^T \rangle \qquad (9) $$

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 $\lambda_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 $\lambda_i$. We can then write

$$BB^T = UWU^T = ULL^T U^T = (UL)(UL^T) \qquad (10)$$

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 \times N_{col}$ elements, this implies that $A$ is a $d_p \times (d_p \times 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 constraints on the value of $N_{col}$, we performed several simulations, for different values of $D$, $L_0$, 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 \le i \le 66$), in order to compute the variance of their coefficients $a_i$. The latter can be derived from the diagonal of the covariance matrix $\langle a_i a_j \rangle$, 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_0$ of 16 meters ($L_0/D = 2$). The ratio of $(D/r_0)$ 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 \times 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 \times 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_0$, we repeated the simulation, this time assuming an outer scale $L_0$ of 64 meters ($L_0/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$.
<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /><figcaption>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_0 = 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.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>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_0 = 64$ meters ($L_0/D = 8$).</figcaption></figure>
This study shows therefore that our method is insensitive to the ($D/L_0$) 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 \times 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_0$. 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 \times 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_0 = 16$ meters ($L_0/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_0 = 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_0$, 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 \times 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_0$. Equation 3 gives us the definition of the **A** matrix: $\mathbf{A} = \langle \mathbf{X}\mathbf{Z}^T \rangle \langle \mathbf{Z}\mathbf{Z}^T \rangle^{-1}$. As discussed in section 2, the generic term of each matrix $\langle \mathbf{X}\mathbf{Z}^T \rangle$ or $\langle \mathbf{Z}\mathbf{Z}^T \rangle^{-1}$ is given by the phase covariance function
<figure><img src="image_9.png" /><figcaption>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_0 = 16$ meters ($L_0/D = 4$)</figcaption></figure>
<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /><figcaption>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_0 = 64$ meters ($L_0/D = 4$)</figcaption></figure>
$C_φ(r)$, whose expression is:

$$C_{\phi}(r) = \left(\frac{L_0}{r_0}\right)^{5/3} \times \frac{\Gamma(11/6)}{2^{5/6}\pi^{8/3}} \left[\frac{24}{5}\Gamma\left(\frac{6}{5}\right)\right]^{5/6} \times \left(\frac{2\pi r}{L_0}\right)^{5/6} K_{5/6}\left(\frac{2\pi r}{L_0}\right) \quad (11)$$

Therefore, when the matrix product between $\langle XZ^T \rangle$ and $\langle ZZ^T \rangle^{-1}$ is performed, the terms $(L_0/r_0)^{5/3}$ cancel. This means that if we assume a constant outer scale $L_0$, A is insensitive to $r_0$.

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

$$BB^T = \langle XX^T \rangle - A \langle ZZ^T \rangle \langle ZZ^T \rangle^{-1} \langle ZX^T \rangle = \langle XX^T \rangle - A \langle ZX^T \rangle \quad (12)$$

Therefore, the matrices $\langle XX^T \rangle$ and $A \langle ZX^T \rangle$ make use of the phase covariance function $C_φ(r)$, and thus the $(L_0/r_0)^{5/3}$ term which can be introduced as a factor in the above equation. Therefore, if we compute the $BB^T$ matrix for $(L_0/r_0) = 1$, multiply the resulting B matrix by the appropriate factor $(L_0/r_0)^{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_0$. If we repeat this operation and change the $r_0$ at each iteration, then this allows us to create a phase screen with a different $r_0$ 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_0$ equal to 4 meters ($L_0/D = 10$), and a Fried parameter $r_0 = 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_0$ value of 15 cm. Next we repeated the operation, but now we assumed a sinusoidal $r_0$ 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_0$. Indeed, for a given ($D/L_0$), the Zernike polynomial variance is only proportional to $(D/r_0)^{5/3}$ [14]. Therefore, we first computed the theoretical variance of the Zernike for ($D/r_0$) = 1. Then we performed a least-square fit of the Zernike coefficient variance to the theoretical one so as to compute the $(D/r_0)^{5/3}$ scaling factor, and hence $r_0$ at each iteration. The evolution of $r_0$ 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_0$ 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_0$ error shows that it can be fitted with a Gaussian distribution, with an average error of 0.34 cm and a standard deviation $\sigma = 0.81$ cm.

We argue that this error is mainly due to the finite number of samples used to compute $r_0$ from the Zernike polynomials variance. Indeed, the bottom panels on the Fig. 5 show the result of the measure of $r_0$ 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_0 = 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_0$ is not static: we find that the $r_0$ error has an average value equal to 0.40 cm, with a standard deviation $\sigma = 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
<figure><img src="image_14.png" /><figcaption>Fig. 5. Comparison (left), error (middle), and histogram of the error (right) between the theoretical and observed Fried parameter $r_0$. The top panels correspond to a phase screen where $r_0$ 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_0 = 15$ cm. The black line on the left figure shows the expected $r_0$, whereas the red line shows the $r_0$ value deduced from Zernike polynomials variance, both as a function of the iteration number. The dashed blue line on the $r_0$ error histogram shows the result of a fit of the error by a Gaussian distribution</figcaption></figure>

of a phase screen with a spatially non uniform Fried parameter $r_0$ is straightforward to implement. The simulation of a non-static outer scale $L_0$ 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_\phi(r)$ does not scale linearly with the outer scale.

**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_0$, 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.

## Acknowledgements

The authors acknowledge the UK PPARC for financial support. François Assémat is grateful to Nazim Bharmal for helpful discussions.