Numerical estimation of wavefront error breakdown in adaptive optics

Adaptive optics (AO) system performance is improved using post-processing techniques, such as point spread function (PSF) deconvolution. The PSF estimation involves characterization of the different wavefront (WF) error sources in the AO system. We propose a numerical error breakdown estimation tool that allows studying AO error source behavior such as their correlations. We also propose a new analytical model for anisoplanatism and bandwidth errors that were validated with the error breakdown estimation tool. This model is the first step for a complete AO residual error model that is expressed in deformable mirror space, leading to practical usage such as PSF reconstruction or turbulent parameters identification. We have developed in the computing platform for adaptive optics systems (COMPASS) code, which is an end-to-end simulation code using graphics processing units (GPU) acceleration, an estimation tool that provides a comprehensive error breakdown by the outputs of a single simulation run. We derive the various contributors from the end-to-end simulator at each iteration step: this method provides temporal buffers of each contributor. Then, we use this tool to validate a new model of anisoplanatism and bandwidth errors including their correlation. This model is based on a statistical approach that computes the error covariance matrices using structure functions. A correlation analysis shows significant correlations between some contributors, especially WF measurement deviation error and bandwidth error due to centroid gain, and the well-known correlation between bandwidth and anisoplanatism errors is also retrieved. The model we propose for the two latter errors shows an SR and EE difference of about one percent compared to the end-to-end simulation, even if some approximations exist.


Introduction
Optical aberrations that are due to turbulence have a huge impact on the image resolution of ground-based large telescopes at visible and infrared wavelengths. Without any compensation, the resolution in the visible is equivalent to a diffraction-limited image of a telescope with a few tens of centimeters in diameter. Adaptive optics (AO) systems have been developed for several years in astronomy (Rousset et al. 1990) to compensate for these aberrations: a wavefront sensor (WFS) measures the wavefront deformation, and a deformable mirror (DM) is controlled in real time to flatten the wavefront for the scientific image. However, the compensation by AO systems is not perfect: there is still a residual wavefront error that has variable impact on the image quality, depending on the system and the observing conditions. Post-processing techniques, such as point spread function (PSF) deconvolution algorithms, have been developed to improve the contrast on the final image (Mugnier et al. 2004). This approach requires estimating the PSF over the scientific field, which involves characterizing the error contributors of the AO system (Véran et al. 1997;Harder & Chelli 2000;Jolissaint et al. 2004;Correia et al. 2011;Gendron et al. 2014;Martin et al. 2016). Another approach is the use of end-to-end system simulation (Gilles et al. 2012). Estimating and distinguishing the various error contributors is a problem because of the propagation and filtering process of the errors in the AO loop (Vidal et al. 2014;Juvenal et al. 2015;Martin et al. 2017).
End-to-end AO simulation tools use a Monte Carlo approach to provide an AO-corrected PSF. It requires several thousand iterations to achieve sufficient convergence on the PSF computation. At the ELT (extrememly large telescope) scale, it is particularly demanding in terms of computing power and data flow. Analytical approaches have been developed to provide analytical expression of the PSF (Jolissaint 2010;Gendron et al. 2014;Neichel et al. 2008). If these models do not have convergence problems, they usually require several simplificating assumptions such as stationarity and statistical independence.
In the first part of this paper, we present the development of a new estimation tool, called ROKET (error breakdown estimation tool), which provides a comprehensive error breakdown as an output of a single run of an end-to-end AO simulation (Ferreira et al. 2016). It is developed in the end-to-end graphics processing units (GPU)-based simulation tool COMPASS (computing platform for adaptive optics systems) (Gratadour et al. 2014). The novelty of ROKET is that all the error terms are determined frame by frame at each single iteration step of the simulation, instead of only assessing their statistical behavior. This is done by using known internal quantities that are computed by the simulation to estimate these errors on the fly without interrupting the main AO loop. Thus, this tool provides an error breakdown Article number, page 1 of 12 arXiv:1804.11071v1 [astro-ph.IM] 30 Apr 2018 A&A proofs: manuscript no. main_aa_corr with no other assumption than those that are made by the endto-end models, which are usually weaker than the assumptions needed by analytical models. The goal of this tool is to quantify a detailed wavefront error breakdown in an AO system design by numerical simulations. In particular, statistical correlations between contributors can be evaluated using ROKET. The code can be use to study the behavior of error breakdown contributors and to validate or disprove strong assumptions made by an analytical model.
In the second part of this paper, we propose an analytical model for estimating the anisoplanatism and bandwidth error terms, including their correlation. For instance, this model can be used in an a posteriori PSF reconstruction algorithm, complementary to a conventional analytical model as proposed by Véran et al. (1997) or Jolissaint (2010). This analytical model is dependent on a reduced number of parameters that might be identified on AO telemetry data and by turbulence profiling tools on site. The novelty of this model is that it does not rely on a Fourier analysis, but on the error covariance matrices directly expressed in the DM space. In addition, this model is a good candidate for parallel implementation, especially using GPU acceleration, leading to efficient computing that is suitable for ELT scale.
We have realized that some models that assume independence of classical AO error terms (bandwidth, anisoplanatism, etc.) sometimes fail to reproduce the PSF obtained with an endto-end model. Therefore, it is quite understandable that mismatches of a reconstructed PSF on a real system exist. If progress has to be made in this field, it has to first address the mismatches between reconstructed PSF on simulation data. Our study intends to identify at the simplest and lowest simulation level how the errors are mixed up and what a PSF reconstruction model should take into account to reproduce end-to-end simulated data with a reliable level of accuracy. Our work in this paper intends to pave the way for a future work on PSF reconstruction that is based on telemetry data of real systems, but it remains preliminary to this future work.
After we describe the error breakdown estimation tool in Section 2, we discuss the error correlations that could be observed in Section 3. In Section 4 we propose a pseudo-analytical model for bandwidth and anisoplanatism errors. We also discuss the results obtained with this model and its limitations in this section. The computing performance is summarized in Section 5.

