Issue 
A&A
Volume 616, August 2018



Article Number  A102  
Number of page(s)  12  
Section  Numerical methods and codes  
DOI  https://doi.org/10.1051/00046361/201832579  
Published online  27 August 2018 
Numerical estimation of wavefront error breakdown in adaptive optics
LESIA, Observatoire de Paris, Université PSL, CNRS, Sorbonne Université, Univ. Paris Diderot, Sorbonne Paris Cité,
5 place Jules Janssen,
92195
Meudon, France
email: florian.ferreira@obspm.fr
Received:
3
January
2018
Accepted:
24
April
2018
Aims. Adaptive optics (AO) system performance is improved using postprocessing 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.
Methods. We have developed in the computing platform for adaptive optics systems (COMPASS) code, which is an endtoend 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 endtoend 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.
Results. On a SPHERElike system, the comparison between a PSF computed from the error breakdown with a PSF obtained from classical endtoend simulation shows that the statistics convergence limits converge very well, with a subpercent difference in terms of Strehl ratio and ensquared energy at 5λ/D separation. A correlation analysis shows significant correlations between some contributors, especially WF measurement deviation error and bandwidth error due to centroid gain, and the wellknown 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 endtoend simulation, even if some approximations exist.
Key words: instrumentation: adaptive optics / methods: numerical
© ESO 2018
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1 Introduction
Optical aberrations that are due to turbulence have a huge impact on the image resolution of groundbased large telescopes at visible and infrared wavelengths. Without any compensation, the resolution in the visible is equivalent to a diffractionlimited 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. Postprocessing 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 endtoend 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).
Endtoend AO simulation tools use a Monte Carlo approach to provide an AOcorrected 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 (Neichel et al. 2008; Jolissaint 2010; Gendron et al. 2014). 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 endtoend AO simulation (Ferreira et al. 2016). It is developed in the endtoend 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 with no other assumption than those that are made by the endtoend 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 endtoend 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 endtoend 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 Sect. 2, we discuss the error correlations that could be observed in Sect. 3. In Sect. 4 we propose a pseudoanalytical model for bandwidth and anisoplanatism errors. We also discuss the results obtained with this model and its limitations in this section. The computing performanceis summarized in Sect. 5.
2 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 tiptilt mirror both conjugated to ground, also called single conjugated adaptive optics (SCAO) system.
2.1 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 Φ _{D M} (x): (1)
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 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 oneframe delay, and using Eqs. (1)–(4), the measurement vector w_{k} of the sensor at the iteration k of the AO loop can be written as (7)
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): (8)
with 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: (9)
where R is the command matrix, which classically computed as the pseudoinverse 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 oneframe delay, we obtain a definition of the error without the fitting term: (10)
Finally, we can derive the error at the iteration k from Eqs. (7), (9), and (10): (11)
2.2 Error breakdown contributors
Writing Eq. (11) in this way allows us to reveal five contributors to ϵ_{k}: (12)
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.
2.2.2 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 nonnull measurement, reconstructed and compensated for by the AO system, introducing an aliased phase component.
2.2.4 Noise measurement error η_{k}
The noise onthe WFS measurements is the term η_{k}, which results from both readout noise and photon noise on the WFS image. Only the lowfrequency 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 whenthe 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 residualphase Φ_{⊥} 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.
2.3 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 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 tiptilt modes;

