© ESO 2018

1 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 (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 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 with no other assumption than those that are made by the end-to-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 end-to-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 Sect. 2, we discuss the error correlations that could be observed in Sect. 3. In Sect. 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 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 tip-tilt 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 one-frame 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 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: (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)

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.

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 non-null 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 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 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 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.

2.4 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 (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 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 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 half-width 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 SPHERE-like 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 X-axis is oriented toward the WFS guide star.

2.4.2 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 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.

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 cross-terms 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 turbulence-related 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 cross-terms, 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 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.

Fig. 1 PSFs in log scale. Axes are expressed in units of . Top panel: PSFC simulated by COMPASS. Middle panel: PSFR reconstructed from ROKET buffers. Bottom panel: ∥PSFC − PSFR∥.

Fig. 2 Cuts of PSFC in red and PSFR in blue. The green curve is ∥PSFC − PSFR∥. Top panel: cut along X-axis. Bottom panel: cuts along Y -axis.

Fig. 3 Correlation matrix between the error contributors. The color scale has been changed to highlight small cross-correlations. The diagonal is still equal to 1.

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 cross-correlation. The diagonal is still equal to 1.

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 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 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⁻¹, 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 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.

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 one-frame 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 multi-layers 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.

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)

leads to (27)

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₀) 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, 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₀, r₀, 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 Cee.

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 .

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 cut-off 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 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 andthe computing performance obtained with this model is promising.

Fig. 7 Cuts of the PSFs in log scale. Top panel: along the X-axis. 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.

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.

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 cut-off frequency depending on the wind speed.

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

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 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.

5 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 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 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 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 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 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 H2020-EU.1.2.2 coordinated in H2020-FETHPC-2014. 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 least-squares fit to piston, so that all modes are orthogonal to the piston (∀i, ∫ B_idS = 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 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 Δ⁻¹ ⋅ 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 τ: (A.3)

Our 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): (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, tip-tilt modes, and normalization, (A.6)

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 : (A.7)

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 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 low-frequency 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₀. 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₀) and : (C.9)

References

All Tables

All Figures

	Fig. 1 PSFs in log scale. Axes are expressed in units of . Top panel: PSFC simulated by COMPASS. Middle panel: PSFR reconstructed from ROKET buffers. Bottom panel: ∥PSFC − PSFR∥.
In the text

	Fig. 2 Cuts of PSFC in red and PSFR in blue. The green curve is ∥PSFC − PSFR∥. Top panel: cut along X-axis. Bottom panel: cuts along Y -axis.
In the text

	Fig. 3 Correlation matrix between the error contributors. The color scale has been changed to highlight small cross-correlations. The diagonal is still equal to 1.
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 cross-correlation. The diagonal is still equal to 1.
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 Cee.
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 .
In the text

	Fig. 7 Cuts of the PSFs in log scale. Top panel: along the X-axis. 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.
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.
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 cut-off frequency depending on the wind speed.
In the text

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