Error breakdown estimation tool
In this section, we present our error breakdown estimation tool. After recalling the equations that have been derived in Ferreira et al. (2016), we validate its outputs against a classical run of COMPASS. In order to make it as simple as possible in a first step, we consider an AO system with a single wavefront sensor and a single deformable mirror coupled to a tip-tilt mirror both conjugated to ground, also called single conjugated adaptive optics (SCAO) system.

AO loop error
A wavefront sensor requires a bright guide star (GS): the observed scientific object is sometimes too faint to allow the wavefront measurement. In that case, a guide star is chosen that lies as close as possible to the scientific object in order to limit the anisoplanatism error (Fried 1982). The residual phase Φ (x, θ) seen by the WFS is the turbulent phase in the GS direction Φ(x, θ), corrected by the DM, which adds a phase Φ DM (x): with θ the GS direction. Hence, we note with Φ(x, 0) the turbulent phase in the scientific direction. These quantities can be expressed on a basis B of the DM: where N is the number of commanded modes and Φ ⊥ (x, θ) is the component of Φ(x, θ) that is orthogonal to the DM space. We note with , a, and v the vectors composed of the coefficients i,θ , a i,θ, and v i , respectively. For clarity, we simplify the notation for the scientific direction: Then, assuming here a linear wavefront sensor, a system with a one-frame delay, and using Eqs. (1) to (4), the measurement vector w k of the sensor at the iteration k of the AO loop can be written as where D is he interaction matrix, and n k is the noise contribution on the measurement vector. u k covers anything that is not noise and that deviates from the theoretical model of ideal wavefront slope measurements, such as centroid gain and truncation effect. We call this term the wavefront measurement deviation. r k is the aliasing contribution in the GS direction (Véran et al. 1997;Rigaut et al. 1998), with M the linear operator that describes the ideal wavefront slope measurement from any input phase.
As the system works in closed loop, we assume that the temporal command law used in real time is an integrator with a gain g: where R is the command matrix, which classically computed as the pseudo-inverse of the interaction matrix D.
Then, if we express the residual phase in the scientific direction on the DM basis using Eq. (1) considering a one-frame delay, we obtain a definition of the error without the fitting term: Finally, we can derive the error at the iteration k from Eqs. (7), (9), and (10):