two of these modes are pure tipandtilt modes.
See Appendix A for details on the computation of this modal basis.
2.4 Validation against endtoend 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 (18)
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 socalled 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 builtin way is naturally produced by COMPASS itself as a simulation output. Thisis 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 halfwidth of . The importance of the error contributors and their possible correlations are then studied in detail in Sect. 3.
2.4.1 Simulation parameters
We consider a telescope with a diameter of 8 m without central obstruction. The parameters used in all the simulations are listed in Table 1. It corresponds to a SPHERElike system (Beuzit et al. 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 s of observation. The sky coordinates and directions are defined in a reference frame where the center is the science target and the Xaxis is oriented toward the WFS guide star.
Simulation parameters.
12layer turbulent profile.
2.4.2 Results
We can directly deduce from the ROKET outputs the conventionnal error breakdown as shown in Table 3. Crossterms are not included in this table and are detailed in Sect. 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 (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 COMPASS. 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 froman averaged phase structure function (Véran et al. 1997) and thus appears smooth and symmetric.
Error breakdown returned by ROKET for the test case.
3 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.
3.1 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 (Jolissaint 2010; Vidal et al. 2014; Martin et al. 2016). A statistical independence between the error contributors is usually assumed, so that Eq. (18) simplifies to (19)
which is, in many cases, a meaningful hypothesis.
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), 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 Sect. 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 nonnull 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 loworder 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.
Fig. 1 PSFs in log scale. Axes are expressed in units of . Top panel: PSF_{C} simulated by COMPASS. Middle panel: PSF_{R} reconstructed from ROKET buffers. Bottom panel: ∥PSF_{C} − PSF_{R}∥. 

Open with DEXTER 
Fig. 2 Cuts of PSF_{C} in red and PSF_{R} in blue. The green curve is ∥PSF_{C} − PSF_{R}∥. Top panel: cut along Xaxis. Bottom panel: cuts along Y axis. 

Open with DEXTER 
Fig. 3 Correlation matrix between the error contributors. The color scale has been changed to highlight small crosscorrelations. The diagonal is still equal to 1. 

Open with DEXTER 
3.2 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 Eq. (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, whichwas 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 remains 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 accountis only a modification of the wavefront measurement models, leading to a new error breakdown.
Fig. 4 Correlation matrix between the error contributors by taking into account the centroid gain. The color scale has been changed to highlight the vanishing of the deviation and bandwidth crosscorrelation. The diagonal is still equal to 1. 

Open with DEXTER 
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).
3.3 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 frozenflow assumption (Taylor 1938), temporal fluctuations of the turbulence can be expressed as spatial fluctuations (Greenwood 1977). Thus, we cannot neglect this correlation, especially in layerbylayer 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 Sect. 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^{−1}, 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 . 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 withthe offaxis 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.
4 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.
4.1 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 delay τ′: (20)
We find the value of τ′ assuming that we can approximate the rejection transfer function of the AO loop for oneframe delay (Demerle et al. 1994) by the transfer function of substraction with pure delay . 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 (21)
The validity of this assumption is discussed 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: (22)
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 multilayers case will come easily under the frozen flow assumption. This will allow usto 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. Equations (22) and (20) both can be rewritten in the target direction:
We emphasize that in both definitions, we consider only the phase Φ _{∥} that belongs to the DM modal space.
Mean and maximum relative error between PSF_{I} compared to PSF_{C} over the 60 simulation runs, in terms of SR and EE in for five values of the wind direction.
4.2 Covariance matrix model
From the previous definitions, we wish to derive an expression of the covariance matrix C_{b b} 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 numbers i and j, respectively. Using the identity (26)
We introduce , the structure function of the turbulent phase restricted to the spatial frequencies of the DM space. We can write it as (28)
Then, considering the separation vector x_{ij} = x_{j} −x_{i} between actuators i and j, Eq. (27) can be written as (29)
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 . We define f_{c}, and the Nyquist spatial frequency of the DM equals to where d is the actuatorpitch. Using a Fourier approach, this structure function limited to the spatial frequencies lower than f_{c} can be written as (32)
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: (33)
The covariance matrix of a complex, multilayer profile can be obtained by summing all the singlelayer matrices together.
Finally, we note with C_{ee} the covariance matrix expressed in the actuator space and filtered from piston and tiptilt modes. The removed tiptilt component is projected on the tiptilt 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).
Fig. 5 Covariance matrices of the global contribution of the anisoplanatism and bandwidth errors. Both matrix diagonals have been nullified to emphasized side structures. Top panel: covariance matrix computed from ROKET. Bottom panel: covariance matrix model C_{ee.} 

Open with DEXTER 
4.3 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, which 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_{e e}, 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. Figure 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 , with a maximum intensity at 0.78 compared to0.77, and an EE at that equals 87.5% compared to 86.8%. Then, the PSF estimation at distances farther than is less accurate.
Fig. 6 Top panel: PSF obtained from ROKET buffers of anisoplanatism and bandwidth errors. Middle panel: PSF obtained from C_{ϵϵ}. Bottom panel: absolute difference between the two PSFs. The actuator pattern is a square array. All figures are in log scale, and axes are expressed in units of . 