Error breakdown contributors
Writing Eq.( 11) in this way allows us to reveal five contributors to k : with 2.2.1. Bandwidth error β k β k is the temporal error that is due to the delay between the time when the command is computed and the time when this command is applied on the DM. During the computation, the turbulent phase is evolving, whereas the DM shape does not. This can be interpreted as the difference between the actual turbulent phase at time k and the command that is applied at the same time. This term only involves the phase in the direction of the scientific object.

Anisoplanatism error τ k
τ k is an anisoplanatism error that is due to the wavefront difference between the scientific direction and the analysis direction. This anisoplanatism error is estimated in the DM space, and is filtered by the temporal response through the AO loop.
2.2.3. Aliasing error ρ k ρ k is the aliasing term: the high frequencies of a turbulent phase are misinterpreted by the WFS and result in a non-null measurement, reconstructed and compensated for by the AO system, introducing an aliased phase component.

Noise measurement error η k
The noise on the WFS measurements is the term η k , which results from both read-out noise and photon noise on the WFS image. Only the low-frequency part of this noise is injected in the command.
2.2.5. WF measurement deviation error µ k µ k describes anything that deviates from an ideal wavefront slope measurement. An example of it is the truncation effect, which occurs when the field of view of subapertures of a Shack-Hartmann is too small compared to the spot size. Another example is the undersampling error that is due to the extended size of the pixel of a Shack-Hartmann. It can be estimated numerically by the difference between the WFS measurements without noise and the measurements obtained by directly computing the phase gradient of the wavefront.
2.2.6. Fitting and filtered mode error Φ ⊥ The residual phase Φ ⊥ includes a term from the filtered modes and a fitting term. These two contributions are evaluated separately.
The filtered mode error is due to the filtering process used to invert the interaction matrix. It can be estimated easily if the subspace of the modes commanded through the system control matrix is orthogonal to the subspace of filtered ones.
The fitting term is derived as the residual phase that is left after the projection onto the DM basis. As this phase is orthogonal to the DM space, it is not expressed on any DM basis. In order to estimate its impact on the Strehl ratio (SR), we compute its spatial variance at each iteration.

Modal basis
In order to exploit the error breakdown estimation for a PSF reconstruction as an example, it is more convenient to express the different contributors as a phase variance spread over a modal basis. Such modal variances can be derived from the error terms obtained and the knowledge of the DM basis or influence function used to control it. To handle it, we built a modal basis B tt from the influence functions with the following properties: -The modes are orthonormal.
-In particular, the subspace of the modes commanded through the system control matrix is orthogonal (same scalar product as above) to the subspace of filtered modes. -The modes are orthogonal to piston and tip-tilt modes.
-Two of these modes are pure tip-and-tilt modes.
See Appendix A for details on the computation of this modal basis.

Validation against end-to-end simulation
As ROKET provides the temporal buffers of each contributor estimated in the DM command space, we are able to compute any covariance matrix between two error terms and the full residual error covariance matrix as From this matrix, we are able to reconstruct the AO PSF using algorithms developed by Véran et al. (1997) or Gendron et al. (2006). These algorithms provide an estimate of the optical transfer function by computing a structure function from the AO residual error covariance matrix t . In this paper, we use the algorithm proposed by Gendron et al. (2006), which uses the so-called V ii functions that can be computed on the fly. This makes it more convenient for a GPU implementation because only a limited amount of memory is available on the device.
Then, we now have two ways of producing a PSF after a simulation run. The built-in way is naturally produced by COM-PASS itself as a simulation output. This is noted PSF C and is regarded as the reference PSF. The other way is computed from ROKET data using a PSF reconstruction algorithm.
We propose to compare the two PSFs and discuss their difference. We have chosen to compare them in terms of maximum value (or SR) and ensquared energy (EE) at a half-width of 5 λ D . The importance of the error contributors and their possible correlations are then studied in detail in Sect. 3.

Simulation parameters
We consider a telescope with a diameter of 8m without central obstruction. The parameters used in all the simulations are listed in Table 1. It corresponds to a SPHERE-like system (Beuzit et al.  Table 2. 12-layer turbulent profile 2008). The parameters regarding the atmospheric conditions are detailed in Table 2. The simulations was run over 20 000 iterations, which is equivalent to 40 seconds of observation. The sky coordinates and directions are defined in a reference frame where the center is the science target and the X-axis is oriented toward the WFS guide star.