Open with DEXTER 
4.4 Model limitations
In this section, we apply the model we developed to the 60 simulations used in Sect. 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 Sect. 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 Sect. 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 Sect. 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 cutoff frequency of the turbulence is given by (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 highfrequency 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 andthe computing performance obtained with this model is promising.
Fig. 7 Cuts of the PSFs in log scale. Top panel: along the Xaxis. Bottom panel: along the Y axis. Blue curves: PSF reconstructed from ROKET buffers. Red curves: PSF reconstructed from the model. Green curves: Absolute difference between the two PSFs. 

Open with DEXTER 
Fig. 8 SR of PSF reconstructed from the model vs. SR of PSF reconstructed from ROKET buffers. Top left panel: loop gain of 0.1. Top right panel: loop gain of 0.2. Bottom left panel: loop gain of 0.3. Bottom right panel: loop gain of 0.4. Blue points: wind speed of 10 m s^{−1}. Green: 15 m s^{−1}. Red: 20 m s^{−1}. The red line is y = x. 

Open with DEXTER 
Fig. 9 Square modulus of the rejection transfer functionand the delay transfer function. Solid line: rejection transfer function, dashed line: delay transfer function, black dot lines: turbulence cutoff frequency depending on the wind speed. 

Open with DEXTER 
Fig. 10 Diagonal of the error covariance matrix. Blue curve: covariance computed from ROKET. Red curve: covariance computed from the model. 

Open with DEXTER 
4.5 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 firstorder 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.
5 Computing performance
The simulations presented in this paper have been run on a Nvidia DGX1 server, with Dual 20core Intel Xeon E52698 v4 @ 2.2 Ghz and 8 Tesla V100SXM2 GPUs. The GPU acceleration enables COMPASS to run the 12layer simulation, without the error estimation feature, at approximately 330 frames per second. Performing the error breakdown estimation, the computation speed decreases by a factor 5 at 60 frames per second.
At the ELT scale, COMPASS is able to run a 35layer SCAO simulation at 73 frames per second, and the error estimation is performed at 15 frames per second.
We have also developed a multiGPU module to perform PSF reconstruction using the algorithm proposed by Gendron et al. (2006). This implementation allows reconstructing a 4096 by 4096 ELT PSF on 5000 modes in 35 s on two GPUs.
The covariance matrix model described in Sect. 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 5000 by 5000 over 35 layers in half a second.
6 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 longexposure 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 offaxis 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 loworder 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 postfacto 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 pro file identification tool from SCAO telemetry data. Another step will be to upscale ROKET and the covariance model for MCAO or MOAO systems.
Acknowledgements
This work is sponsored through a grant from project #671662, a.k.a. Green Flash, funded by European Commissionunder program H2020EU.1.2.2 coordinated in H2020FETHPC2014. We also wish to thank the anonymous referee for their meaningful comments on this article.
Appendix A Modal basis computation
As specified in Sect. 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 ( with theintegral defined over the pupil area);

they are orthogonal, in the sense that the scalar product over the pupil area is null ( 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 tiptilt and are associated with the tiptilt mirror, while the tiptilt that can be produced by the DM was suppressed.
To ensure that our modal basis 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 (A.1)
Then, the dimensions of this basis are N_{Φ} × N_{actus}, where N_{Φ} is the number of points in the pupil area and N_{actus} is the number of DM actuators. The geometric covariance matrix Δ (Gaffard & Boyer 1987) is then defined as (A.2)
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 tiptilt. 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, tiptilt modes, and store them in a matrix τ: (A.3)
Our is to be orthogonal to the piston mode, and this includes pure tiptilt modes. Then, we have to generate a set of generators G from IF that cannot produce those modes (Gendron 1995): (A.4)
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, tiptilt modes, and normalization, (A.6)
we obtain an orthonormal basis B such that B^{t} ⋅ Δ ⋅ B = I.
Finally, we have to add the pure tiptilt modes in the basis B to obtain the modal basis : (A.7)
where μ are the eigenvalues of the geometric covariance matrix of the tiptilt 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 aneffect on the equations described in Ferreira et al. (2016) and Sect. 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 (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 (B.2)
Finally, the residual error ϵ_{k} becomes (B.3)
Hence, the equations for the estimation of the error breakdown contributors described in Sect. 2.2 become (B.4)
Appendix C Calculation of the lowfrequency structure function
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 of the DM.
Starting from the Kolmogorov spectrum, (C.1)
the structure function that we search for can be written as
With the variable change u = 2πkr, we obtain
However, this expression does not take into account the outer scale of the turbulence L_{0}. To include this, we note that (C.7)
From this result and Eq. (C.6), we derive an expression of the structure function limited to the frequencies higher than f_{c}: (C.8)
Moreover, we know that the effect of the outer scale on the spectrum is a saturation effect on low frequencies. Hence, the expression of is not affected since . Then, we can obtain the final as the difference between the complete structure function D_{ϕ}(r, L_{0}) and : (C.9)
References
 Beuzit, J.L., Feldt, M., Dohlen, K., et al. 2008, Groundbased and Airborne Instrumentation for Astronomy II, Proc. SPIE, 7014, 701418 [CrossRef] [Google Scholar]
 Conan, R. 2000, PhD Thesis, Sciences de l’univers Nice [Google Scholar]
 Conan, J.M., Rousset, G., & Madec, P.Y. 1995, J. Opt. Soc. Am. A, 12, 1559 [NASA ADS] [CrossRef] [Google Scholar]
 Correia, C., Véran, J.P., Ellerbroek, B., Gilles, L., & Wang, L. 2011, Second International Conference on Adaptive Optics for Extremely Large Telescopes, 70 [Google Scholar]
 Demerle, M., Madec, P. Y., & Rousset, G. 1994, in NATO Advanced Science Institutes (ASI) Series C, eds. D. M., Alloin, & J. M., Mariotti, 423, 73 [Google Scholar]
 Ferreira, F., Gendron, E., Rousset, G., & Gratadour, D. 2016, in Adaptive Optics Systems V, Proc. SPIE, 9909, 990979 [CrossRef] [Google Scholar]
 Fried, D. L. 1982, J. Opt. Soc. Am., 72, 52 [NASA ADS] [CrossRef] [Google Scholar]
 Gaffard, J. P.,& Boyer, C. 1987, Appl. Opt., 26, 3772 [NASA ADS] [CrossRef] [Google Scholar]
 Gendron, E. 1995, Theses, Université Denis Diderot (Paris 7) [Google Scholar]
 Gendron, E., Clénet, Y., Fusco, T., & Rousset, G. 2006, A&A, 457, 359 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gendron, É., Charara, A., Abdelfattah, A., et al. 2014, in Adaptive Optics Systems IV, Proc. SPIE, 9148, 91486L [Google Scholar]
 Gilles, L., Correia, C., Véran, J.P., Wang, L., & Ellerbroek, B. 2012, Appl. Opt., 51, 7443 [NASA ADS] [CrossRef] [Google Scholar]
 Gratadour, D., Puech, M., Vérinaud, C., et al. 2014, in Adaptive Optics Systems IV, Proc. SPIE, 9148, 91486O [CrossRef] [Google Scholar]
 Greenwood, D. P. 1977, J. Opt. Soc. Am., 67, 390 [NASA ADS] [CrossRef] [Google Scholar]
 Harder, S., & Chelli, A. 2000, A&AS, 142, 119 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hudgin, R. 1977, J. Opt. Soc. Am., 67, 393 [NASA ADS] [CrossRef] [Google Scholar]
 Jolissaint, L. 2010, J. Eur. Opt. Soc., 5, [arXiv:1009.1581] [CrossRef] [Google Scholar]
 Jolissaint, L., Veran, J.P., & Marino, J. 2004, in Advancements in Adaptive Optics, eds. D., Bonaccini Calia, B. L., Ellerbroek, & R., Ragazzoni, Proc. SPIE, 5490, 151 [NASA ADS] [CrossRef] [Google Scholar]
 Jolissaint, L., Véran, J.P., & Conan, R. 2006, J. Opt. Soc. Am. A, 23, 382 [NASA ADS] [CrossRef] [Google Scholar]
 Juvenal, R., Kulcsar, C., Raynaud, H.F., Conan, J.M., & Sivo, G. 2015, in Adaptive Optics for Extremely Large Telescopes IV (AO4ELT4), E63 [Google Scholar]
 Martin, O. A., Correia, C. M., Gendron, E., et al. 2016, J. Astron. Telesc. Instrum. Syst., 2, 048001 [NASA ADS] [CrossRef] [Google Scholar]
 Martin, O. A., Gendron, É., Rousset, G., et al. 2017, A&A, 598, A37 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mugnier, L. M., Fusco, T., & Conan, J.M. 2004, J. Opt. Soc. Am. A, 21, 1841 [NASA ADS] [CrossRef] [Google Scholar]
 Neichel, B., Fusco, T., & Conan, J.M. 2008, J. Opt. Soc. Am. A, 26, 219 [NASA ADS] [CrossRef] [Google Scholar]
 Rigaut, F. J., Veran, J.P., & Lai, O. 1998, in Adaptive Optical System Technologies, eds. D., Bonaccini & R. K., Tyson, Proc. SPIE, 3353, 1038 [NASA ADS] [CrossRef] [Google Scholar]
 Rousset, G., Fontanella, J. C., Kern, P., Gigan, P., & Rigaut, F. 1990, A&A, 230, L29 [NASA ADS] [Google Scholar]
 Taylor, G. I. 1938, Math. Phys. Sci. Proc. Roy. Soc. London Ser. A, 164, 476 [NASA ADS] [CrossRef] [Google Scholar]
 Véran, J.P., & Herriot, G. 2000, J. Opt. Soc. Am. A, 17, 1430 [NASA ADS] [CrossRef] [Google Scholar]
 Véran, J.P., Rigaut, F., Maitre, H., & Rouan, D. 1997, J. Opt. Soc. Am. A, 14, 3057 [NASA ADS] [CrossRef] [Google Scholar]
 Vidal, F., Gendron, E., & Rousset, G. 2010, J. Opt. Soc. Am. A, 27, A253 [NASA ADS] [CrossRef] [Google Scholar]
 Vidal, F., Gendron, É., Rousset, G., et al. 2014, A&A, 569, A16 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
All Tables
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).
Mean and maximum relative error between PSF_{I} compared to PSF_{C} over the 60 simulation runs, in terms of SR and EE in for five values of the wind direction.
All Figures
Fig. 1 PSFs in log scale. Axes are expressed in units of . Top panel: PSF_{C} simulated by COMPASS. Middle panel: PSF_{R} reconstructed from ROKET buffers. Bottom panel: ∥PSF_{C} − PSF_{R}∥. 

Open with DEXTER  
In the text 
Fig. 2 Cuts of PSF_{C} in red and PSF_{R} in blue. The green curve is ∥PSF_{C} − PSF_{R}∥. Top panel: cut along Xaxis. Bottom panel: cuts along Y axis. 

Open with DEXTER  
In the text 
Fig. 3 Correlation matrix between the error contributors. The color scale has been changed to highlight small crosscorrelations. The diagonal is still equal to 1. 

Open with DEXTER  
In the text 
Fig. 4 Correlation matrix between the error contributors by taking into account the centroid gain. The color scale has been changed to highlight the vanishing of the deviation and bandwidth crosscorrelation. The diagonal is still equal to 1. 

Open with DEXTER  
In the text 
Fig. 5 Covariance matrices of the global contribution of the anisoplanatism and bandwidth errors. Both matrix diagonals have been nullified to emphasized side structures. Top panel: covariance matrix computed from ROKET. Bottom panel: covariance matrix model C_{ee.} 

Open with DEXTER  
In the text 
Fig. 6 Top panel: PSF obtained from ROKET buffers of anisoplanatism and bandwidth errors. Middle panel: PSF obtained from C_{ϵϵ}. Bottom panel: absolute difference between the two PSFs. The actuator pattern is a square array. All figures are in log scale, and axes are expressed in units of . 

Open with DEXTER  
In the text 
Fig. 7 Cuts of the PSFs in log scale. Top panel: along the Xaxis. Bottom panel: along the Y axis. Blue curves: PSF reconstructed from ROKET buffers. Red curves: PSF reconstructed from the model. Green curves: Absolute difference between the two PSFs. 

Open with DEXTER  
In the text 
Fig. 8 SR of PSF reconstructed from the model vs. SR of PSF reconstructed from ROKET buffers. Top left panel: loop gain of 0.1. Top right panel: loop gain of 0.2. Bottom left panel: loop gain of 0.3. Bottom right panel: loop gain of 0.4. Blue points: wind speed of 10 m s^{−1}. Green: 15 m s^{−1}. Red: 20 m s^{−1}. The red line is y = x. 

Open with DEXTER  
In the text 
Fig. 9 Square modulus of the rejection transfer functionand the delay transfer function. Solid line: rejection transfer function, dashed line: delay transfer function, black dot lines: turbulence cutoff frequency depending on the wind speed. 

Open with DEXTER  
In the text 
Fig. 10 Diagonal of the error covariance matrix. Blue curve: covariance computed from ROKET. Red curve: covariance computed from the model. 

Open with DEXTER  
In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.