Results
We can directly deduce from the ROKET outputs the conventionnal error breakdown as shown in Table 3. Cross-terms are not included in this table and are detailed in Section 3. The error breakdwon obtained is dominated by the anisoplanatism and bandwidth terms, and it is compliant with common approximations used to estimate the contributors, such as the fitting error estimated from (Hudgin 1977) using

Contributors
(with d the inter actuator distance), leading here to 46 nm rms. The variance of aliasing, set to 33% of the fitting (Rigaut et al. 1998), gives 31 nm rms.
To validate this error breakdown, we reconstructed the PSF R by directly computing the residual error covariance matrix from the sum of the error buffers returned by ROKET just as in Eq. (18), and we compared it to the PSF C simulated by COM-PASS. Figure 1 shows these two PSFs and their absolute difference in log scale. Cuts along the X-and Y-axes are shown in Fig. 2.
We note the perfect agreement between the two PSFs, with an estimated SR of 73.6% by ROKET compared to 73.8% computed by COMPASS, and an EE of 83.2% compared to 83.3%, respectively. This demonstrates that the estimation of ROKET is accurate, as we can reliably retrieve the PSF from it. Of course, the PSF C exhibits a speckled aspect because the simulation does not converge, while the PSF R is built from an averaged phase structure function (Véran et al. 1997) and thus appears smooth and symmetric.

Correlation analysis
This section takes advantage of ROKET capacities to highlight the correlations that exist between some of the error breakdown contributors. We also study the conditions that lead to these correlations.

Correlation between error terms
Developing Eq. 18 leads to 5 main terms that are the covariance matrices of each error breakdown contributor, and to 20 crossterms that are often neglected in many studies (Martin et al. 2016;Jolissaint 2010;Vidal et al. 2014). A statistical independence between the error contributors is usually assumed, so that Eq. (18) simplifies to However, we cannot infer with certainty that this assumption is always valid. Obviously, some contributors of the error breakdown are effectively independent (such as noise with respect to turbulent terms, e.g.), but we have no guarantee that turbulencerelated terms are independent among themselves. This assumption could have an impact in the PSF estimation as it directly affects the error covariance matrix. With ROKET, we will be able to assess the validity of assuming that the error term is statistically independent, and determine what to do if this is not true.
We are able to compute all covariance terms, including crossterms, from the temporal buffers of each error contributor returned by ROKET. Figure 3 shows the correlation coefficients obtained between all error contributors for the simulated case described in Section 2.4.1. They are displayed in matrix form, with lines and rows associated with a given contributor. Thus, for two contributors x and y, the value of the cell that intersects line x and column y is the correlation coefficient between these terms. We note a strong anticorrelation of -0.43 between the WF measurement deviation and bandwidth. Other non-null correlations are a WF measurement deviation with aliasing (-0.03), anisoplanatism (0.03), noise (-0.02), and also bandwidth and anisoplanatism (0.03). The WF measurement deviation affects the measurements of the low-order modes, so that it is naturally correlated to the temporal and the aliasing error, which are both made of low orders. ROKET shows a correlation between WF measurement deviation and noise. The Shack-Hartmann image profile affects the centroiding noise n k and also the centroid gain, which is part of the deviation term. Consequently, it is impossible to argue that these terms cannot be correlated, as they both come from the centroiding process of the image.

WF measurement deviation and bandwidth correlation
The correlation between deviation and bandwidth is first due to a centroid gain that can be modeled by a multiplicative scalar factor γ on the WFS measurement. By construction of ROKET, the effect of the centroid gain is counted as part of u k -an odd behavior. On an other hand, the effect of a centroid gain is often regarded as the equivalent to an alteration of the loop gain and is expected to primarily affect the bandwidth term, and to a lesser extent the propagation of noise, aliasing, and other terms (Véran & Herriot 2000). It would make more sense to identify this centroid gain as such to establish a new error breakdown. We have modified the wavefront measurement models described in Ferreira et al. (2016) for this purpose. For the sake of clarity, the new equations are presented in Appendix B, Eq. (B.4). We highlight one point: the summation of the five equations obtained in (B.4) will always result in the same k , by definition. If the decomposition is affected by γ, the total remains unchanged because the AO system is not retuned with respect to γ. As a result, different decompositions using different gamma will lead to exactly the same PSF. The γ parameter only changes the breakdown by transferring variances and covariances from one item to another. Now, using the particular value of γ that matches the system centroid gain has some particular property: it causes the covariances to vanish. Of course, a new simulation run needs to be performed because the centroid gain γ also needs to be applied to the phase gradient WFS model. To compute the correct value of γ, ROKET performs a linear regression between the WFS measurement without noise and the phase gradient model measurement, in a first run. The value of γ is then estimated as the average of the regression coefficients over the iterations. Finally, we need to launch a new ROKET simulation run with this γ value in order to compute the error breakdown as given in Appendix B. Figure 4 shows the correlation coefficients obtained with these equations. We note that the WF measurement deviation term is significantly reduced and the correlation between WF measurement deviation and bandwidth has been zeroed. We also note the disappearance of the correlation between WF measurement deviation and noise, which was then effectively due to the existence of this centroid gain. A correlation remains with aliasing. Introducing this centroid gain γ in the equations has moved variances and covariances around, as shown in Table 4. It is worth noticing that the error breakdown returned by ROKET re-  Table 4. Variances and covariances obtained with γ = 1 (i.e., without taking the centroid gain into account) and with γ = 0.95 (value computed after the first run of ROKET). σ 2 η is thenoise variance, σ 2 µ is the WF measurement deviation variance, σ 2 ρ is the aliasing variance, σ 2 β is the bandwidth variance, σ 2 τ is the anisoplanatism variance, and σ 2 µβ is the covariance between the WF measurement deviation and bandwidth. All are expressed in 10 −3 µm 2 mains valid in any case, since the sum of the time buffers remains unchanged. Similarly, the sum of all elements of the covariance tables in Figs. 3 and 4 leads to the same total variance, regardless of the value of γ. Taking γ into account is only a modification of the wavefront measurement models, leading to a new error breakdown.

Anisoplanatism and bandwidth correlation
In most cases, bandwidth and anisoplanatism are the main contributors of the error breakdown of a SCAO system. Even if the correlation between anisoplanatism and bandwidth errors appears to be low, it is well known that it strongly depends on wind conditions. Under the frozen-flow assumption (Taylor 1938), temporal fluctuations of the turbulence can be expressed as spatial fluctuations (Greenwood 1977). Thus, we cannot neglect this correlation, especially in layer-by-layer error modeling as in (Jolissaint et al. 2006).
To illustrate this purpose and prepare further model validation, we ran 60 ROKET simulations with the same conditions as described in Section 2.4.1, but with only a single layer at 5 km altitude and various wind conditions: a windspeed between 5 and 20 m/s, and a wind direction between 0 • and 180 • . We compare then the PSF R computed by ROKET, and the PSF I computed by neglecting anisoplanatism and bandwidth correlation. The comparison is presented in terms of SR and EE at ±5 λ D . Table 5 summarizes the results obtained for each wind direction in terms of relative error on the SR and the EE. As expected, the anisoplanatism and bandwidth correlation depends on the wind direction. It is strongest when the wind is aligned with the off-axis GS direction. Neglecting it in this case leads to huge errors in the PSF estimate. Conversely, the correlation is negligible when the wind direction is orthogonal to that of the GS.

Bandwidth and anisoplanatism models
This work prepares some tools for PSF reconstruction based on the telemetry data of a real system. We propose an analytical model of the covariance matrix between bandwidth and anisoplanatism to be retrieved from these telemetry and turbulence profiling tools. Similar models already exist that aim at PSF computation Jolissaint et al. (2006); Neichel et al. (2008). In this section, we propose a different approach that does not rely on a Fourier analysis, but on an expression of covariance matrices of the anisoplanatism and bandwidth errors, including their possible correlation, directly in the DM space. Applications and limitations of this approach are discussed.

Definitions
We define the bandwidth error as the temporal error due to the delay between the time when the wavefront is measured and the time when this command is applied on the DM. Then, we write the bandwidth error b as a phase difference in the direction of the target between two moments that are separated by some de- We find the value of τ assuming that we can approximate the rejection transfer function of the AO loop for one-frame delay H(p) = 1 1+g Fs p e −τp (Demerle et al. 1994) by the transfer function of substraction with pure delay H (p) = 1 − e −τ p . The variable p is the Laplace variable, F s is the sampling frequency of the AO system, and τ = 1/F s is the sampling period. Considering the Taylor series of both expressions around p = 0 at first order, they will match for the value The validity of this assumption is discuss later. In a real AO system, we should include the filtering by the exposure time on the WFS and by the blocker of the DM command. In any case, however, our pure delay identification procedure can still be applied. We define the anisoplanatism error as the phase difference between wavefront sensor and the target directions: where θ is the angular separation between the WFS guide star and the target. Now, we assume the case of a phase in a single layer. The multi-layers case will come easily under the frozen flow assumption. This will allow us to transform angular coordinates of observation into spatial coordinates at the layer surface, knowing the turbulent layer altitude h. We also assume a frozen flow hypothesis, which allows us to turn time into space using the wind velocity vector v. Eqs. (22) and (20) can both be rewritten in the target direction: We emphasize that in both definitions, we consider only the phase Φ that belongs to the DM modal space.

Covariance matrix model
From the previous definitions, we wish to derive an expression of the covariance matrix C bb of the bandwidth error, of the covariance matrix C aa of the anisoplanatism error, and of the covariance matrix C ba between these errors. We first compute these covariance matrices on the DM actuators. Components (i, j) of the covariance matrix C ba can be written as (25) where x i and x j are the position vectors in the pupil of the DM actuator number i and j, respectively. Using the identity (27) We introduce D low φ (r), the structure function of the turbulent phase restricted to the spatial frequencies of the DM space. We can write it as Then, considering the separation vector x ij = x j − x i between actuator i and j, Eq. (27) can be written as The same approach stands for covariance matrices C aa and C bb and leads to the very similar expressions These expressions are simple and fully analytical as we are able to compute the structure function D low φ (r). We define f c , and the Nyquist spatial frequency of the DM equals to 1 2d, where d is the actuator pitch. Using a Fourier approach, this structure function limited to the spatial frequencies lower than f c can be written as The D φ (r, L 0 ) expression can be found in Conan (2000). This structure function is an approximation that considers circular symmetry and not the actual actuator pattern. The detailed calculation and computation methods are reported in Appendix C. Then, we obtain the full contribution of the anisoplanatism and bandwidth errors by adding these covariance matrices: The covariance matrix of a complex, multi-layer profile can be obtained by summing all the single-layer matrices together. Finally, we note with C ee the covariance matrix expressed in the actuator space and filtered from piston and tip-tilt modes. The removed tip-tilt component is projected on the tip-tilt mirror actuators.
This model only needs a few parameters from the atmospheric turbulence (L 0 , r 0 , turbulent profile, and wind) and from the AO system (loop frequency, actuator positions, and guide star position).

Model result
To assess the validity of our model, we generated PSFs following two methods. The first method used the error buffers from ROKET of bandwidth and anisoplanatism, wich are summed together and converted into a covariance matrix that will feed the PSF reconstruction algorithm. The second method used Eq. (33) to compute the C ee , which also feeds the PSF reconstruction algorithm. Figure 5 shows both covariance matrices expressed in the DM actuators space. We note that the model reproduces the features of the covariance matrix well. Fig. 6 shows the PSFs obtained through the two methods and their difference. Cuts of these PSFs along the X-and Y-axes are shown in Fig. 7.
We note a good accuracy in the PSF reconstruction up to 5 λ D , with a maximum intensity at 0.78 compared to 0.77, and an EE at ±5 λ D that equals 87.5 % compared to 86.8 %. Then, the PSF estimation at distances farther than 10 λ D is less accurate.

Model limitations
In this section, we apply the model we developed to the 60 simulations used in Section 3.3. The results obtained highlight some limitations of the model that we explain below. We now apply the same methods as we used in Section 4.3 to compute PSFs for the same 60 simulation cases. Figure 8 shows the SR obtained from a PSF reconstructed from C versus the one obtained from PSF reconstructed from ROKET buffers of anisoplanatism and bandwidth errors. Globally, the faster the wind and the lower the loop gain, the less efficient the model. This is explained by the approximation we made in Section 4.1 concerning the rejection transfer function. This assumption is made to obtain a simple model that could be computed fast and accurate in most cases, as shown in Section 4. However, the validity of this assumption depends on the loop gain, the loop frequency, and the wind speed. Figure 9 compares the modulus of the rejection transfer function and the modulus of the delay transfer function for various gains. The higher the gain, the closer the transfer functions, and so the more valid the assumption. Moreover, the approximation appears to be better for low frequencies than for high frequencies. As the cut-off frequency of the turbulence is given by f c ≈ 0.3 v d (Conan et al. 1995), where d is the subaperture diameter, the approximation, and so the model that we have developed, is better for low wind speed. The results obtained with Fig. 8 confirm this behavior.
It also leads to a poor estimation of the high-frequency components of the error. Figure 10 shows the diagonal of the covariance matrix C in the modal basis space, compared to the covariance matrix computed from ROKET. The model is accurate for modes from 0 to 800, then the covariance for higher modes is underestimated.
However, the results obtained with the model encourage us to keep the assumptions we made. For most cases, the accuracy is relevant and the computing performance obtained with this model is promising.

Model applications
As described above, models for anisoplanatism and bandwidth errors are available that use a power spectral density analysis of the AO residual phase Jolissaint et al. (2006). The model that we propose is based on a slightly different approach that will be extended in future works to model the other error breakdown contributors. It will provide another method for estimating the PSF of an AO system. It will allow quick first-order AO performance evaluation even at the ELT scale.
Moreover, this approach can easily take into account the DM actuator geometry such as the hexagonal pattern and the shape of the influence functions within the limit of a modal basis that is restricted to a circular frequency domain. It could also be used in a PSF reconstruction algorithm. It might also be useful in practical cases since it is based on an estimation of covariance matrices in the DM actuator space, which is a space available from AO loop telemetry. For example, a possible application could be the identification of turbulent parameters from AO loop data using a method similar to the one developed by Vidal et al. (2010) based on a fitting algorithm that minimizes the difference between the modeled covariance matrix with the matrix computed from telemetry data. Thus, retrieving the turbulence parameters from an SCAO system will be possible.

Computing performance
The simulations presented in this paper have been run on a Nvidia DGX-1 server, with Dual 20-core Intel Xeon E5-2698 v4 @ 2.2 Ghz and 8 Tesla V100-SXM2 GPUs. The GPU acceleration enables COMPASS to run the 12-layer simulation, without the error estimation feature, at approximately 330 frames per At the ELT scale, COMPASS is able to run a 35-layer SCAO simulation at 73 frames per second, and the error estimation is performed at 15 frames per second. We have also developed a multi-GPU module to perform PSF reconstruction using the algorithm proposed by Gendron et al. (2006). This implementation allows reconstructing a 4 096 by 4 096 ELT PSF on 5 000 modes in 35 seconds on two GPUs. The covariance matrix model described in Section 4 is highly parallelizable as each element of the matrix can be computed independently of the other elements. Hence, we have developed a GPU module dedicated to the computation of those matrices that lead to a computation of an ELT scale covariance matrix 5 000 by 5 000 over 35 layers in half a second.

Conclusion and perspective
We have developed a numerical tool, called ROKET, which is included in the COMPASS AO simulation code. It produces a comprehensive breakdown of the wavefront error for an AO loop. The estimated contributors are the temporal error, the anisoplanatism, the aliasing, the WF measurement deviation including centroid gain, the fitting, and the filtered mode error. The tool allows us to evaluate not only the variances, but the possible correlations between the different contributors to the error breakdown.
Our results show that the error breakdown estimation is accurate and fast. The PSF, reconstructed from a ROKET error breakdown, is very similar to the empirical long-exposure PSF computed during the simulation run by COMPASS. The computed error breakdowns show that correlations exist between WF measurement deviation, noise, aliasing, temporal error, and anisoplanatism. We focused on the bandwidth and anisoplanatism errors, which are two main contributors in off-axis WFS. Their correlation often cannot be neglected.
We have developed an analytical model of these two main contributors and their correlation. It allows us to compute the resulting error covariance matrix that can be directly used by classical PSF reconstruction algorithms. Our results show that this model is accurate, even if the accuracy reached for the highest order controlled modes is lower than the one for low-order modes.
This model only requires a few parameters of the AO system and of the atmospheric conditions to compute the covariance matrix. The novelty compared to other analytical models is that it is expressed directly on the DM actuator space, allowing then convenient implementation for post-facto PSF reconstruction, for instance. In addition, the computing efficiency can be greatly improved by a GPU implementation, which allows easily handling an ELT scale.
Future works will use ROKET as a reference tool to develop the same type of models for the other error breakdown contributors. With such models, fast and accurate PSF estimation of an AO system will be possible, even at the ELT scale. We plan to use ROKET for the validation by simulation of a turbulent profile identification tool from SCAO telemetry data. Another step will be to upscale ROKET and the covariance model for MCAO or MOAO systems. F. Ferreira et al.: Error breakdown estimation for adaptive optics systems Appendix A: Modal basis B t t computation As specified in Section 2.3 , the error breakdown needs to be estimated on a modal basis with specific properties: the modes span the full DM space the modes are normalized ( 1 S B 2 i dS = 1µm 2 , ∀i with the integral defined over the pupil area) they are orthogonal, in the sense that the scalar product over the pupil area is null ( 1 S B i B j dS = 0 for i j) the subspace of the modes commanded through the system control matrix is orthogonal (same scalar product as above) to the subspace of filtered ones one of the filtered modes is constructed as the best DM leastsquares fit to piston, so that all modes are orthogonal to the piston (∀i, B i dS = 0). This property is important for deriving phase variances that do not include any component along the piston term. -Two of these modes are pure tip-tilt and are associated with the tip-tilt mirror, while the tip-tilt that can be produced by the DM was suppressed.
To ensure that our modal basis B tt spans the full DM space, it is computed from the influence functions of the DM. Let IF the basis of influence functions of the DM, that is, any phase produced by the DM can be decomposed as As IF is a basis, ∆ is invertible and ∆ −1 . IF t projects a phase onto the IF basis. Let T p be the matrix containing the phase corresponding to a piston and pure tip-tilt. The dimensions of this matrix are therefore N Φ × 3. Using the projection matrix defined above, we can retrieve the coefficients a i on the basis IF that fits the piston, tip-tilt modes, and store them in a matrix τ: Our B tt is to be orthogonal to the piston mode, and this includes pure tip-tilt modes. Then, we have to generate a set of generators G from IF that cannot produce those modes (Gendron 1995): where I is the identity matrix. Now, the diagonalization of G provides a basis B : (A.5) where λ are the eigenvalues of G t . ∆ . G . After truncation of the three last columns of B corresponding to piston, tip-tilt modes, and normalization, we obtain an orthonormal basis B such that B t . ∆ . B = I.
Finally, we have to add the pure tip-tilt modes in the basis B to obtain the modal basis B tt , where µ are the eigenvalues of the geometric covariance matrix of the tip-tilt modes, corresponding to its diagonal as these two modes composed a basis.

Appendix B: AO loop error with centroid gain
A centroid gain γ will have an effect on the equations described in Ferreira et al. (2016) and Section 2. To take it into account in ROKET, we have to rewrite the equations with a centroid gain. It can be modeled as a gain γ on the WFS measurement, so that Eq. (7) becomes w k = γDa k,θ + γDv k−1 + γr k + n k + u k . (B.1) In this expression, u k is no longer the same vector as in Eq. (7) as it does not include the centroid gain error. Then, the command vector v k can be written as Finally, the residual error k becomes Hence, the equations for the estimation of the error breakdown contributors described in Section 2.2 become β k = (1 − γgRD)β k−1 + (a k − a k−1 ) τ k = (1 − γgRD)τ k−1 + γgRD(a k−1 − a k−1,θ ) ρ k = (1 − γgRD)ρ k−1 − γgRr k−1 η k = (1 − γgRD)η k−1 − gRn k−1 µ k = (1 − γgRD)µ k−1 − gRu k−1 .
(B.4) Appendix C: Calculation of the low-frequency structure function D low φ To compute a model for the covariance between anisoplanatism and bandwidth error as a sum of structure functions, we have to compute these functions. However, as we are searching for errors made on the DM command, we need to calculate the structure function that excludes spatial frequencies higher than the Nyquist frequency f c = 1 2d of the DM. Starting from the Kolmogorov spectrum, W(k) = 0.023 r −5/3 0 (k) −11/3 , (C.1) the structure function that we search for can be written as