© ESO, 2017

1. Introduction

New scientific challenges will become achievable thanks to the future 30–40 m class telescopes such as the ESO European Extremely Large Telescope (E-ELT; McPherson et al. 2012). In particular, it will be possible to track properties of first galaxies at high red-shifts (Evans et al. 2008). Such distant objects are smaller than the seeing disk, and so faint that they may require hours of integration time. They must be observable by integral field spectroscopy and statistics on these objects will be provided by a large multiplex system. Such a statistical analysis requires the observation of the greatest possible number of galaxies distributed over a very large field of view (FoV). Thanks to the large light-collecting power of the E-ELT combined with an adaptive optics (AO) system having high sky coverage, such observations will be feasible.

To perform the multi-object correction, an AO concept has been proposed (Hammer et al. 2002): multi-object adaptive optics (MOAO). A MOAO system is composed of individual optical trains that split the FoV into different scientific directions. Each of these optical trains includes a single deformable mirror (DM) which corrects the turbulence over the small sub-field (of the order of an arc-second) encircling the galaxy. In addition in MOAO, objects to be observed are too faint to provide the number of photons required to measure the wave front (WF). The WF sensing is thus achieved thanks to other natural guide stars (NGS) within the field and laser guide stars (LGS). Unfortunately, the turbulence measured by WF sensors (WFS) decorrelates with angular distance because of the anisoplanatism effect (Fried 1982). The mitigation of this effect requires a tomographic approach: the DMs are driven from a combination of the off-axis WF measurements taking into account their angular correlation. Using several NGS and LGS in the FoV, MOAO estimates the wave-front in the scientific directions by a tomographic reconstruction to mitigate the anisoplanatism effect.

Eagle (Cuby et al. 2010) is a MOAO-assisted multi-object spectrograph (MOS), that has been proposed for the E-ELT for the statistical study of the formation process of distant galaxies. The concept of Eagle is to cover a 7 arcmin diameter FoV and features 20 MOAO spectroscopic channels, working in parallel, and delivering a corrected wave-front for achieving an ensquared energy (EE) greater than 30% in a resolution element of 75 mas. The AO study (Rousset et al. 2010) has demonstrated that the tomographic reconstruction must be based on a minimum of six LGS WFS and four to five NGS WFS in order to be able to reach the correction performance and the sky coverage requirement. Today, a new multi-object spectrograph, Mosaic, possibly equipped with MOAO, is entering in phase A study for the E-ELT (Hammer et al. 2014).

To demonstrate the technical feasibility of MOAO sky observations, the pathfinder Canary (Myers et al. 2008) has been developed by an international consortium. It began in 2007, with first light in 2010 at the William Herschel Telescope (WHT) on La Palma, in the Canary islands (Spain) (Gendron et al. 2011). Canary is also a platform where innovative techniques and algorithms can be tested on-sky for AO control techniques (Sivo et al. 2014; Osborn et al. 2014; Bitenc et al. 2015) and Shack-Hartmann (SH) centroiding (Basden et al. 2012, 2014) for instance.

In this paper, we present the analysis of on-sky results acquired by Canary in 2013. We first describe the MOAO control and tomographic reconstruction with LGS in Sect. 2, as used with Canary. In particular, we discuss the implementation of the mixed NGS/LGS tomographic reconstruction.

In Sect. 3, we introduce the tomographic error in LGS MOAO and describe how its evaluation is achieved for Canary. In Sect. 4, we present an analytic LGS MOAO error breakdown and describe how each term is derived. The variance of the residual phase taken in summing all terms, assumed to be independent, is compared to what we estimate from the truth sensor (TS). We show in Sect. 5 the two methods lead to the same evaluation of the residual phase, with a 99% correlation over a large set of 4500 on-sky datasets. We also present an overview of the on-sky Canary performance in single conjugate AO (SCAO), MOAO, and open-loop ground layer AO (GLAO).

Finally in Sect. 6, we focus on particular results acquired by Canary running successively in SCAO, MOAO and GLAO modes. Using tools introduced in Sect. 4, we evaluate the Canary error breakdown in these modes and analyze the improvement achieved using LGS/NGS tomography. We conclude in Sect. 7.

2. Canary design: phase B

2.1. Phase B design

The Canary pathfinder is deployed on one of the Nasmyth focii of the 4.2 m WHT. The project has been planned with several phases of increasing system complexity (Myers et al. 2008) to lead to a comprehensive demonstration of the MOAO configuration as foreseen for Eagle on the E-ELT (Rousset et al. 2010).

The Canary optical design during phase B (Morris et al. 2013) includes up to eight SH WFS of 36 valid sub-apertures each (7 × 7), three being off-axis NGS WFS, four being off-axis Rayleigh LGS WFS and one being an on-axis NGS WFS. They are installed behind an optical de-rotator (K mirror) placed at the entrance of the bench to compensate for the field rotation at the Nasmyth focus. The three off-axis NGS WFS are installed on a target acquisition system (TAS) in the telescope focal plane which allows us to point at directions anywhere within a FoV of 2.5′ diameter. The four Rayleigh LGS WFS have a selectable range gate (to tune the guide star height and depth), set to 21 km for the data presented in this paper, and are optically combined on a unique electronically time-gated CCD (Morris et al. 2013, 2014). The on-axis TS measures the WF corrected by the DM and tip-tilt (TT) mirror for calibration procedures and for performance diagnostics. A dichroic plate splits the corrected WF beam into a visible part (400–900 nm) for the TS and an infra-red part (900–2500 nm) for an infra-red imaging camera, called Camicaz, equipped with a NICMOS detector and a filter wheel delivering J, H and K-band images (Gratadour et al. 2013a).

All the WFS cameras are triggered with a master clock derived from a sub-frequency of the laser pulses, leading to synchronous WFS measurements sampled at 150 Hz for all cameras. The real-time computer (RTC) DARC (Basden et al. 2010; Basden et al. 2016) acquires synchronously the WFS data and provides the user with time series of wave-front slopes and DM and TT commands.

2.2. Wave-front sensor (WFS) slope vectors

The RTC provides for each WFS a vector S_p of n_s (see Table 1) slopes (n_s = 2 × 36) every WFS frame. The set of local WF slopes along x and y directions, measured by the pth WFS with n_s/ 2 sub-apertures, is collected into a global slope vector S_p: $\begin{array}{ccc}{{S}}_{\mathit{p}}\mathrm{=}\left(\hspace{0.17em}\begin{array}{c}\begin{array}{c}\end{array}{\mathit{S}}_{{\mathit{x}}_{\mathrm{1}}}\begin{array}{c}\end{array}\\ \begin{array}{c}\end{array}\begin{array}{c}\mathit{.}\\ \mathit{.}\\ \mathit{.}\end{array}\\ \begin{array}{c}\end{array}{\mathit{S}}_{{\mathit{x}}_{{\mathit{n}}_{\mathrm{s}}\mathit{/}\mathrm{2}}}\end{array}\right)\mathrm{-}{{S}}_{\mathit{p}}^{\mathrm{ref}}\mathit{,}& & \end{array}$(1)where ${{S}}_{\mathit{p}}^{\mathrm{ref}}$ is the vector of the reference slopes, previously calibrated in the laboratory (Vidal et al. 2014). To remove tip-tilt signal from LGS measurements, we subtract the average slopes on both the x and y axis: $\hspace{0.17em}\begin{array}{c}\begin{array}{c}\end{array}{{S}}_{\mathrm{LGSp}}\mathrm{=}\left({{I}}_{{\mathit{n}}_{\mathrm{s}}}\mathrm{-}{F}\right)\mathrm{\times }{{S}}_{\mathit{p}}\mathit{,}\end{array}$(2)where I_ns is the identity matrix of size n_s and F is defined as: ${F}\mathrm{=}\frac{\mathrm{2}}{{\mathit{n}}_{\mathrm{s}}}\left[\begin{array}{c}\\ \mathrm{[}{1}{\mathrm{]}}_{{\mathit{n}}_{\mathrm{s}}\mathit{/}\mathrm{2}}& \mathrm{[}{0}{\mathrm{]}}_{{\mathit{n}}_{\mathrm{s}}\mathit{/}\mathrm{2}}\\ \mathrm{[}{0}{\mathrm{]}}_{{\mathit{n}}_{\mathrm{s}}\mathit{/}\mathrm{2}}& \mathrm{[}{1}{\mathrm{]}}_{{\mathit{n}}_{\mathrm{s}}\mathit{/}\mathrm{2}}\end{array}\right]\mathit{,}$(3)where [1] _{ns/ 2} and [0] _{ns/ 2} are squared matrices of size n_s/ 2 filled respectively with ones and zeros (see Table 1).

We denote S_off as the array resulting from the concatenation of slope measurements in the off-axis directions. This therefore includes the slopes coming from 4 LGS and 3 NGS, and is of size (2 × 36 × 7) × n_f. Likewise, we denote S_on as the array of size n_s × n_f that concatenates the TS measurements on the on-axis source. We have: $\begin{array}{ccc}{{S}}_{\mathrm{off}}\mathrm{=}\left(\hspace{0.17em}\begin{array}{c}\begin{array}{c}\end{array}{{S}}_{\mathrm{LGS}\mathrm{1}}\mathrm{(}\mathit{t}\mathrm{=}\mathrm{1}\mathrm{)}\mathrm{·}\mathrm{·}\mathrm{·}{{S}}_{\mathrm{LGS}\mathrm{1}}\mathrm{(}\mathit{t}\mathrm{=}{\mathit{n}}_{\mathrm{f}}\mathrm{)}\begin{array}{c}\end{array}\\ \begin{array}{c}\end{array}{{S}}_{\mathrm{LGS}\mathrm{2}}\mathrm{(}\mathit{t}\mathrm{=}\mathrm{1}\mathrm{)}\mathrm{·}\mathrm{·}\mathrm{·}{{S}}_{\mathrm{LGS}\mathrm{2}}\mathrm{(}\mathit{t}\mathrm{=}{\mathit{n}}_{\mathrm{f}}\mathrm{)}\\ \begin{array}{c}\end{array}{{S}}_{\mathrm{LGS}\mathrm{3}}\mathrm{(}\mathit{t}\mathrm{=}\mathrm{1}\mathrm{)}\mathrm{·}\mathrm{·}\mathrm{·}{{S}}_{\mathrm{LGS}\mathrm{3}}\mathrm{(}\mathit{t}\mathrm{=}{\mathit{n}}_{\mathrm{f}}\mathrm{)}\\ \begin{array}{c}\end{array}{{S}}_{\mathrm{LGS}\mathrm{4}}\mathrm{(}\mathit{t}\mathrm{=}\mathrm{1}\mathrm{)}\mathrm{·}\mathrm{·}\mathrm{·}{{S}}_{\mathrm{LGS}\mathrm{4}}\mathrm{(}\mathit{t}\mathrm{=}{\mathit{n}}_{\mathrm{f}}\mathrm{)}\end{array}\right)& & \end{array}$(4)${{S}}_{\mathrm{on}}\mathrm{=}\left({{S}}_{\mathrm{TS}}\mathrm{(}\mathit{t}\mathrm{=}\mathrm{1}\mathrm{)}\mathrm{·}\mathrm{·}\mathrm{·}{{S}}_{\mathrm{TS}}\mathrm{(}\mathit{t}\mathrm{=}{\mathit{n}}_{\mathrm{f}}\mathrm{)}\right)\mathit{,}$(5)where the S_NGSp(t) and S_TS(t) are acquired according to Eq. (1), while S_LGSp comes from Eq. (2). In all the following, S_off and S_on will be WFS slopes measured with a disengaged loop, which are measurements of uncompensated phase.

We define ${{S}}_{\mathrm{on}}^{\mathrm{Eng}}$ as the TS measurements acquired when loop is engaged. This represents a measurement of the residual phase. Each time series of ${{S}}_{\mathrm{on}}^{\mathrm{Eng}}\mathrm{\left(}\mathit{t}\mathrm{\right)}$ gathers TS slopes following Eq. (1).

2.3. NGS/LGS-based reconstructors

2.3.1. Ground layer compensation (GLAO) reconstructor

Canary has a ground layer compensation mode (GLAO mode) with open-loop correction. It is comparable to the MOAO mode, except the reconstructor is simplified. The GLAO reconstructor averages the WFS signals, ignoring the tip and tilt from the lasers. It is computed as a sum of a tip-tilt (TT) reconstructor ${{R}}_{\mathrm{TT}}^{\mathrm{Glao}}$ that averages the slopes on NGS WFS only, and a high order reconstructor ${{R}}_{\mathrm{HO}}^{\mathrm{Glao}}$ that filters out the tip-tilt on all measurements and averages the slopes. Mathematically we have: ${{R}}^{\mathrm{Glao}}\mathrm{=}{{R}}_{\mathrm{HO}}^{\mathrm{Glao}}\mathrm{+}{{R}}_{\mathrm{TT}}^{\mathrm{Glao}}\mathit{,}$(6)where ${{R}}_{\mathrm{TT}}^{\mathrm{Glao}}$ is the n_s × (n_s × n_wfs): ${{R}}_{\mathrm{HO}}^{\mathrm{Glao}}\mathrm{=}\frac{\mathrm{1}}{{\mathit{n}}_{\mathrm{wfs}}}\begin{array}{c}\underset{\mathrm{􏽼}0\mathrm{􏽻􏽺}0\mathrm{􏽽}}{\left[\left[{{I}}_{{\mathit{n}}_{\mathrm{s}}}\mathrm{-}{F}\right]\mathrm{·}\mathrm{·}\mathrm{·}\left[{{I}}_{{\mathit{n}}_{\mathrm{s}}}\mathrm{-}{F}\right]\right]}\\ {\mathit{n}}_{\mathrm{wfs}}\end{array}\mathit{,}$(7)and ${{R}}_{\mathrm{TT}}^{\mathrm{Glao}}$ the matrix averaging the slopes on NGS WFS only is designated as: ${{R}}_{\mathrm{TT}}^{\mathrm{Glao}}\mathrm{=}\frac{\mathrm{1}}{{\mathit{n}}_{\mathrm{ngs}}}\left[\begin{array}{c}\underset{\mathrm{􏽼}0\mathrm{􏽻􏽺}0\mathrm{􏽽}}{\mathrm{\left[}{0}{\mathrm{]}}_{{\mathit{n}}_{\mathrm{s}}}\mathrm{·}\mathrm{·}\mathrm{·}\mathrm{\right[}{0}{\mathrm{]}}_{{\mathit{n}}_{\mathrm{s}}}}\\ {\mathit{n}}_{\mathrm{lgs}}\end{array}\mathit{,}\begin{array}{c}\underset{\mathrm{􏽼}0\mathrm{􏽻􏽺}0\mathrm{􏽽}}{{{F}}_{}\mathrm{·}\mathrm{·}\mathrm{·}{{F}}_{}}\\ {\mathit{n}}_{\mathrm{ngs}}\end{array}\right]\mathit{.}$(8)Calibration of references slopes, DM offsets, interaction and command matrices are described by Vidal et al. (2014).

2.3.2. Minimum mean square error (MMSE) reconstructor

The on-axis phase is determined by minimizing the distance between the on-axis measurements and the linear combination of off-axis measurements propagated through the reconstructor R: $\mathit{\epsilon }\mathrm{=}\mathrm{⟨}{\parallel {{S}}_{\mathrm{on}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{-}⟨{{S}}_{\mathrm{on}}⟩\mathrm{-}{R}\mathrm{\left(}{{S}}_{\mathrm{off}}\mathrm{\right(}\mathit{t}\mathrm{)}\mathrm{-}⟨{{S}}_{\mathrm{off}}⟩\mathrm{)}\parallel }^{\mathrm{2}}\mathrm{⟩}\mathit{,}$(9)where ⟨...⟩ denotes the temporal average over a time series. Ensuring ∂ε/∂R = 0 with S_off and S_on assumed to be Gaussian stochastic processes, we get: ${{R}}_{\mathrm{MMSE}}\mathrm{=}\underset{{R}}{\mathrm{argmin}}\mathrm{\left(}\mathit{\epsilon }\mathrm{\right)}\mathrm{=}{\mathrm{𝒞}}_{\mathrm{onoff}}{\mathrm{𝒞}}_{\mathrm{offoff}}^{†}$(10)with the n_s × (n_s × n_wfs) spatial cross-covariance matrix between on/off measurements, while denotes the (n_s × n_wfs) × (n_s × n_wfs) off-axis WFS spatial covariance matrix. They are defined by: $\hspace{0.17em}\begin{array}{c}\begin{array}{c}\end{array}{\mathrm{𝒞}}_{\mathrm{onoff}}\mathrm{=}\mathrm{⟨}\mathrm{(}{{S}}_{\mathrm{on}}\mathrm{-}⟨{{S}}_{\mathrm{on}}⟩\mathrm{)}\mathrm{(}{{S}}_{\mathrm{off}}\mathrm{-}⟨{{S}}_{\mathrm{off}}⟩{\mathrm{)}}^{\mathit{t}}\mathrm{⟩}\begin{array}{c}\end{array}\\ \begin{array}{c}\end{array}{\mathrm{𝒞}}_{\mathrm{offoff}}\mathrm{=}\mathrm{⟨}\mathrm{(}{{S}}_{\mathrm{off}}\mathrm{-}⟨{{S}}_{\mathrm{off}}⟩\mathrm{)}\mathrm{(}{{S}}_{\mathrm{off}}\mathrm{-}⟨{{S}}_{\mathrm{off}}⟩{\mathrm{)}}^{\mathit{t}}\mathrm{⟩}\mathit{.}\end{array}$(11)We note that the matrix ${\mathrm{𝒞}}_{\mathrm{offoff}}^{†}$ in Eq. (10) is the generalized inverse of the covariance matrix .

2.3.3. Filtering of pupil average slope

Laser guide stars (LGS) WFS do not properly measure the pupil average slope (the overall angle of arrival), and we have explained that we subtract this from their measurements, as stated by Eq. (2). The same filtering operation on the lasers needs to be taken into account in the tomographic reconstructor in order to get rid of undetermined modes. We use a global approach (as opposed to techniques known as split tomography (Gilles & Ellerbroek 2008) where the reconstruction is performed independently for LGS and NGS).

Before applying Eq. (10), we rework the covariance matrix and apply on each of its blocks between two wave-front sensors p and p′ the following filtering operation: ${\mathit{A}}^{\mathit{t}}\mathrm{[}{\mathrm{𝒞}}_{\mathrm{offoff}}{\mathrm{]}}_{\mathit{p}{\mathit{p}}^{\mathrm{\prime }}}\hspace{0.17em}\mathit{A}\prime$(12)where A = I_ns when p is a NGS wave-front sensor, and A = I_ns−F in the case of a LGS wave-front sensor (and similarly for A′ and p′). After filtering, the covariance matrix should be understood as the covariance of measurements excluding the undetermined modes.

As a comparison with other studies addressing the same problem, we find it interesting here to point out the method used by Tallon et al. (2010), who modifies the diagonal wave-front sensor blocks of the covariance matrix using $\mathrm{[}{\mathrm{𝒞}}_{\mathrm{offoff}}{\mathrm{]}}_{\mathit{pp}}\mathrm{+}{B}$(13)with B the covariance matrix of undetermined modes measurements for a given LGS wave-front sensor p. This is equivalent to adding noise to the measured undetermined modes, which will automatically result in their filtering by the inversion process.

The inversion process is interesting and can be compared between the two methods. In the case of Tallon, the addition does not modify the rank of , while in our case filtering two degrees of freedom (tip and tilt) for each of our four LGS WFS creates eight null eigenvalues (associated to eight eigenmodes that span the subspace of our filtered modes). Therefore our matrix cannot be inverted directly. Instead, we perform a singular value decomposition (SVD, which reduces to a diagonalization in this particular case of a square, symmetric, positive definite matrix), invert all non-zero eigenvalues, and set the inverse of null ones to zero. This is how the average slope of LGS wave-front gets filtered out in our case. We believe this approach is equivalent to Tallon’s approach with the coefficients of matrix B tending towards infinity.

We have chosen this method because the coefficients on B would have been difficult to adjust experimentally. The positioning of our LGS suffers from drifts, possibly from rotations of the asterism, that translate into long-term drifts of the average slopes that are quite empirical to turn into a covariance noise. For this reason, we prefer to completely eliminate of those modes.

In addition, we should mention that, in any case, the large dispersion of the eigenvalues of (even before filtering of tilt) makes a direct inversion numerically risky, and that is it preferable to filter out the weakest eigenmodes, instead of attempting to take the inverse of their eigenvalues. This still requires an SVD.

2.4. Reconstructor identification

The matrix can not be determined on-sky when observing because the lack of photons forbids the use a TS in the scientific directions. In addition, empirical covariance matrices suffer from statistical convergence errors (Martin et al. 2012) that make the tomographic reconstructor too specific to the particular time series from which it has been computed.

For these reasons, we fit the on-sky measured to a model of this matrix detailed in Appendix A.1. It provides analytic values of and from a list of parameters:

• n_l: number of discrete layers, set ad hoc;

• h_l: altitude of layers;

• r₀(h_l): strength of layers;

• L₀(h_l): outer scale of layers;

• telescope tracking error and vibration: adding of an isoplanatic tip-tilt over all NGS WFS;

• WFS angular positions: they are sufficiently well known using encoders information from the target acquisition system;

• system mis-registration (pupil shifts, rotation and magnification): well calibrated on bench;

• centroid gain variations;

• noise slope variance for each WFS, separately identified on the S_off time series by computing its temporal autocorrelation.

In the previous list, only the altitude, strength and outer scale of layers are retrieved, along with the tracking error. Other parameters are either set ad hoc or assumed to be known well enough.

The retrieval of these parameters is ensured thanks to the Learn&Apply algorithm (Vidal et al. 2010). It is a least-square minimization of the distance between the empiric covariance and our model of that. The algorithm operates iteratively using a Levenberg-Marquardt method. For determining the turbulent profile on an engaged set of data, we estimate the disengaged TS measurements, S_on from the measured ones ${{S}}_{\mathrm{on}}^{\mathrm{Eng}}$ (see Appendix B.1). We thus have an additional WFS, centred on-axis, that allows us to increase the maximum FoV covered by tomography.

Previously we have used a model (Vidal et al. 2010, 2014) based on a Fourier transformation of the Von-Kármán spectrum to get the bi-dimensional (2D) spatial covariance of the wave-front slopes. This algorithm was particularly efficient when coupled to wave-front sensors with fixed sub-aperture sizes, because the computation using fast Fourier transforms (FFT) could be used to extract all covariance values from the same 2D map between two wave-front sensors with arbitrary separation, without requiring any interpolation of the map and provided the sub-aperture size being a multiple of the FFT sampling.

However, the cone effect inherent to LGS induces layer-projected sub-aperture sizes and pupil sizes that vary with altitude, which breaks the efficiency of the implementation when it’s a matter of computing covariances between LGS and NGS. An interpolation procedure could be implemented in order to access values that lay between covariance samples of the FFT, but unfortunately this kind of approximation does not go integrate with a fitting procedure that internally requires the computation of derivatives as finite differences of the covariance function. For this reason, we have developed a new analytic way to compute the covariances as it is presented in Appendix A.1. This method allows us to get each coefficient of the matrix independently, which is a great advantage for faster computation for parallel implementation (Gendron et al. 2014a).

3. Evaluation of the tomographic error

3.1. Raw tomographic error

Tomographic reconstruction introduces a tomographic error when compared to single conjugate AO, to be included in the MOAO residual phase variance. In addition, the calibrated reconstructor, applied on a given data set on-sky, is usually different from the MMSE one given by Eq. (10), because the turbulence may evolve between the calibration time and the acquisition time. Moreover, the identification of the model parameters assumes the turbulence is composed by a few discrete layers, which introduces additional errors.

We define e(t) as the error between the on-axis measurement and its tomographic reconstruction from the off-axis measurements: ${e}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}{{S}}_{\mathrm{on}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{-}{R}{{S}}_{\mathrm{off}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathit{,}$(14)with R the on-sky tomographic reconstructor. The tomographic error can be derived by computing the trace of the error covariance matrix C_ee: ${\mathit{C}}_{\mathrm{ee}}\mathrm{=}\mathrm{⟨}\mathrm{\left(}{e}\mathrm{\right(}\mathit{t}\mathrm{)}\mathrm{-}⟨{e}⟩\mathrm{)}\mathrm{\left(}{e}\mathrm{\right(}\mathit{t}\mathrm{)}\mathrm{-}⟨{e}⟩{\mathrm{)}}^{\mathit{t}}\mathrm{⟩}\mathit{.}$(15)By replacing e in Eq. (15) by its expression given in Eq. (14), the covariance matrix C_ee can be split into four terms as follows: ${\mathit{C}}_{\mathrm{ee}}\mathrm{=}{\mathrm{𝒞}}_{\mathrm{onon}}\mathrm{-}{\mathrm{𝒞}}_{\mathrm{onoff}}{{R}}^{\mathit{t}}\mathrm{-}{R}{\mathrm{𝒞}}_{\mathrm{onoff}}^{\mathit{t}}\mathrm{+}{R}{\mathrm{𝒞}}_{\mathrm{offoff}}{R}t$(16)where is the n_s × n_s spatial covariance matrix for the on-axis direction: ${\mathrm{𝒞}}_{\mathrm{onon}}\mathrm{=}\mathrm{⟨}\mathrm{\left(}{{S}}_{\mathrm{on}}\mathrm{\right(}\mathit{t}\mathrm{)}\mathrm{-}⟨{{S}}_{\mathrm{on}}⟩\mathrm{)}\mathrm{(}{{S}}_{\mathrm{on}}\mathrm{-}⟨{{S}}_{\mathrm{on}}⟩{\mathrm{)}}^{\mathit{t}}\mathrm{⟩}\mathit{.}$(17)In Eq. (11), we define and . We see in Eq. (16) that the reduction of the tomographic error is given by the two terms involving the covariance matrix between the on-axis and the off-axis measurements. If the correlation is very low between these measurements, in the case of strong anisoplanatism, the on-axis error is not reduced by these two terms but instead increased by the term : the tomographic reconstruction should be avoided. This may appear when considering the highest altitude turbulence layers. On the contrary, for low altitude layers, WFS measurements are strongly correlated and we are able to significantly reduce the on-axis WF variance by the tomographic reconstruction.

For the purpose of the WF error breakdown analysis in this paper, we will reconstruct the slopes on the Zernike basis. For Canary, we consider a reconstruction of the polynomials 2 to 36 (Vidal et al. 2014). We denote M_rz the corresponding reconstruction matrix. The WF error is then defined by the quadratic sum of the Zernike coefficients. We get the WF error ${\mathit{\sigma }}_{\mathrm{2}}^{}\mathrm{\left(}{\mathit{C}}_{\mathrm{ee}}\mathrm{\right)}$, denoted as the raw tomographic error, from the trace of the Zernike coefficient covariance matrix given by: ${\mathit{\sigma }}_{\mathrm{2}}^{}\mathrm{\left(}{\mathit{C}}_{\mathrm{ee}}\mathrm{\right)}\mathrm{=}\mathrm{tr}\mathrm{(}{{M}}_{\mathrm{rz}}{\mathit{C}}_{\mathrm{ee}}{{{M}}_{\mathrm{rz}}^{\mathit{t}}}^{\mathrm{)}}\mathit{.}$(18)

3.2. Separating the aliasing error from the tomographic one

In Eq. (16), we compute the raw tomographic error covariance matrix which may include a number of error terms like noise, aliasing and tomography. In this section, we will focus on the two spatial error terms related to the turbulence: the first linked to the low spatial frequencies of the turbulence well measured by the WFS and well compensated by the DM and the second linked to the high frequencies aliased in the measurements because of the spatial sampling by the WFS. A detailed analysis of aliasing in tomography can be found in Quirós-Pacheco et al. (2010). To build the error breakdown, we want to separate the aliasing error from the tomographic one. For that purpose, we compute the aliasing error covariance matrix as: ${\mathrm{𝒞}}_{\mathrm{ee}}^{\mathrm{Alias}}\mathrm{=}{\mathrm{𝒞}}_{\mathrm{onon}}^{\mathrm{Alias}}\mathrm{-}{\mathrm{𝒞}}_{\mathrm{onoff}}^{\mathrm{Alias}}{{R}}^{\mathit{t}}\mathrm{-}{R}\mathrm{(}{\mathrm{𝒞}}_{\mathrm{onoff}}^{\mathrm{Alias}}{\mathrm{)}}^{\mathit{t}}\mathrm{+}{R}{\mathrm{𝒞}}_{\mathrm{offoff}}^{\mathrm{Alias}}{R}t$(19)where the covariance matrices, ${\mathrm{𝒞}}_{\mathrm{onon}}^{\mathrm{Alias}}$, ${\mathrm{𝒞}}_{\mathrm{onoff}}^{\mathrm{Alias}}$ and ${\mathrm{𝒞}}_{\mathrm{offoff}}^{\mathrm{Alias}}$, are derived by filtering in the turbulence spectrum all the low spatial frequencies below the DM cut-off frequency 1/2d with d the DM actuators pitch, as described in Appendix A.2. In Eq. (19), the first term is related to the WFS placed on axis, the TS for Canary. The fourth term is the error due to the aliasing of the off-axis WFS propagated through the tomographic reconstruction. The two other terms are the correlation between the aliasing in the TS and in the off-axis WFS. This correlation is only significant for the ground layer where the same turbulence is sampled and measured by all the WFS. Hence these two terms reduce the aliasing error in the measurement of the TS. We will make use of Eq. (19) to derive the aliasing contributors in the error breakdown established in the next Section.

We now derive the tomographic error subtracted by the aliasing error. We define the tomographic error covariance matrix expressed for the low spatial frequencies of the turbulence (the WF low orders controlled by the DM) as: ${\mathit{C}}_{\mathrm{ee}}^{\mathrm{\parallel }}\mathrm{=}{\mathit{C}}_{\mathrm{ee}}\mathrm{-}{\mathrm{𝒞}}_{\mathrm{ee}}^{\mathrm{Alias}}\mathit{.}$(20)And using Eqs. (16) and (19), we find: ${\mathit{C}}_{\mathrm{ee}}^{\mathrm{\parallel }}\mathrm{=}{\mathrm{𝒞}}_{\mathrm{onon}}^{\mathrm{\parallel }}\mathrm{-}{\mathrm{𝒞}}_{\mathrm{onoff}}^{\mathrm{\parallel }}{{R}}^{\mathit{t}}\mathrm{-}{R}\mathrm{(}{\mathrm{𝒞}}_{\mathrm{onoff}}^{\mathrm{\parallel }}{\mathrm{)}}^{\mathit{t}}\mathrm{+}{R}{\mathrm{𝒞}}_{\mathrm{offoff}}^{\mathrm{\parallel }}{{R}}^{\mathit{t}}\mathit{.}$(21)All these new covariance matrices C^∥ are deduced from the subtraction of the slope covariance matrices, defined in Sect. 3.1 and computed as presented in Appendix A.1, by the aliasing covariance matrices defined above. Hence ${\mathit{C}}_{\mathrm{ee}}^{\mathrm{\parallel }}$ can be fully computed from the identified model of atmosphere, as presented in Sect. 2.4, here excluding the noise contributions. Using Eq. (18), we find the tomographic error of interest: ${\mathit{\sigma }}_{\mathrm{Tomography}}^{\mathrm{2}}\mathrm{=}{\mathit{\sigma }}_{\mathrm{2}}^{}\mathrm{\left(}{\mathit{C}}_{\mathrm{ee}}^{\mathrm{\parallel }}\mathrm{\right)}\mathit{.}$(22)To further describe the tomographic error decomposition, Gendron et al. (2014b) has introduced the vertical error distribution (VED) as the decomposition of the WF tomographic error with altitude. Thanks to the statistical independence of the turbulent layers, the total tomographic error is the sum of the tomographic error on each turbulent layer at altitude h_l: $\begin{array}{ccc}\hspace{0.17em}\begin{array}{c}\begin{array}{c}\end{array}{\mathit{\sigma }}_{}^{\mathrm{2}}\mathrm{\left(}{\mathit{C}}_{\mathrm{ee}}^{\mathrm{\parallel }}\mathrm{\right)}\mathrm{=}\sum _{\mathit{l}\mathrm{=}\mathrm{1}}^{{\mathit{n}}_{\mathit{l}}}\mathrm{tr}\mathrm{(}{{M}}_{\mathrm{rz}}{\mathit{C}}_{\mathrm{ee}}^{\mathrm{\parallel }}\mathrm{\left(}{\mathit{h}}_{\mathrm{l}}\mathrm{\right)}{{{M}}_{\mathrm{rz}}^{\mathit{t}}}^{\mathrm{)}}\begin{array}{c}\end{array}\\ \begin{array}{c}\end{array}\mathrm{=}\sum _{\mathit{l}\mathrm{=}\mathrm{1}}^{{\mathit{n}}_{\mathit{l}}}{\mathit{\sigma }}_{}^{\mathrm{2}}\mathrm{\left(}{\mathit{C}}_{\mathrm{ee}}^{\mathrm{\parallel }}\mathrm{\right(}{\mathit{h}}_{\mathrm{l}}\mathrm{\left)}\mathrm{\right)}\mathit{.}\end{array}& & \end{array}$(23)Such a decomposition allows us to identify which altitude layer contributes mostly to the tomographic error, and gives us the ability to quantify the robustness of reconstructor to profile variability (Gendron et al. 2014b).

4. Error breakdown in MOAO

The design of future AO-assisted MOS for ELTs will require detailed numerical simulations of AO systems with a very large number of degrees of freedom. To prepare those simulations, it is necessary to validate the current modelling of the MOAO systems and therefore to build an accurate WF error breakdown. Canary provides two sources of performance diagnostics: the residual phase as measured by the TS and the sky point spread function (PSF) as delivered by the imaging camera Camicaz.

We propose two methods to evaluate the WF variance ${\mathit{\sigma }}_{\mathit{\epsilon }}^{\mathrm{2}}$ as observed by the science camera. The first one as the independent terms (IT) method gives ${\mathit{\sigma }}_{\mathit{\epsilon }}^{\mathrm{2}}$ from Canary telemetry data by splitting the residual phase into a sum of terms assumed to be uncorrelated. We compare the IT approach with the TS method, as the direct evaluation of ${\mathit{\sigma }}_{\mathit{\epsilon }}^{\mathrm{2}}$ based on the TS measurements. Our goal is then to compare the two methods on a large number of data sets acquired on-sky by Canary and to demonstrate that each of them leads to a very close evaluation of ${\mathit{\sigma }}_{\mathit{\epsilon }}^{\mathrm{2}}$.

4.1. Analytic method: the IT method

The IT method consists of splitting the residual phase variance ${\mathit{\sigma }}_{{\mathit{\epsilon }}_{\mathrm{IT}}}^{\mathrm{2}}$ in a sum of error terms assumed to be statistically independent and coming from the AO system performance known limitations: $\begin{array}{ccc}\hspace{0.17em}\begin{array}{c}\begin{array}{c}\end{array}{\mathit{\sigma }}_{{\mathit{\epsilon }}_{\mathrm{IT}}}^{\mathrm{2}}\mathrm{=}\hspace{0.17em}{\mathit{\sigma }}_{\mathrm{Tomography}}^{\mathrm{2}}\mathrm{+}{\mathit{\sigma }}_{\mathrm{Aliasing}}^{\mathrm{2}}\mathrm{+}{\mathit{\sigma }}_{\mathrm{Noise}}^{\mathrm{2}}\mathrm{+}{\mathit{\sigma }}_{\mathrm{Servo}}^{\mathrm{2}}\begin{array}{c}\end{array}\\ \begin{array}{c}\end{array}\mathrm{+}{\mathit{\sigma }}_{\mathrm{Fitting}}^{\mathrm{2}}\mathrm{+}{\mathit{\sigma }}_{\mathrm{Go}\mathrm{-}\mathrm{to}}^{\mathrm{2}}\mathrm{+}{\mathit{\sigma }}_{\mathrm{Static}}^{\mathrm{2}}\mathrm{+}{\mathit{\sigma }}_{\mathrm{NCPA}}^{\mathrm{2}}\mathit{,}\end{array}& & \end{array}$(24)where:

• ${\mathit{\sigma }}_{\mathrm{Tomography}}^{\mathrm{2}}$ is the tomographic error on the WF low orders compensated by the system;

• ${\mathit{\sigma }}_{\mathrm{Aliasing}}^{\mathrm{2}}$ is the aliasing error from the off-axis WFS propagated through R;

• ${\mathit{\sigma }}_{\mathrm{Noise}}^{\mathrm{2}}$ is the noise error from the off-axis WFS propagated through R and the MOAO loop;

• ${\mathit{\sigma }}_{\mathrm{Servo}}^{\mathrm{2}}$ is the temporal error due to the system transfer function including all delay;

• ${\mathit{\sigma }}_{\mathrm{Fitting}}^{\mathrm{2}}$ is the uncorrected high order error;

• ${\mathit{\sigma }}_{\mathrm{Go}\mathrm{-}\mathrm{to}}^{\mathrm{2}}$ is the go-to error of the DM;

• ${\mathit{\sigma }}_{\mathrm{Static}}^{\mathrm{2}}$ is the on-axis uncontrolled and uncalibrated static, or quasi-static error;

• ${\mathit{\sigma }}_{\mathrm{NCPA}}^{\mathrm{2}}$ is the non common path aberrations (NCPA) residuals.

Each of those variances involved in Eq. (24) are mostly computed using analytical equations, listed in the following subsections, depending on parameters deduced from the calibration data and the RTC telemetry data, which are S_off(t) and ${{S}}_{\mathrm{on}}^{\mathrm{Eng}}\mathrm{\left(}\mathit{t}\mathrm{\right)}$ when the loop is engaged.

4.1.1. Tomographic error ${\mathit{\sigma }}_{\mathrm{Tomography}}^{\mathrm{2}}$

The tomographic error is the error on the WF low orders compensated by the system, linked to the tomographic process, the used reconstructor R and the ${\mathit{C}}_{\mathit{n}}^{\mathrm{2}}\mathrm{\left(}\mathit{h}\mathrm{\right)}$ profile encountered during the observation. This error ${\mathit{\sigma }}_{\mathrm{Tomography}}^{\mathrm{2}}$ is given by Eq. (22). The interest of our approach is to evaluate it without limiting statistic convergence effects specific to the particular acquired time series involved. Compared to Vidal et al. (2014), the TS measurement is no longer required. It is only needed to compute the error covariance matrix for the low order modes, ${\mathit{C}}_{\mathrm{ee}}^{\mathrm{\parallel }}$, from the identified model of the atmosphere free from any noise and aliasing contribution (see Sect. 3.2).

4.1.2. Aliasing error ${\mathit{\sigma }}_{\mathrm{Aliasing}}^{\mathrm{2}}$

To evaluate the aliasing term, we only need to consider the off-axis aliasing error propagation through the reconstructor on the WF observed by the IR camera. From Eq. (19), we can extract this term, i.e. ${R}{\mathrm{𝒞}}_{\mathrm{offoff}}^{\mathrm{Alias}}{{R}}^{\mathit{t}}$. The aliasing error variance is then given by: ${\mathit{\sigma }}_{\mathrm{Aliasing}}^{\mathrm{2}}\mathrm{=}\mathrm{tr}\mathrm{(}{{M}}_{\mathrm{rz}}{R}{\mathrm{𝒞}}_{\mathrm{offoff}}^{\mathrm{Alias}}{{R}}^{\mathit{t}}{{{M}}_{\mathrm{rz}}^{\mathit{t}}}^{\mathrm{)}}\mathit{.}$(25)The off-axis aliasing is also time-filtered by the system, which is filtering out the high temporal frequencies. It is reasonable, for a first approximation, to neglect the impact of the filtering on the aliasing error variance, which means to neglect the temporal high frequencies from the aliasing with regards to the others terms of the error breakdown.

4.1.3. Noise error ${\mathit{\sigma }}_{\mathrm{Noise}}^{\mathrm{2}}$

As the aliasing, the noise in the residual WF error is due to the noise of the off-axis WFS propagated through the reconstructor R and the system. Considering the noise as a white process, the propagation through the loop can be readily taken into account. We denote ${\mathit{\sigma }}_{\mathrm{Noise}}^{\mathrm{2}}$ as the variance noise contribution in the residual error variance. ${\mathit{\sigma }}_{\mathrm{Noise}}^{\mathrm{2}}$ is derived by: ${\mathit{\sigma }}_{\mathrm{Noise}}^{\mathrm{2}}\mathrm{=}\mathit{\eta }\mathrm{\times }\mathrm{tr}\mathrm{(}{{M}}_{\mathrm{rz}}{R}{\mathrm{𝒞}}_{\mathrm{offoff}}^{\mathrm{Noise}}{{R}}^{\mathit{t}}{{{M}}_{\mathrm{rz}}^{\mathit{t}}}^{\mathrm{)}}\mathit{,}$(26)with: $\mathit{\eta }\mathrm{=}\frac{\mathit{g}}{\mathrm{2}\mathrm{-}\mathit{g}}\mathrm{(}\mathrm{1}\mathrm{-}\mathrm{2}\mathit{g}{\mathrm{\Delta }}_{\mathit{t}}\mathrm{(}\mathrm{1}\mathrm{-}{\mathrm{\Delta }}_{\mathit{t}}\mathrm{\left)}\mathrm{\right)}\mathit{,}$(27)where g is the loop gain and Δ_t the system latency additional to the WFS exposure-time, assumed to be varying between 0 and 1 frame. From bench calibration, we get Δ_t = 0.45 frame.

The noise covariance matrix ${\mathrm{𝒞}}_{\mathrm{offoff}}^{\mathrm{Noise}}$ is assumed to be diagonal since the LGS spots elongation is negligible on Canary (Morris et al. 2013). It is identified slope by slope computing the temporal autocorrelation on the slope time series. We extract the noise variance as a Dirac value at null delay through a parabolic fit of the turbulence contribution.

4.1.4. Servo-lag error ${\mathit{\sigma }}_{\mathrm{Servo}}^{\mathrm{2}}$

Vidal et al. (2014) has proposed to base the calculation of the servo-lag error on a end-to-end approach, particularly to take into account a fractional system delay, but it does not include the propagation of the off-axis WFS measurements through the reconstructor R.

We start by modelling the time-filtering process of the MOAO loop using the correction transfer function $\begin{array}{c}{}_{\mathrm{˜}}\\ {\mathit{h}}_{\mathrm{cor}}\end{array}\mathrm{\left(}\mathit{z}\mathrm{\right)}$ derived by a z-transform (see Appendix B.1): $\begin{array}{c}\mathrm{˜}\\ {\mathit{h}}_{\mathrm{cor}}\end{array}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{=}\mathrm{1}\mathrm{-}\mathit{g}\frac{{\mathrm{\Delta }}_{\mathit{t}}\mathrm{+}\mathrm{(}\mathrm{1}\mathrm{-}{\mathrm{\Delta }}_{\mathit{t}}\mathrm{)}\mathit{z}}{\mathit{z}\mathrm{(}\mathit{z}\mathrm{-}\mathrm{1}\mathrm{+}\mathit{g}\mathrm{)}}$(28)with z = e^2πiν/ν_e where ν_e is the time frequency sampling. The digital filter $\begin{array}{c}{}_{\mathrm{˜}}\\ {\mathit{h}}_{\mathrm{cor}}\end{array}\mathrm{\left(}\mathit{z}\mathrm{\right)}$ operates on the temporal spectrum ${R}{\begin{array}{c}{}_{{˜}}\\ {S}\end{array}}_{\mathrm{off}}$. We note a_off the atmospheric off-axis modes filtered by the tomographic reconstructor and the MOAO loop. Its temporal spectrum is derived by: ${\begin{array}{c}{˜}\\ {a}\end{array}}_{\mathrm{off}}\mathrm{\left(}\mathit{\nu }\mathrm{\right)}\mathrm{=}\begin{array}{c}\mathrm{˜}\\ {\mathit{h}}_{\mathrm{cor}}\end{array}\mathrm{(}\mathit{z}\mathrm{=}{\mathrm{e}}^{\mathrm{2}\mathrm{i}\mathit{\pi \nu }\mathit{/}{\mathit{\nu }}_{\mathrm{e}}}\mathrm{)}\mathrm{\times }{{M}}_{\mathrm{rz}}\mathit{.}{R}\mathit{.}{\begin{array}{c}{˜}\\ {S}\end{array}}_{\mathrm{off}}\mathrm{\left(}\mathit{\nu }\mathrm{\right)}\mathit{.}$(29)In Eq. (29), ${\begin{array}{c}{}_{{˜}}\\ {S}\end{array}}_{\mathrm{off}}\mathrm{\left(}\mathit{\nu }\mathrm{\right)}$ includes a noise contribution from the off-axis WFS through the reconstructor and the correction transfer function $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathit{h}}_{\mathrm{cor}}\end{array}$. For de-noising the variance of a_off, we subtract the variance of the off-axis noise, $\mathrm{tr}{}^{\mathrm{(}}{M}_{\mathrm{rz}}{R}{\mathrm{𝒞}}_{\mathrm{offoff}}^{\mathrm{Noise}}{{R}}^{\mathit{t}}{{{M}}_{\mathrm{rz}}^{\mathit{t}}}^{\mathrm{)}}$, multiplied by a coefficient which is the modulus of the MOAO correction transfer function, integrated along temporal frequencies. It has been computed by Vidal et al. (2014) to be 1 + η (see Eq. (27)).

The servo-lag WF error is then derived by integrating the modal spectrum over temporal frequencies domain and modes, until the 36th one considered as the higher order to be correctable by the system. We have: $\begin{array}{ccc}\hspace{0.17em}\begin{array}{c}\begin{array}{c}\end{array}{\mathit{\sigma }}_{\mathrm{Servo}}^{\mathrm{2}}\mathrm{=}{\begin{array}{c}\mathrm{\parallel }\\ \mathrm{\parallel }\\ \mathrm{\parallel }\\ \mathrm{\parallel }\\ \mathrm{\parallel }\end{array}}^{\mathrm{\int }}\mathit{\nu }{\begin{array}{c}{˜}\\ {a}\end{array}}_{\mathrm{off}}\mathrm{\left(}\mathit{\nu }\mathrm{\right)}\mathrm{d}\mathit{\nu }\begin{array}{c}\mathrm{\parallel }\\ \mathrm{\parallel }\\ \mathrm{\parallel }\\ \mathrm{\parallel }\\ \mathrm{\parallel }\end{array}2\begin{array}{c}\end{array}\\ \begin{array}{c}\end{array} \mathrm{-}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{\eta }\mathrm{)}\mathrm{\times }\mathrm{tr}\mathrm{(}{{M}}_{\mathrm{rz}}{R}{\mathrm{𝒞}}_{\mathrm{offoff}}^{\mathrm{Noise}}{{R}}^{\mathit{t}}{{{M}}_{\mathrm{rz}}^{\mathit{t}}}^{\mathrm{)}}\mathit{.}\end{array}& & \end{array}$(30)In Eq. (30) we do not extract the contribution of the off-axis aliasing. This term comes from the vector S_off and is propagated through the tomographic reconstructor and then through the MOAO correction transfer function. This later behaves as a high-pass filter. It means we neglect the contribution from the high temporal frequencies of the aliasing. For a Taylor-like turbulence, these frequencies are related to high spatial frequencies. Considering the turbulence spatial PSD is falling as k^{− 11/3}, the high temporal frequency contribution of the aliasing we get in Eq. (30) is negligible when compared to the parallel modes.

4.1.5. Fitting error ${\mathit{\sigma }}_{\mathrm{Fitting}}^{\mathrm{2}}$

From end-to-end simulations, Vidal et al. (2014) has established for the DM of Canary the fitting error expression versus the DM actuators pitch d and the Fried’s parameter r₀: ${\mathit{\sigma }}_{\mathrm{Fitting}}^{\mathrm{2}}\mathrm{=}\mathrm{0.3125}{\left(\frac{\mathit{d}}{{\mathit{r}}_{\mathrm{0}}}\right)}^{\mathrm{5}\mathit{/}\mathrm{3}}\mathrm{·}$(31)

4.1.6. Go-to error ${\mathit{\sigma }}_{\mathrm{Go}\mathrm{-}\mathrm{to}}^{\mathrm{2}}$

Go-to errors are produced because of the open-loop control of DMs. Generally, the DM behaviour is modelled using an interaction matrix. Many ways exist to perform this calibration, and this has been investigated for the Canary case (Kellerer et al. 2012). In closed-loop AO systems, the DN mis-registration or inaccuracies in the interaction matrix calibration are compensated iteratively by the AO loop, thanks to the feedback between WFS and DM. However, this is not the case for open-loop systems such as Canary.

We therefore have an additional WF error ${\mathit{\sigma }}_{\mathrm{Go}\mathrm{-}\mathrm{to}}^{\mathrm{2}}$ that contributes to the residual phase variance. In (Kellerer et al. 2012), this go-to error is derived as a function of the commanded WF for several methods of computation of the interaction matrix.

In the Canary case, this matrix has been computed using a sinusoidal method and, from results of Kellerer et al. (2012), we propose to define the go-to error in taking 5% of the DM command. We have thus: ${\mathit{\sigma }}_{\mathrm{Go}\mathrm{-}\mathrm{to}}^{\mathrm{2}}\mathrm{=}{\mathrm{0.05}}^{\mathrm{2}}\hspace{0.17em}\mathrm{tr}\mathrm{(}{{M}}_{\mathrm{rz}}{R}{\mathrm{𝒞}}_{\mathrm{offoff}}{{R}}^{\mathit{t}}{{{M}}_{\mathrm{rz}}^{\mathit{t}}}^{\mathrm{)}}\mathit{.}$(32)This error is very small compared to the other error terms as we will see in Sect. 6.

4.1.7. Static aberration errors ${\mathit{\sigma }}_{\mathrm{Static}}^{\mathrm{2}}$

Differential quasi static aberrations may appear during an observation. A main reason is the optical beams to the wave-front sensors and to the IR camera do not exactly cross the same areas of the optical surfaces of, both the telescope and the Canary bench. They include the field aberrations of the telescope and of the Canary de-rotator and the creeping effect. This later is a drifting of the DM shape with time (Kellerer et al. 2012).

In order to mitigate the effects of static aberration error, we could use a model of the telescope field aberrations and de-rotator optical aberrations, as well as a figure sensor to monitor the DM shape; however this approach was deemed too challenging to be applied in Canary. Instead, we perform on-sky a number of calibrations and related pre-compensations of these static aberrations for the different channels both before and during the observation. The procedure is given in Vidal et al. (2014).

Nevertheless in Canary, the calibration is not perfect and we are not able to update the static aberration compensation during the observation. Hence we have to include, in the WF error breakdown, a static aberration error term ${\mathit{\sigma }}_{\mathrm{Static}}^{\mathrm{2}}$. We evaluate post-facto this error by computing: ${\mathit{\sigma }}_{\mathrm{Static}}^{\mathrm{2}}\mathrm{=}{\begin{array}{c}\mathrm{\parallel }\\ \mathrm{\parallel }\\ \mathrm{\parallel }\\ \mathrm{\parallel }\end{array}{\mathrm{⟨}{{M}}_{\mathrm{rz}}{{S}}_{\mathrm{on}}^{\mathrm{Eng}}\mathrm{(}\mathit{t}{\mathrm{)}}^{\mathrm{⟩}}}_{\mathit{t}}\begin{array}{c}\mathrm{\parallel }\\ \mathrm{\parallel }\\ \mathrm{\parallel }\\ \mathrm{\parallel }\end{array}}^{\mathrm{2}}\mathit{.}$(33)

4.1.8. NCPA residual errors ${\mathit{\sigma }}_{\mathrm{NCPA}}^{\mathrm{2}}$

Non-common path aberrations (NCPA) are calibrated in closed loop on the TS and the bench artificial sources using a phase diversity algorithm processing the IR camera images (Gratadour et al. 2013a). The measured NCPA are then applied on the DM offset voltages during the whole duration of the observation as explained by Vidal et al. (2014). The NCPA calibration errors are mainly dominated to the high spatial frequencies of the DM surface quality not correctable by the DM. They lead to a 80% best Strehl ratio (SR) as measured in the 1550 nm image at the end of the phase diversity procedure. This corresponds to an error of 115 nm root-mean square (rms).

4.2. Empirical method: the TS method

The empirical method consists in considering the TS as the primary source of information for determining the system performance. We denote ${\mathit{\sigma }}_{\mathrm{TS}}^{\mathrm{2}}$ the variance of the TS WF engaged measurements that includes both static and dynamic terms: ${\mathit{\sigma }}_{\mathrm{TS}}^{\mathrm{2}}\mathrm{=}{⟨{\begin{array}{c}\mathrm{\parallel }\\ \mathrm{\parallel }\\ \mathrm{\parallel }\\ \mathrm{\parallel }\end{array}{{M}}_{\mathrm{rz}}{{S}}_{\mathrm{on}}^{\mathrm{Eng}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\begin{array}{c}\mathrm{\parallel }\\ \mathrm{\parallel }\\ \mathrm{\parallel }\\ \mathrm{\parallel }\end{array}}^{\mathrm{2}}⟩}_{\mathit{t}}\mathit{.}$(34)We expect the TS method to give an accurate performance estimation when compared to the science image, since only a dichroic plate separates the TS from the science path. The differences between the science wave-front and that measured by the TS are the TS noise, the TS aliasing, the fitting errors due to the wave-front high frequencies and NCPA residual errors. We want to estimate the variance of the science wave-front, and it can be written as: $\begin{array}{ccc}\hspace{0.17em}\begin{array}{c}\begin{array}{c}\end{array}{\mathit{\sigma }}_{{\mathit{\epsilon }}_{\mathrm{TS}}}^{\mathrm{2}}\mathrm{=}{\mathit{\sigma }}_{\mathrm{TS}}^{\mathrm{2}}\mathrm{-}{\mathit{\sigma }}_{\mathrm{NoiseTS}}^{\mathrm{2}}\mathrm{-}{\mathit{\sigma }}_{\mathrm{AliasingTS}}^{\mathrm{2}}\begin{array}{c}\end{array}\\ \begin{array}{c}\end{array} \mathrm{+}{\mathit{\sigma }}_{\mathrm{Fitting}}^{\mathrm{2}}\mathrm{+}{\mathit{\sigma }}_{\mathrm{NCPA}}^{\mathrm{2}}\end{array}& & \end{array}$(35)where ${\mathit{\sigma }}_{\mathrm{NoiseTS}}^{\mathrm{2}}$ is the variance of the noise and ${\mathit{\sigma }}_{\mathrm{AliasingTS}}^{\mathrm{2}}$ the variance of the aliasing, in the TS measurements.

The parameter ${\mathit{\sigma }}_{\mathrm{NoiseTS}}^{\mathrm{2}}$ is estimated using the noise covariance matrix ${\mathrm{𝒞}}_{\mathrm{onon}}^{\mathrm{Noise}}$, identified using the method based on the temporal autocorrelation of the time series of ${{S}}_{\mathrm{on}}^{\mathrm{Eng}}$. Then, the TS noise variance term is given by: ${\mathit{\sigma }}_{\mathrm{NoiseTS}}^{\mathrm{2}}\mathrm{=}\mathrm{tr}\mathrm{(}{{M}}_{\mathrm{rz}}{\mathrm{𝒞}}_{\mathrm{onon}}^{\mathrm{Noise}}{{{M}}_{\mathrm{rz}}^{\mathit{t}}}^{\mathrm{)}}\mathit{.}$(36)The aliasing included into the TS measurements is firstly given by the conventional term linked to ${\mathrm{𝒞}}_{\mathrm{onon}}^{\mathrm{Alias}}$ and is partially correlated to the off-axis aliasing propagated through the tomographic reconstructor. As presented in Sect. 4.1.1, the TS aliasing variance term is given by: ${\mathit{\sigma }}_{\mathrm{AliasingTS}}^{\mathrm{2}}\mathrm{=}\mathrm{tr}\mathrm{\left(}{{M}}_{\mathrm{rz}}\mathrm{\right(}{\mathrm{𝒞}}_{\mathrm{onon}}^{\mathrm{Alias}}\mathrm{-}{\mathrm{𝒞}}_{\mathrm{onoff}}^{\mathrm{Alias}}{{R}}^{\mathit{t}}\mathrm{-}{R}\mathrm{\left(}{\mathrm{𝒞}}_{\mathrm{onoff}}^{\mathrm{Alias}}{\mathrm{)}}^{\mathit{t}}\mathrm{\right)}{{{M}}_{\mathrm{rz}}^{\mathit{t}}}^{\mathrm{)}}\mathit{.}$(37)In this expression, the correlated part of the aliasing due to the ground layer is directly taken into account without any additional computation, compared to the previous approach given in Vidal et al. (2014).

4.3. SCAO error breakdown for CANARY

Error breakdown in SCAO has been investigated for over 20 yr (see Rigaut et al. 1991; Gendron & Léna 1994). In SCAO, the TS is obviously available and it makes sense to estimate the residual phase variance using the TS measurements. We can directly use Eq. (35), but here not subtracting the term ${\mathit{\sigma }}_{\mathrm{AliasingTS}}^{\mathrm{2}}$ which has to be kept in ${\mathit{\sigma }}_{\mathrm{TS}}^{\mathrm{2}}$.

In addition, we want to have the same breakdown as in Eq. (24), but specific to SCAO. In SCAO, we will not have tomographic, static aberrations and go-to errors. The residual phase variance can be expressed as follows in SCAO mode: $\begin{array}{ccc}\hspace{0.17em}\begin{array}{c}\begin{array}{c}\end{array}{\mathit{\sigma }}_{{\mathit{\epsilon }}_{\mathrm{IT}}}^{\mathrm{2}}\mathrm{=}{\mathit{\sigma }}_{\mathrm{AliasingTS}}^{\mathrm{2}}\mathrm{+}{\mathit{\sigma }}_{\mathrm{NoiseTS}}^{\mathrm{2}}\mathrm{+}{\mathit{\sigma }}_{\mathrm{Servo}}^{\mathrm{2}}\begin{array}{c}\end{array}\\ \begin{array}{c}\end{array} \mathrm{+}{\mathit{\sigma }}_{\mathrm{Fitting}}^{\mathrm{2}}\mathrm{+}{\mathit{\sigma }}_{\mathrm{NCPA}}^{\mathrm{2}}\mathit{.}\end{array}& & \end{array}$(38)Values of the fitting error ${\mathit{\sigma }}_{\mathrm{Fitting}}^{\mathrm{2}}$ and NCPA error ${\mathit{\sigma }}_{\mathrm{NCPA}}^{\mathrm{2}}$ are discussed respectively in Sects. 4.1.5 and 4.1.8.

We compute the aliasing error term ${\mathit{\sigma }}_{\mathrm{AliasingTS}}^{\mathrm{2}}$ by projecting the covariance matrix of the aliasing in open-loop, derived from known turbulence parameters, onto the Zernike basis: ${\mathit{\sigma }}_{\mathrm{AliasingTS}}^{\mathrm{2}}\mathrm{=}\mathrm{tr}\mathrm{(}{{M}}_{\mathrm{rz}}{\mathrm{𝒞}}_{\mathrm{onon}}^{\mathrm{Alias}}{{{M}}_{\mathrm{rz}}^{\mathit{t}}}^{\mathrm{)}}\mathit{.}$(39)In Eq. (39), we neglect the propagation of the aliasing through the closed loop and the cross-correlation with parallel modes as well. But once again, a posteriori over the large amount of SCAO data we have processed, this approximation is not critical.

The noise error ${\mathit{\sigma }}_{\mathrm{NoiseTS}}^{\mathrm{2}}$ is retrieved by the temporal auto-correlation function of the TS measurements. To make this estimation more accurate, we handle the estimation of the TS slopes in loop disengaged, S_on from the engaged ones ${{S}}_{\mathrm{on}}^{\mathrm{Eng}}$ (see Appendix B.2). The value of the noise we estimate is then multiplied by a propagation factor 1 + g/ (2−g)(1−2gΔ_t(1−Δ_t)) (Vidal et al. 2014) to take into account both the noise level at the sequence t, and the one from the sequence t−1 that has been propagating through the loop to be measured again by the TS.

Finally, the servo-lag error ${\mathit{\sigma }}_{\mathrm{Servo}}^{\mathrm{2}}$ is derived identically to that in Eq. (30): ${\mathit{\sigma }}_{\mathrm{Servo}}^{\mathrm{2}}\mathrm{=}{\mathrm{\int }}_{\mathit{\nu }}{\parallel {{a}}_{\mathrm{on}}\parallel }^{\mathrm{2}}\mathrm{\left(}\mathit{\nu }\mathrm{\right)}\mathrm{d}\mathit{\nu }\mathrm{-}{\mathit{\sigma }}_{\mathrm{NoiseTS}}^{\mathrm{2}}\mathit{,}$(40)where ${\begin{array}{c}{\mathrm{˜}}_{}\\ {a}\end{array}}_{\mathrm{on}}$ is estimated from the TS measurements: ${\begin{array}{c}{˜}\\ {a}\end{array}}_{\mathrm{on}}\mathrm{\left(}\mathit{\nu }\mathrm{\right)}\mathrm{=}\begin{array}{c}\mathrm{˜}\\ {\mathit{h}}_{\mathrm{cor}}\end{array}\mathrm{\left(}\mathit{\nu }\mathrm{\right)}\mathrm{\times }{{M}}_{\mathrm{rz}}{\begin{array}{c}{˜}\\ {S}\end{array}}_{\mathrm{on}}\mathrm{\left(}\mathit{\nu }\mathrm{\right)}\mathit{,}$(41)where $\begin{array}{c}{}_{\mathrm{˜}}\\ {\mathit{h}}_{\mathrm{cor}}\end{array}\mathrm{\left(}\mathit{\nu }\mathrm{\right)}$ is the correction transfer function given by Gendron & Léna (1994) for SCAO systems.

4.4. Evaluation of the Strehl ratio

SR is known to be related to the phase variance by the Maréchal approximation: ${\begin{array}{c}\mathrm{􏽣}\\ \mathrm{SR}\end{array}}_{\mathrm{mar}}\mathrm{=}\exp \mathrm{(}\mathrm{-}\mathrm{(}\mathrm{2}\mathit{\pi }\mathit{/}\mathit{\lambda }{\mathrm{)}}^{\mathrm{2}}{\mathit{\sigma }}_{\mathit{\epsilon }}^{\mathrm{2}}\mathrm{)}\mathit{,}$(42)where λ is the wavelength of the science image and ${\mathit{\sigma }}_{\mathit{\epsilon }}^{\mathrm{2}}$ can be either ${\mathit{\sigma }}_{{\mathit{\epsilon }}_{\mathrm{IT}}}^{\mathrm{2}}$ or ${\mathit{\sigma }}_{{\mathit{\epsilon }}_{\mathrm{TS}}}^{\mathrm{2}}$. Such an approximation is valid for small variance $\mathrm{(}{\mathit{\sigma }}_{\mathit{\epsilon }}^{\mathrm{2}}\mathit{<}\mathrm{1}{\mathrm{rad}}^{\mathrm{2}}$), or high SR greater than 10–20%. MOAO is not designed to get very high performance and the Maréchal approximation is not be accurate enough for our purposes. We propose to use the Parenti & Sasiela (1994) heuristic formula: ${\begin{array}{c}\mathrm{􏽣}\\ \mathrm{SR}\end{array}}_{\mathrm{par}}\mathrm{=}\frac{{\mathrm{e}}^{\mathrm{-}\mathrm{(}\mathrm{2}\mathit{\pi }\mathit{/}\mathit{\lambda }{\mathrm{)}}^{\mathrm{2}}{\mathit{\sigma }}_{{\mathit{\epsilon }}_{\mathrm{HO}}}^{\mathrm{2}}}}{\mathrm{1}\mathrm{+}\mathrm{(}\mathrm{2}\mathit{\pi }\mathit{/}\mathit{\lambda }{\mathrm{)}}^{\mathrm{2}}{\mathit{\sigma }}_{{\mathit{\epsilon }}_{\mathrm{TT}}}^{\mathrm{2}}}\mathrm{+}\frac{\mathrm{1}\mathrm{-}{\mathrm{e}}^{\mathrm{-}\mathrm{(}\mathrm{2}\mathit{\pi }\mathit{/}\mathit{\lambda }{\mathrm{)}}^{\mathrm{2}}{\mathit{\sigma }}_{{\mathit{\epsilon }}_{\mathrm{HO}}}^{\mathrm{2}}}}{\mathrm{1}\mathrm{+}\mathrm{(}\mathit{D}\mathit{/}{\mathit{r}}_{\mathrm{0}}\mathrm{(}\mathit{\lambda }\mathrm{)}{\mathrm{)}}^{\mathrm{2}}}$(43)where ${\mathit{\sigma }}_{{\mathit{\epsilon }}_{\mathrm{TT}}}^{\mathrm{2}}$ and ${\mathit{\sigma }}_{{\mathit{\epsilon }}_{\mathrm{HO}}}^{\mathrm{2}}$ are the residual WF error on respectively the tip-tilt only and the higher orders (tip-tilt removed).

Contrary to the Maréchal approach, the SR estimated from Eq. (43) includes two terms. The first one is the contribution of the central core of the partially AO-corrected PSF but broadened by the residual TT variance, while the second one is the contribution of the PSF broad halo. The Parenti approximation, as we will see, gives a better estimation for SR lower than 30% than the Maréchal one, which is exactly the range of performance we have encountered with Canary.

5. On-sky validation of MOAO

5.1. Observation conditions

We now focus on the Canary on-sky results acquired during the observing runs of phase B. Our purpose is firstly to prove the feasibility of managing a mixed LGS plus NGS MOAO system, secondly to evaluate the contribution of the tomographic reconstruction in the final performance of Canary. We have compared SR and error breakdown between several AO modes: SCAO (closed-loop on TS), NGS MOAO, NGS plus LGS MOAO, TT plus LGS MOAO and GLAO (ground layer compensation only).

We have simultaneously recorded system telemetry (2048 frames at 150 Hz) and H-band images. Each science image is the sum of 15 individual exposures of one second, with no recentring in order to keep tip-tilt errors between individual exposures.

During the phase B runs, we have mostly observed three different NGS asterisms, that are presented in Table 2. In addition to the NGS, we also made use of four Rayleigh LGS. They are placed on a square asterism centred on-axis. Each LGS is located at 22.6′′ off-axis from the TS position. The distance of LGS, as defined by the range gate on the LGS WFS, is 21 km, regardless of the airmass, for all nights. The gate width corresponds to a depth of 1.5 km in the atmosphere. The number of detected photons per LGS beacon is of the order of 300 photons by sub-aperture and by frame (Morris et al. 2013).

The resolution in turbulence layer altitude and the maximum distance vary between the data sets depending on the asterism. In the phase B configuration, the minimum altitude resolution was varying between 1 and 1.5 km, while the maximum distance ranged from 18 to 28 km. This maximum distance is unaffected by the LGS asterism because of the cone effect and their relatively low altitude, but mostly determined by the off-axis distance of the NGS.

We present in Fig. 1 all the retrieved layers on-sky giving for each point their altitude layers versus their local seeing at 500 nm. Each sample is a result from the L&A algorithm applied on a single data set, over 4500 sets acquired by Canary in 2013. This figure shows the distribution of the turbulence along the altitude as a probability density. It makes a strong ground layer appear, which is mainly concentrated in the first kilometer. There is a significant strength of the turbulence between 1 and 5 km. We also find a group of high altitude layers, around 12 km spreading over 4 km. A strong layer is observed around 16 km, then the strength of the turbulence decreases at higher altitudes.

We present in Fig. 2 histogram of the seeing for the ground layer (altitude below 1 km), for the altitude layer and for the total atmosphere. The ground is defined by the sum of the contributions of the layers between 0 and 1 km. We get an altitude seeing that is relatively stable, around 0.21′′ median with a standard deviation of 0.09′′, while at the ground, the median seeing reaches 0.59′′ with a much larger standard deviation of 0.34′′. The median total seeing is 0.66′′ with a standard deviation of 0.33′′.

We conclude that the variation of seeing during our nights was mostly dominated by the variation of the ground layer strength. The seeing in altitude remains quite constant during the night while the ground evolves quickly and strongly. It changes the weights between altitude and ground layers in terms of strength. This introduces an additional error propagating through the MOAO reconstructor that can be mitigated only by updating this as often as possible. We note the GLAO reconstructor, that averages directly off-axis WFS measurements (see Eqs. (6) to (8)), is not required to be updated.

Above 20 km, we are not able to identify the altitude layers properly. In Fig. 1, we have gathered all layers retrieved above 20 km at this maximal altitude. We can determine that we have around 0.089′′ seeing in median above 20 km, that corresponds to 240 nm on a 4.2 m telescope. However, we observe very high altitude layers only 34% of the time in all the phase B observations. It means such processes, quite strong but rare, exist above the WHT and directly impact the tomographic performance of Canary.

5.2. CANARY on-sky results

We report in Fig. 3H-band image SR measured with Canary running in SCAO, MOAO and GLAO. We notice on SR three different performance regimes: globally SCAO performs better than MOAO that performs better than GLAO.

Figure 3 makes a large scattering of results in any modes appear. In SCAO, this scattering involves mainly the wind speed variability and stars magnitude difference between asterisms.

In MOAO and GLAO, the scattering is mainly due to the turbulence profile variability with time. For a given seeing, the turbulence could be characterized by very different layer relative strengths. Since GLAO is compensating the ground layer only, for a given seeing, the SR will drop down if the turbulence in altitude become stronger.Conversely, for the same seeing, the SR will increase if the ground layer dominates the turbulence.

In MOAO, this scattering is less expanded than in GLAO, as it is confirmed by Table 3, and this is what we should expect from a MOAO system: if the tomographic reconstruction is properly achieved, the MOAO system must be robust to turbulent profile variability and provides the best reachable correction for a given ${\mathit{C}}_{\mathit{n}}^{\mathrm{2}}\mathrm{\left(}\mathit{h}\mathrm{\right)}$ profile.

In Table 3, we give median values of SR for various ranges of seeing to illustrate what we have previously said. For each range, we get the best performance in SCAO and the worst in GLAO. In addition, results obtained with GLAO are more variable than MOAO and SCAO.

Below 0.6′′ of seeing, we were working with good to exceptionally good conditions. The turbulence was not dominated by the ground layer and the tomographic reconstruction has operated well in comparing MOAO and GLAO. The difference between SCAO and MOAO comes essentially from the tomographic resolution we have. At this level of seeing, it becomes very important to identify accurately the ${\mathit{C}}_{\mathit{n}}^{\mathrm{2}}\mathrm{\left(}\mathit{h}\mathrm{\right)}$ profile to take into account weak layers, that are no more dominated by the most probable layers identified in Fig. 1.

Between 0.6 and 0.8′′, we were observing in nominal observation conditions. Table 3 gives nominal performance of Canary on-sky. For seeing larger than 0.8 arcsec, the SR is dropping down drastically and for very bad conditions, seeing larger than 1.2 arcsec, SCAO, MOAO and GLAO perform at similar levels. It directly comes from our previous discussion on the seeing distribution: when the seeing is so bad, it means the turbulence is dominated by the ground layer. The tomographic reconstruction can not provide a real improvement on performance as it does for better conditions.

5.3. Statistical comparison with the IT and TS methods

In Fig. 4, we present the residual WF error σ_εIT (see Eq. (24)) versus σ_εTS (see Eq. (35)). By comparing these two quantities, we quantify the accuracy of our error breakdown proposed in Sect. 4. We get about 99% of correlation in all modes and less than 1% of the total points are further than one sigma from the y = x line. We get a larger scattering in MOAO and GLAO since the error breakdown evaluation depends directly on the ${\mathit{C}}_{\mathit{n}}^{\mathrm{2}}\mathrm{\left(}\mathit{h}\mathrm{\right)}$ profile estimation. We accumulate more potential sources of error.

A posteriori, the very high correlation between analytic calculations and measurements shows all the assumptions we have made in Sect. 4 are not limiting results at this scale of atmosphere compensation.

On top of that, Fig. 4 illustrates the identification of the turbulence has been done properly, with enough accuracy and precision to get those fitting between calculations and observations. We do not observe any bias whatever the range of residual WF error as well. We are thus able to determine the turbulence characteristics, at a same level of accuracy, whatever the seeing conditions.

These results make us confident in future developments for preparing the future MOS design simulations, and especially for Mosaic the next multi-object instrument proposed for the E-ELT.

Figures 5 and 6 show the H-band SR estimated using respectively Eqs. (42) and (43), versus the sky image SR. Both of Maréchal and Parenti approximations are using ${\mathit{\sigma }}_{{\mathit{\epsilon }}_{\mathrm{IT}}}^{\mathrm{2}}$ in Eq. (24) to estimate the SR.

The relation between the SR and the residual WF variance is not bijective: for a given variance ${\mathit{\sigma }}_{\mathit{\epsilon }}^{\mathrm{2}}$, several values of SR are possible. For instance, the tip-tilt modes make the PSF core spread. For a given residual variance, if there is more residual tip-tilt in a situation than in another, the first one will comes with a worse SR. The Parenti formulation proposed in Eq. (6) tries to take into account this modal weighting, but the real relation between SR and variance cannot be given by a simple formula.

In Fig. 5, we observe a correlation between analytic and image SR for SR larger than 20%, but find a bias for lower values of image SR. For both methods, we get very similar Pearson coefficient, but the Parenti approximation makes the under-estimation of low SR disappear as shown in Fig. 6. The scattering we observe on Figs. 5 and 6 comes from the inaccuracy of Maréchal and Parenti approximations, since we do not observe such a scattering in Fig. 4.

6. Analysis of the joint NGS plus LGS tomography performed on-sky

We focus in this section on the Canary performance to target the improvement of the atmosphere compensation, achieved thanks to both tomographic reconstruction and LGS. We compare on-sky results acquired during what we call a script. A script is five data recordings of the same sequence made of three successive AO modes observation. A sequence lasts 45 s, 15 s per AO mode. This duration has been chosen to enable us to compare Canary performance, running in different AO modes, in similar atmospheric conditions.

We compare on-sky performance achieved by SCAO, GLAO and MOAO with and without LGS. We use the letters “N” and “L” following the AO modes to define respectively NGS and LGS sensors. For instance, “MOAO4L3N” means Canary was compensating the turbulence in MOAO mode using the four LGS and three NGS.

In the following, we will present three of those scripts, numbered 274, 292 and 275. In Table 4, we report observation conditions we have identified during those scripts. The seeing and outer scale values comes from the Learn&Apply algorithm.

The wind speed has been retrieved from the full width at half maximum of the slope measurement temporal auto-correlation function, after subtraction of the noise contribution. This method provides an integrated value of the wind speed along the turbulence profile.

These scripts have been selected because of their relatively close observing conditions. The observed asterism was A47 for scripts 274 and 275 and A53 for 292 (Vidal et al. 2014).

6.1. Script 274: SCAO/MOAO 4L3N/GLAO 4L3N

In Fig. 7, we show three PSFs averaged over five H-band images acquired during script 274 in SCAO, MOAO 4L3N and GLAO 4L3N, the September 13th 2013 night. Thanks to the tomographic reconstruction, we make the SR improving from 12.9% in GLAO to 21.0% in MOAO, corresponding to a wave-front error reduction by an order of 190 nm rms between the two modes, using Eq. (42). The best SR is obtained by SCAO with a value of 28.3%.

In Fig. 8, we highlight the evolution of both analytic and IR image SR with time. Analytic SR have been derived from Eq. (43) using the residual phase variance developed in Eq. (24) (IT method). The figure displays the five samples per AO mode, acquired successively and interlaced between each other.

By interlacing different AO modes, we minimize the impact of the possible turbulence non-stationarity during the script. If a particular event happens during the script, we should observe it with each mode if it lasts for at least 45 s. This ensures that we are able to compare AO modes with nearly similar observational conditions.

Figure 8 shows that we are in the most able to reproduce analytically the SR in a satisfactory way. The evolution of the IR image SR with time is well reproduced, as the absolute value and the difference between modes. However, we have an under-estimation of the SR that becomes more severe when the SR becomes higher, especially for SCAO, thus implied by the Parenti approximation according to Fig. 6. Note from the same figure, the Maréchal approximation would give worse reproduction of low SR.

For further performance analysis, we present in Fig. 9 top, in the left panel, a comparison between the turbulent profile on which the sky reconstructor has been calibrated and the post-retrieved one on the script 274 data sets. These results are from an average of fifteen profiles, each of which has been retrieved on five layers on each data set. We obtain (by post processing) 75% of turbulence at the ground (below the first kilometer), against 70% expected by the sky reconstructor (computed 20 min. before the script was run). In addition, the 11 km altitude layer observed during the script was quite well identified in the reconstructor. However, the observed 16 km altitude layer was partially identified thanks to the high altitude layers, between 14 and 18 km, as expected by the reconstructor. Finally, it seems that our MMSE reconstructor was still relevant for the turbulence profile encountered during the script.

In the middle, left panel of Fig. 9, we present the full error breakdown for each AO mode. It clearly highlights the improvement coming from the tomographic reconstruction: the term ${\mathit{\sigma }}_{\mathrm{Tomography}}^{\mathrm{2}}$ fell from 365 nm (GLAO) to 185 nm (MOAO), while the other terms were very similar. The difference between SCAO and MOAO is also explained by the tomographic error and the static terms as well.

The fluctuations of the fitting, aliasing, servo-lag and noise errors we observe are relatively small between AO modes and mainly due to the variance on the seeing estimation. They also depend on the system temporal transfer function including the loop gain and the reconstructor used.

In both MOAO and GLAO, we get 140 nm of static errors, despite the calibrations that have been done before the script (see Sect. 4.1.7). It shows these calibrations are mandatory to operate the turbulence tomographic compensation properly.

In a further analysis of the tomographic error, we show the tomographic VED in Fig. 9 bottom, left panel. It becomes clear that the MOAO reconstructor allows better compensation of the altitude layers than the GLAO reconstructor. Note we reduce the WF error by 170 nm on the 11 km layer by using MOAO. However, the GLAO reconstructor is better than the MOAO one at compensating the very low altitude layers, below 2 km, near to be the tomographic resolution with this asterism. This is due to the robustness of MMSE reconstructor: the more layers it predicts, the less efficiently it is able to compensate individual layers (Gendron et al. 2014b). This is for the same reason that the GLAO reconstructor is the only one to get 0 nm rms at 0 km, predicting this unique layer.

6.2. Script 292: SCAO/MOAO 4L3N/GLAO 4L3N

We discuss here the results from script 292. In the same way as script 274, we compared SCAO, MOAO and GLAO, but with a mis-calibration of the MMSE reconstructor. We illustrate in Fig. 9 the turbulence profile used to compute the MOAO reconstructor, compared to the one retrieved during the script 292. According to Table 4, there was one hour between the calibration and the script. During this time, strong altitude layers were appeared at 10.5 and 13.5 km, accounting together for 35% of the total turbulence and these are not expected by the MOAO reconstructor.

The global observation conditions, between script 292 and 274, was quite similar as given in Table 4. In addition, Fig. 9 allows us to compare error breakdown between these two scripts, to highlight that we obtained the same level on each error term except for tomography. We get 285 nm rms during script 292 against 185 during script 274 for MOAO. These values are comparable since the turbulence profiles during observations, as shown in Fig. 9, were largely similar, with 60–70% of ground turbulence and the rest of the energy distributed into two altitude layers. The 215 nm rms difference in tomographic error is thus due to the mis-calibration of the tomographic reconstructor. It has led to a decrease of the IR image SR according to Table 4.

Vertical error distribution plots given in Fig. 9 also illustrate that the MOAO reconstructor has been mis-calibrated. However it still performs much better than the GLAO one, especially at high altitudes. For instance at the 11 km layer, we get similar performance in GLAO during the two scripts 274 and 292, with 275 nm and 285 nm of error respectively. We underline the tomographic errors of MOAO in script 292 are very significant for the two high altitude layers at 11 and 13.5 km with 140 nm and 235 nm respectively. Much poorer performance is shown when compared to the high altitude layers of script 274.

The MOAO reconstructor has to be calibrated several times each night. It is thus mandatory to propose innovative techniques to compute and update a MMSE reconstructor for large number of degrees of freedom in order to prepare the E-ELT. Fast computation and data transfer are required. Several works are on-going on that, in particular using GPUs (Gratadour et al. 2013b; Gratadour et al. 2012).

6.3. Script 275: SCAO/MOAO 4L3N/MOAO 3N

Script 275 was obtained directly after script 274. The turbulent profile did not evolve significantly between the two scripts. We were comparing SCAO with MOAO, with and without LGS. We aim to evaluate the impact of using LGS on the tomographic reconstruction for the Canary case. We report in Fig. 10 the error breakdown we get during this script. In MOAO 4L3N, we get the same errors compared to script 274 in Fig. 9, as we expect since the observational conditions were similar according to Table 4 and the turbulent profile was similar as well. In addition, we get very similar error values between MOAO 4L3N and MOAO 3N. The main difference comes from the tomographic error with 260 nm rms against 185 nm rms respectively in MOAO 3N and MOAO 4L3N. Overall, we see an increase of nearly 4% in IR image SR (see Table 4) when introducing LGS.

Figure 11 gives the VED of the tomographic error. Thanks to LGS, we reduce the error on every layer, even at the ground. The most substantial gain is obtained at 11km because of the relatively compact LGS asterism (LGS only 22.6′′ off-axis). While at 16 km, the low altitude of the LGS (21 km) limits the gain. The VED curve, as explained in Gendron et al. (2014b), shows the distribution of the tomographic error along the altitude. In using more WFS, we should decrease the error level in this curve at any altitude. But, because of the cone effect, in adding LGS, the decrease of the tomographic error is not homogeneous over all altitudes, we so get a better relative improvement for low altitude layers than for high altitude ones. It is for this reason that we get such a gap for the 0 and 2 km layers between MOAO 4L3N and MOAO 3N. By using sodium LGS, we could expect to reduce this performance gap at higher altitudes.

7. Conclusions

This paper gives a detailed analysis of Canary sky performance in its phase B configuration, for data sets acquired in 2013 at the WHT. In particular, we gave our method to manage the LGS tip-tilt filtering and how to compute mixed NGS plus LGS GLAO and MOAO reconstructors. We have showed the tomographic error can be computed from only a reconstructor and a covariance matrix of slopes, measured on uncompensated wave-fronts. This matrix can be produced either empirically or analytically (see Appendix A.1). We have also proposed an analytic formulation of the residual phase variance. It is split into several error terms, which are assumed to be independent. We have detailed the calculation of each terms based on either the off-axis measurements or analytic formulas. We gave also another method to get the residual phase variance directly from the TS measurements.

Using 4500 data sets, we gave statistics of Canary performance and presented average values of SR for different ranges of seeing, with 30.1, 21.4 and 17.1% of SR respectively in SCAO, MOAO and GLAO using 4 LGS and 3 NGS for median seeing.

In addition, we have characterized the turbulence over the WHT during Canary on-sky observations. This was mostly dominated by the ground turbulence (h< 1 km), which reaches 0.59′′ of seeing (500 nm) with a standard-deviation of 0.34′′, compared with 0.21′′ and 0.09′′ rms for combined higher altitude layers. The total turbulence has reached 0.66′′ of median seeing and 0.33′′ in standard-deviation. In addition, we have 0.09′′ of seeing above 20 km where Canary is not able to reconstruct the turbulence phase, but only for about 34% of the time.

We have demonstrated that our analytic error breakdown computations follow the IR image measured SR fairly closely. Between the two error computation approaches (analytic and TS based), we get a correlation of 99%, for all AO modes. We have detailed the Canary performance on particular data sets. We have showed the LGS can increase the SR by up to 4% compared to the using only the NGS, while the MOAO tomographic reconstructor allows a gain of 8.5 percentage points of SR, compared to GLAO compensation.

We have also discussed the VED of the tomographic error, evaluating the gain brought by MOAO in altitude compared to GLAO, while the GLAO nulls the error in the ground layer. Moreover, we have highlighted the static aberrations in MOAO that must be carefully calibrated and monitored during the observations.

Acknowledgments

This work was supported by CNRS, INSU, Observatoire de Paris and Université Paris Diderot-Paris 7 in France, Science and Technology Facilities Council and Univ. of Durham in UK, and European Commission (Fp7 Infrastructures 2012-1, OPTICON Grant 312430, WP1). These projects have also received the support of the A*MIDEX project (no. ANR-11-IDEX-0001-02) funded by the “Investissements d’Avenir” French Government programme, managed by the French National Research Agency (ANR). This work was also funded by the UK Science and Technology Facilities Council, grants ST/K003569/1 and ST/L00075X/1. The raw data used in this publication is available from the authors.

References

Appendix A: Slopes covariance modelling

Appendix A.1: Parallel modes covariance

We consider the ith sub-aperture of the pth WFS. In the pupil plane, we note r_i the position of the ith sub-aperture than does not depend on the WFS that are all conjugated at the pupil plane. Then, α_p will be the angular position of the pth WFS from the center of the FoV. We will finally denote ${\mathit{h}}_{\mathrm{p}}^{\mathrm{LGS}}$ the altitude focus of guide stars. We denote δ_ipl the coordinate of this sub-aperture at the lth layer in the Cartesian coordinate system of the pupil plane. We have: $\begin{array}{ccc}{{\delta }}_{\mathit{i}\mathrm{pl}}\mathrm{=}\{\hspace{0.17em}\begin{array}{c}\begin{array}{c}\end{array}{{\alpha }}_{\mathrm{p}}{\mathit{h}}_{\mathrm{l}}\mathrm{+}{{r}}_{\mathit{i}}\left(\mathrm{1}\mathrm{-}\frac{{\mathit{h}}_{\mathrm{l}}}{{\mathit{h}}_{\mathrm{p}}^{\mathrm{LGS}}}\right)\begin{array}{c}\end{array}\\ \begin{array}{c}\end{array}\mathrm{foraLGSWFSand}{\mathit{h}}_{\mathrm{l}}\mathrm{\le }{\mathit{h}}_{\mathrm{p}}^{\mathrm{LGS}}\mathit{.}\\ \begin{array}{c}\end{array}{{\alpha }}_{\mathrm{p}}{\mathit{h}}_{\mathrm{l}}\mathrm{+}{{r}}_{\mathit{i}}\mathrm{else}\end{array}& & \end{array}$(A.1)For the following theoretical developments, we define x_ipl and y_ipl the projection of δ_ipl on the pupil plane reference: ${{\delta }}_{\mathit{i}\mathrm{pl}}\mathrm{=}{\mathit{x}}_{\mathit{i}\mathrm{pl}}{x}\mathrm{+}{\mathit{y}}_{\mathit{i}\mathrm{pl}}{y}\mathit{,}$(A.2)where x and y form an orthonormal basis along directions x and y. For a Shack-Hartmann WFS, the slopes 2D map comes from the average phase gradient over the lenslet: ${S}\mathrm{\left(}{{\delta }}_{\mathit{i}\mathrm{pl}}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{{\mathit{d}}_{\mathit{i}}}{”}_{\mathrm{𝒜}}{\nabla }{\phi }\mathrm{\left(}{\mathit{u}}_{\mathit{i}}\mathit{,}{\mathit{v}}_{\mathit{i}}\mathrm{\right)}\mathrm{d}{\mathit{u}}_{\mathit{i}}\mathrm{d}{\mathit{v}}_{\mathit{i}}\mathit{,}$(A.3)where d_i is the ith sub-aperture pitch and is the square integration domain comprised between −d_i/ 2 and d_i/ 2 in x and y directions. The integration variables u_iy_i are related to the ith sub-aperture and are given by u_i = x−x_ipl and v_i = y−y_ipl.

We define the x-axis and y-axis slopes as the scalar product between S(δ_ipl) and respectively x and y. In projecting the slopes along x or y directions, the integral on sub-aperture surface in Eq. (A.3) becomes integrals along sub-aperture side: $\begin{array}{ccc}\hspace{0.17em}\begin{array}{c}\begin{array}{c}\end{array}{S}\mathrm{\left(}{{\delta }}_{\mathit{i}\mathrm{pl}}\mathrm{\right)}\mathit{.}{x}\mathrm{=}\frac{\mathrm{1}}{{\mathit{d}}_{\mathit{i}}}\underset{\mathrm{-}{\mathit{d}}_{\mathit{i}}\mathit{/}\mathrm{2}}{\overset{{\mathit{d}}_{\mathit{i}}\mathit{/}\mathrm{2}}{\mathrm{\int }}}\mathrm{d}\mathit{v}i\begin{array}{c}\end{array}\\ \begin{array}{c}\end{array} \mathrm{\times }\mathrm{(}\mathit{\phi }\mathrm{(}{\mathit{x}}_{\mathit{i}\mathrm{pl}}\mathrm{+}{\mathit{d}}_{\mathit{i}}\mathit{/}\mathrm{2}\mathit{,}{\mathit{v}}_{\mathit{i}}\mathrm{)}\mathrm{-}\mathit{\phi }\mathrm{(}{\mathit{x}}_{\mathit{i}\mathrm{pl}}\mathrm{-}{\mathit{d}}_{\mathit{i}}\mathit{/}\mathrm{2}\mathit{,}{\mathit{v}}_{\mathit{i}}{\mathrm{)}}^{\mathrm{)}}\\ \begin{array}{c}\end{array}{S}\mathrm{\left(}{{\delta }}_{\mathit{i}\mathrm{pl}}\mathrm{\right)}\mathit{.}{y}\mathrm{=}\frac{\mathrm{1}}{{\mathit{d}}_{\mathit{i}}}\underset{\mathrm{-}{\mathit{d}}_{\mathit{i}}\mathit{/}\mathrm{2}}{\overset{{\mathit{d}}_{\mathit{i}}\mathit{/}\mathrm{2}}{\mathrm{\int }}}\mathrm{d}\mathit{u}i\\ \begin{array}{c}\end{array} \mathrm{\times }\mathrm{(}\mathit{\phi }\mathrm{(}{\mathit{u}}_{\mathit{i}}\mathit{,}{\mathit{y}}_{\mathit{i}\mathrm{pl}}\mathrm{+}{\mathit{d}}_{\mathit{i}}\mathit{/}\mathrm{2}\mathrm{)}\mathrm{-}\mathit{\phi }\mathrm{(}{\mathit{u}}_{\mathit{i}}\mathit{,}{\mathit{y}}_{\mathit{i}\mathrm{pl}}\mathrm{-}{\mathit{d}}_{\mathit{i}}\mathit{/}\mathrm{2}{\mathrm{)}}^{\mathrm{)}}\mathit{.}\end{array}& & \end{array}$(A.4)We then define the spatial covariance G_xx, G_yy and G_xy as the following functions: $\begin{array}{ccc}\hspace{0.17em}\begin{array}{c}\begin{array}{c}\end{array}{\mathit{G}}_{\mathit{xx}}\mathrm{\left(}{{\Delta }}_{\mathit{ij}\mathrm{pql}}\mathrm{\right)}\mathrm{=}\mathrm{⟨}\mathrm{\left(}{S}\mathrm{\right(}{{\delta }}_{\mathit{i}\mathrm{pl}}\mathrm{\left)}\mathit{.}{x}\mathrm{\right)}\mathrm{\times }\mathrm{\left(}{S}\mathrm{\right(}{{\delta }}_{\mathit{j}\mathrm{ql}}\mathrm{)}\mathit{.}{x}{\mathrm{)}}^{\mathrm{⟩}}\begin{array}{c}\end{array}\\ \begin{array}{c}\end{array}{\mathit{G}}_{\mathit{yy}}\mathrm{\left(}{{\Delta }}_{\mathit{ij}\mathrm{pql}}\mathrm{\right)}\mathrm{=}\mathrm{⟨}\mathrm{\left(}{S}\mathrm{\right(}{{\delta }}_{\mathit{i}\mathrm{pl}}\mathrm{\left)}\mathit{.}{y}\mathrm{\right)}\mathrm{\times }\mathrm{\left(}{S}\mathrm{\right(}{{\delta }}_{\mathit{j}\mathrm{ql}}\mathrm{)}\mathit{.}{y}{\mathrm{)}}^{\mathrm{⟩}}\\ \begin{array}{c}\end{array}{\mathit{G}}_{\mathit{xy}}\mathrm{\left(}{{\Delta }}_{\mathit{ij}\mathrm{pql}}\mathrm{\right)}\mathrm{=}\mathrm{⟨}\mathrm{\left(}{S}\mathrm{\right(}{{\delta }}_{\mathit{i}\mathrm{pl}}\mathrm{\left)}\mathit{.}{x}\mathrm{\right)}\mathrm{\times }\mathrm{\left(}{S}\mathrm{\right(}{{\delta }}_{\mathit{j}\mathrm{ql}}\mathrm{)}\mathit{.}{y}{\mathrm{)}}^{\mathrm{⟩}}\mathit{,}\end{array}& & \end{array}$(A.5)where Δ_ijpql = δ_ipl−δ_jql is the separation between the ith sub-aperture of pth WFS and the gth sub-aperture of qth WFS, at altitude h_l. We develop here the method to get G_xx(Δ_ijpql) only, both G_yy(Δ_ijpql) and G_xy(Δ_ijpql) are derived by the same way. The averaged gradient is the difference of phase along each side of the sub-aperture. We then get: $\begin{array}{ccc}\hspace{0.17em}\begin{array}{c}\begin{array}{c}\end{array}{\mathit{G}}_{\mathit{xx}}\mathrm{\left(}{{\Delta }}_{\mathit{ij}\mathrm{pql}}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{{\mathit{d}}_{\mathit{i}}{\mathit{d}}_{\mathit{j}}}”\mathrm{d}{\mathit{v}}_{\mathit{i}}\mathrm{d}\mathit{v}j\begin{array}{c}\end{array}\\ \begin{array}{c}\end{array} \mathrm{\times }\mathrm{⟨}{\mathrm{(}\mathit{\phi }\mathrm{(}{\mathit{x}}_{\mathit{i}\mathrm{pl}}\mathrm{+}{\mathit{d}}_{\mathit{i}}\mathit{/}\mathrm{2}\mathit{,}{\mathit{v}}_{\mathit{i}}\mathrm{+}{\mathit{y}}_{\mathit{i}\mathrm{pl}}\mathrm{)}\mathrm{-}\mathit{\phi }\mathrm{(}{\mathit{x}}_{\mathit{i}\mathrm{pl}}\mathrm{-}{\mathit{d}}_{\mathit{i}}\mathit{/}\mathrm{2}\mathit{,}{\mathit{v}}_{\mathit{i}}\mathrm{+}{\mathit{y}}_{\mathit{i}\mathrm{pl}}{\mathrm{)}}^{\mathrm{)}}}^{}\\ \begin{array}{c}\end{array} \mathrm{\times }{\mathrm{(}\mathit{\phi }\mathrm{(}{\mathit{x}}_{\mathit{j}\mathrm{ql}}\mathrm{+}{\mathit{d}}_{\mathit{j}}\mathit{/}\mathrm{2}\mathit{,}{\mathit{v}}_{\mathit{j}}\mathrm{+}{\mathit{y}}_{\mathit{j}\mathrm{ql}}\mathrm{)}\mathrm{-}\mathit{\phi }\mathrm{(}{\mathit{x}}_{\mathit{j}\mathrm{ql}}\mathrm{-}{\mathit{d}}_{\mathit{j}}\mathit{/}\mathrm{2}\mathit{,}{\mathit{v}}_{\mathit{j}}\mathrm{+}{\mathit{y}}_{\mathit{j}\mathrm{ql}}{\mathrm{)}}^{\mathrm{)}}}^{\mathrm{⟩}}\mathit{.}\end{array}& & \end{array}$(A.6)The trick to handle Eq. (A.6) is to use the following remarkable identity: $\begin{array}{ccc}\hspace{0.17em}\begin{array}{c}\begin{array}{c}\end{array}\mathrm{2}\mathrm{(}\mathit{A}\mathrm{-}{\mathit{a}}^{\mathrm{)}}\mathrm{(}\mathit{B}\mathrm{-}{\mathit{b}}^{\mathrm{)}}\mathrm{=}\mathrm{-}\mathrm{(}\mathit{A}\mathrm{-}\mathit{B}{\mathrm{)}}^{\mathrm{2}}\mathrm{+}\mathrm{(}\mathit{A}\mathrm{-}\mathit{b}\mathrm{)}2\begin{array}{c}\end{array}\\ \begin{array}{c}\end{array}\mathrm{+}\mathrm{(}\mathit{a}\mathrm{-}\mathit{B}{\mathrm{)}}^{\mathrm{2}}\mathrm{-}\mathrm{(}\mathit{a}\mathrm{-}\mathit{b}{\mathrm{)}}^{\mathrm{2}}\mathit{,}\end{array}& & \end{array}$(A.7)that allows us to rewrite Eq. (A.6) as follows: $\begin{array}{ccc}\hspace{0.17em}\begin{array}{c}\begin{array}{c}\end{array}{\mathit{G}}_{\mathit{xx}}\mathrm{\left(}{{\Delta }}_{\mathit{ij}\mathrm{pql}}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}{\mathit{d}}_{\mathit{i}}{\mathit{d}}_{\mathit{j}}}”\mathrm{d}{\mathit{v}}_{\mathit{i}}\mathrm{d}\mathit{v}j\begin{array}{c}\end{array}\\ \begin{array}{c}\end{array} \mathrm{\times }\mathrm{⟨}\mathrm{-}\mathrm{\left(}\mathit{\phi }\mathrm{\right(}{\mathit{x}}_{\mathit{i}\mathrm{pl}}\mathrm{+}{\mathit{d}}_{\mathit{i}}\mathit{/}\mathrm{2}\mathit{,}{\mathit{v}}_{\mathit{i}}\mathrm{+}{\mathit{y}}_{\mathit{i}\mathrm{pl}}\mathrm{)}\mathrm{-}\mathit{\phi }\mathrm{(}{\mathit{x}}_{\mathit{j}\mathrm{ql}}\mathrm{+}{\mathit{d}}_{\mathit{j}}\mathit{/}\mathrm{2}\mathit{,}{\mathit{v}}_{\mathit{j}}\mathrm{+}{\mathit{y}}_{\mathit{j}\mathrm{ql}}\mathrm{)}{\mathrm{)}}^{\mathrm{2}}\\ \begin{array}{c}\end{array} \mathrm{-}\mathrm{\left(}\mathit{\phi }\mathrm{\right(}{\mathit{x}}_{\mathit{i}\mathrm{pl}}\mathrm{-}{\mathit{d}}_{\mathit{i}}\mathit{/}\mathrm{2}\mathit{,}{\mathit{v}}_{\mathit{i}}\mathrm{+}{\mathit{y}}_{\mathit{i}\mathrm{pl}}\mathrm{)}\mathrm{-}\mathit{\phi }\mathrm{(}{\mathit{x}}_{\mathit{j}\mathrm{ql}}\mathrm{-}{\mathit{d}}_{\mathit{j}}\mathit{/}\mathrm{2}\mathit{,}{\mathit{v}}_{\mathit{j}}\mathrm{+}{\mathit{y}}_{\mathit{j}\mathrm{ql}}\mathrm{\left)}\mathrm{\right)}2\\ \begin{array}{c}\end{array} \mathrm{+}\mathrm{\left(}\mathit{\phi }\mathrm{\right(}{\mathit{x}}_{\mathit{i}\mathrm{pl}}\mathrm{+}{\mathit{d}}_{\mathit{i}}\mathit{/}\mathrm{2}\mathit{,}{\mathit{v}}_{\mathit{i}}\mathrm{+}{\mathit{y}}_{\mathit{i}\mathrm{pl}}\mathrm{)}\mathrm{-}\mathit{\phi }\mathrm{(}{\mathit{x}}_{\mathit{j}\mathrm{ql}}\mathrm{-}{\mathit{d}}_{\mathit{j}}\mathit{/}\mathrm{2}\mathit{,}{\mathit{v}}_{\mathit{j}}\mathrm{+}{\mathit{y}}_{\mathit{j}\mathrm{ql}}\mathrm{\left)}\mathrm{\right)}2\\ \begin{array}{c}\end{array} \mathrm{+}\mathrm{\left(}\mathit{\phi }\mathrm{\right(}{\mathit{x}}_{\mathit{i}\mathrm{pl}}\mathrm{-}{\mathit{d}}_{\mathit{i}}\mathit{/}\mathrm{2}\mathit{,}{\mathit{v}}_{\mathit{i}}\mathrm{+}{\mathit{y}}_{\mathit{j}\mathrm{ql}}\mathrm{)}\mathrm{-}\mathit{\phi }\mathrm{(}{\mathit{x}}_{\mathit{j}\mathrm{ql}}\mathrm{+}{\mathit{d}}_{\mathit{j}}\mathit{/}\mathrm{2}\mathit{,}{\mathit{v}}_{\mathit{j}}\mathrm{+}{\mathit{y}}_{\mathit{j}\mathrm{ql}}\mathrm{)}{\mathrm{)}}^{\mathrm{2}}\mathrm{⟩}\mathit{.}\end{array}& & \end{array}$(A.8)Considering D_φ(ρ) = ^⟨(φ(r)−φ(r + ρ))^2⟩ as the phase Structure Function (SF), we have the final expression of the x-axis slopes covariance: $\begin{array}{ccc}\hspace{0.17em}\begin{array}{c}\begin{array}{c}\end{array}{\mathit{G}}_{\mathit{xx}}\mathrm{\left(}{{\Delta }}_{\mathit{ij}\mathrm{pql}}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}{\mathit{d}}_{\mathit{i}}{\mathit{d}}_{\mathit{j}}}”\mathrm{d}{\mathit{v}}_{\mathit{i}}\mathrm{d}\mathit{v}j\begin{array}{c}\end{array}\\ \begin{array}{c}\end{array} \mathrm{\times }(\mathrm{-}\mathrm{2}{\mathit{D}}_{\mathit{\phi }}\left({{\Delta }}_{\mathit{ij}\mathrm{pql}}\mathrm{+}\frac{{\mathit{d}}_{\mathit{i}}\mathrm{-}{\mathit{d}}_{\mathit{j}}}{\mathrm{2}}{x}\mathrm{+}\mathrm{(}{\mathit{v}}_{\mathit{i}}\mathrm{-}{\mathit{v}}_{\mathit{j}}\mathrm{)}{y}\right)\\ \begin{array}{c}\end{array} \mathrm{+}{\mathit{D}}_{\mathit{\phi }}\left({{\Delta }}_{\mathit{ij}\mathrm{pql}}\mathrm{+}\frac{{\mathit{d}}_{\mathit{i}}\mathrm{+}{\mathit{d}}_{\mathit{j}}}{\mathrm{2}}{x}\mathrm{+}\mathrm{(}{\mathit{v}}_{\mathit{i}}\mathrm{-}{\mathit{v}}_{\mathit{j}}\mathrm{)}{y}\right)\\ \begin{array}{c}\end{array} \mathrm{+}{\mathit{D}}_{\mathit{\phi }}\left({{\Delta }}_{\mathit{ij}\mathrm{pql}}\mathrm{-}\frac{{\mathit{d}}_{\mathit{i}}\mathrm{+}{\mathit{d}}_{\mathit{j}}}{\mathrm{2}}{x}\mathrm{+}\mathrm{(}{\mathit{v}}_{\mathit{i}}\mathrm{-}{\mathit{v}}_{\mathit{j}}\mathrm{)}{y}\right))\mathit{.}\end{array}& & \end{array}$(A.9)We now consider W_φ(k) as the Von-Kármán spatial power spectral density (PSD) given by: ${\mathit{W}}_{\mathit{\phi }}\mathrm{\left(}{k}\mathrm{\right)}\mathrm{=}\mathrm{0.023}{\mathit{r}}_{\mathrm{0}}^{\mathrm{-}\mathrm{5}\mathit{/}\mathrm{3}}{\mathrm{(}{\mathit{k}}^{\mathrm{2}}\mathrm{+}\mathrm{1}\mathit{/}{{\mathit{L}}_{\mathrm{0}}^{\mathrm{2}}}^{\mathrm{)}}}^{\mathrm{-}\mathrm{11}\mathit{/}\mathrm{6}}\mathit{,}$(A.10)The phase structure function (SF) is related to the spatial PSD W_φ(k) by: ${\mathit{D}}_{\mathit{\phi }}\left({\rho }\right)\mathrm{=}\mathrm{2}{”}_{R^{\mathrm{2}}}{\mathit{W}}_{\mathit{\phi }}\mathrm{\left(}{k}\mathrm{\right)}\left(\mathrm{1}\mathrm{-}\cos \mathrm{\left(}\mathrm{2}\mathit{\pi }{k}{\rho }\mathrm{\right)}\right)\mathrm{d}{k}\mathit{.}$(A.11)We commonly model the phase SF using the Von-Kármán expression: $\begin{array}{ccc}\hspace{0.17em}\begin{array}{c}\begin{array}{c}\end{array}{\mathit{D}}_{\mathit{\phi }}\mathrm{\left(}{\rho }\mathrm{\right)}\mathrm{=}{\mathit{k}}_{\mathrm{1}}{\left(\frac{{\mathit{L}}_{\mathrm{0}}\mathrm{\left(}{\mathit{h}}_{\mathrm{l}}\mathrm{\right)}}{{\mathit{r}}_{\mathrm{0}}\mathrm{\left(}{\mathit{h}}_{\mathrm{l}}\mathrm{\right)}}\right)}^{\mathrm{5}\mathit{/}\mathrm{3}}\mathrm{[}{\mathit{k}}_{\mathrm{2}}\mathrm{-}{\mathit{\rho }}^{\mathrm{5}\mathit{/}\mathrm{6}}{\mathit{K}}_{\mathrm{5}\mathit{/}\mathrm{6}}\mathrm{(}\mathit{\rho }{\mathrm{)}}^{\mathrm{]}}\mathit{,}\begin{array}{c}\end{array}\\ \begin{array}{c}\end{array}{\mathit{k}}_{\mathrm{1}}\mathrm{=}\frac{{\mathrm{2}}^{\mathrm{1}\mathit{/}\mathrm{6}}\mathrm{\Gamma }\mathrm{(}\mathrm{11}\mathit{/}\mathrm{6}\mathrm{)}}{{\mathit{\pi }}^{\mathrm{8}\mathit{/}\mathrm{3}}}{\left[\frac{\mathrm{24}}{\mathrm{5}}\mathrm{\Gamma }\mathrm{(}\mathrm{6}\mathit{/}\mathrm{5}\mathrm{)}\right]}^{\mathrm{5}\mathit{/}\mathrm{6}}\\ \begin{array}{c}\end{array}{\mathit{k}}_{\mathrm{2}}\mathrm{=}\frac{\mathrm{\Gamma }\mathrm{(}\mathrm{5}\mathit{/}\mathrm{6}\mathrm{)}}{{\mathrm{2}}^{\mathrm{1}\mathit{/}\mathrm{6}}}\mathit{,}\end{array}& & \end{array}$(A.12)with ρ = 2π|Δ_ijpql| /L₀(h_l). By the same way, we get: $\begin{array}{ccc}\hspace{0.17em}\begin{array}{c}\begin{array}{c}\end{array}{\mathit{G}}_{\mathit{yy}}\mathrm{\left(}{{\Delta }}_{\mathit{ij}\mathrm{pql}}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}{\mathit{d}}_{\mathit{i}}{\mathit{d}}_{\mathit{j}}}”\mathrm{d}{\mathit{u}}_{\mathit{i}}\mathrm{d}\mathit{u}j\begin{array}{c}\end{array}\\ \begin{array}{c}\end{array} \mathrm{\times }(\mathrm{-}\mathrm{2}{\mathit{D}}_{\mathit{\phi }}\left({{\Delta }}_{\mathit{ij}\mathrm{pql}}\mathrm{+}\frac{{\mathit{d}}_{\mathit{i}}\mathrm{-}{\mathit{d}}_{\mathit{j}}}{\mathrm{2}}{y}\mathrm{+}\mathrm{(}{\mathit{u}}_{\mathit{i}}\mathrm{-}{\mathit{u}}_{\mathit{j}}\mathrm{)}{x}\right)\\ \begin{array}{c}\end{array} \mathrm{+}{\mathit{D}}_{\mathit{\phi }}\left({{\Delta }}_{\mathit{ij}\mathrm{pql}}\mathrm{+}\frac{{\mathit{d}}_{\mathit{i}}\mathrm{+}{\mathit{d}}_{\mathit{j}}}{\mathrm{2}}{y}\mathrm{+}\mathrm{(}{\mathit{u}}_{\mathit{i}}\mathrm{-}{\mathit{u}}_{\mathit{j}}\mathrm{)}{x}\right)\\ \begin{array}{c}\end{array} \mathrm{+}{\mathit{D}}_{\mathit{\phi }}\left({{\Delta }}_{\mathit{ij}\mathrm{pql}}\mathrm{-}\frac{{\mathit{d}}_{\mathit{i}}\mathrm{+}{\mathit{d}}_{\mathit{j}}}{\mathrm{2}}{y}\mathrm{+}\mathrm{(}{\mathit{u}}_{\mathit{i}}\mathrm{-}{\mathit{u}}_{\mathit{j}}\mathrm{)}{x}\right))\mathit{,}\end{array}& & \end{array}$(A.13)and: $\begin{array}{ccc}\hspace{0.17em}\begin{array}{c}\begin{array}{c}\end{array}{\mathit{G}}_{\mathit{xy}}\mathrm{\left(}{{\Delta }}_{\mathit{ij}\mathrm{pql}}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}{\mathit{d}}_{\mathit{i}}{\mathit{d}}_{\mathit{j}}}”\mathrm{d}{\mathit{u}}_{\mathit{j}}\mathrm{d}\mathit{v}i\begin{array}{c}\end{array}\\ \begin{array}{c}\end{array} \mathrm{\times }(\mathrm{-}{\mathit{D}}_{\mathit{\phi }}\left({{\Delta }}_{\mathit{ij}\mathrm{pql}}\mathrm{-}\left({\mathit{u}}_{\mathit{j}}\mathrm{-}\frac{{\mathit{d}}_{\mathit{j}}}{\mathrm{2}}\right){x}\mathrm{+}\left({\mathit{v}}_{\mathit{i}}\mathrm{-}\frac{{\mathit{d}}_{\mathit{i}}}{\mathrm{2}}\right){y}\right)\\ \begin{array}{c}\end{array} \mathrm{-}{\mathit{D}}_{\mathit{\phi }}\left({{\Delta }}_{\mathit{ij}\mathrm{pql}}\mathrm{-}\left({\mathit{u}}_{\mathit{j}}\mathrm{+}\frac{{\mathit{d}}_{\mathit{j}}}{\mathrm{2}}\right){x}\mathrm{+}\left({\mathit{v}}_{\mathit{i}}\mathrm{+}\frac{{\mathit{d}}_{\mathit{i}}}{\mathrm{2}}\right){y}\right)\\ \begin{array}{c}\end{array} \mathrm{+}{\mathit{D}}_{\mathit{\phi }}\left({{\Delta }}_{\mathit{ij}\mathrm{pql}}\mathrm{-}\left({\mathit{u}}_{\mathit{j}}\mathrm{-}\frac{{\mathit{d}}_{\mathit{j}}}{\mathrm{2}}\right){x}\mathrm{+}\left({\mathit{v}}_{\mathit{i}}\mathrm{+}\frac{{\mathit{d}}_{\mathit{i}}}{\mathrm{2}}\right){y}\right)\\ \begin{array}{c}\end{array} \mathrm{+}{\mathit{D}}_{\mathit{\phi }}\left({{\Delta }}_{\mathit{ij}\mathrm{pql}}\mathrm{-}\left({\mathit{u}}_{\mathit{j}}\mathrm{+}\frac{{\mathit{d}}_{\mathit{j}}}{\mathrm{2}}\right){x}\mathrm{+}\left({\mathit{v}}_{\mathit{i}}\mathrm{-}\frac{{\mathit{d}}_{\mathit{i}}}{\mathrm{2}}\right){y}\right))\mathit{.}\end{array}& & \end{array}$(A.14)Integrals in Eqs. (A.9), (A.13) and (A.14) must be discretized along the sub-aperture side. To speed up the computation, we assume a Hudgin-like model: the phase gradient is calculated in taking only the side-to-side phase difference at the middle of the sub-apertures. We have thus u_i,u_j,v_i and v_j forced to zero and the spatial covariance is obtained in summing the phase structure functions under integrals (Gendron et al. 2014b; Martin et al. 2012).

Appendix A.2: Aliasing covariance

The covariance of the aliased phase can be derived using Eqs. (A.9), (A.13) and (A.14). These equations involve an expression of the phase SF, that depends on the turbulence characteristics and the geometry, according to Eq. (A.12). To get the aliasing contribution in the covariance, we have to split the phase into a parallel part, which is compensated by the system, and an orthogonal part: $\mathit{\phi }\mathrm{=}{\mathit{\phi }}^{\mathrm{\parallel }}\mathrm{+}{\mathit{\phi }}^{\mathrm{\perp }}\mathit{.}$(A.15)The orthogonal phase contains only spatial frequencies higher than 1/2d, the DM cut-off frequency. Considering the quantity L₀ is large compared to 2d, the aliasing SF ${\mathit{D}}_{\mathit{\phi }}^{\mathrm{\perp }}\mathrm{\left(}{\rho }\mathrm{\right)}$ is given by: ${\mathit{D}}_{\mathit{\phi }}^{\mathrm{\perp }}\mathrm{\left(}{\rho }\mathrm{\right)}\mathrm{=}\mathrm{0.046}{\mathit{r}}_{\mathrm{0}}^{\mathrm{-}\mathrm{5}\mathit{/}\mathrm{3}}\underset{\mathrm{1}\mathit{/}\mathrm{2}\mathit{d}}{\overset{ \mathrm{\infty }}{”}}{\mathit{k}}^{\mathrm{-}\mathrm{11}\mathit{/}\mathrm{3}}\mathrm{(}\mathrm{1}\mathrm{-}\cos \mathrm{(}\mathrm{2}\mathit{\pi }{k}{\rho }\mathrm{\left)}\mathrm{\right)}\mathrm{d}{k}\mathit{.}$(A.16)Using Eq. (A.16) to replace the SF expression in Eqs. (A.9), (A.13) and (A.14), we get an analytic expression of the spatial covariance of the aliasing.

Appendix B: MOAO control

Appendix B.1: MOAO control law

The residual modes in science directions are split into a compensated, or parallel part by the DM, ε_∥ and high frequencies non-reproducible, or orthogonal, by the DM, a_⊥: ${\epsilon }\mathrm{=}{{\epsilon }}_{{\parallel }}\mathrm{+}{{a}}_{{\perp }}\mathit{.}$(B.1)The turbulence compensation achieved thanks to the DM results in that direction from the subtraction between the incoming parallel modes and the mirror modes m_DM: ${{\epsilon }}_{{\parallel }}\mathrm{=}{{a}}_{{\parallel }}\mathrm{-}{{m}}_{\mathrm{DM}}\mathit{.}$(B.2)The MOAO controller determines the mirror modes m_DM to apply at a sequence t from the previous sequence t−1, the off-axis measurements S_off and M_c the modal command matrix used on-sky. This later includes both tomographic reconstruction and DM/WFS calibrations (Vidal et al. 2014). In order to get a system fractional delay 1 + Δ_t, we combine off-axis measurements acquired at t−2 and t−1 in this way: $\begin{array}{ccc}\hspace{0.17em}\begin{array}{c}\begin{array}{c}\end{array}{{m}}_{\mathrm{DM}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{g}\mathrm{)}\mathrm{\times }{{m}}_{\mathrm{DM}}\mathrm{(}\mathit{t}\mathrm{-}\mathrm{1}\mathrm{)}\begin{array}{c}\end{array}\\ \begin{array}{c}\end{array} \mathrm{+}\mathit{g}{{M}}_{\mathrm{c}}{R}\left({\mathrm{\Delta }}_{\mathit{t}}{{S}}_{\mathrm{off}}\mathrm{(}\mathit{t}\mathrm{-}\mathrm{2}\mathrm{)}\mathrm{+}\mathrm{(}\mathrm{1}\mathrm{-}{\mathrm{\Delta }}_{\mathit{t}}\mathrm{)}{{S}}_{\mathrm{off}}\mathrm{(}\mathit{t}\mathrm{-}\mathrm{1}\mathrm{)}\right)\mathit{.}\end{array}& & \end{array}$(B.3)Applying the z-transform to Eq. (B.3), we get: $\begin{array}{c}{˜}\\ {{m}}_{\mathrm{DM}}\end{array}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{=}\begin{array}{c}\mathrm{˜}\\ {\mathit{h}}_{\mathrm{ol}}\end{array}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{\times }{{M}}_{\mathrm{c}}{R}{\begin{array}{c}\mathrm{˜}\\ {S}\end{array}}_{\mathrm{off}}\mathrm{\left(}\mathit{z}\mathrm{\right)}$(B.4)where z = e^{− 2iπν/ν_e} and ν_e the sampling frequency. The transfer function $\begin{array}{c}{}_{\mathrm{˜}}\\ {\mathit{h}}_{\mathrm{ol}}\end{array}\mathrm{\left(}\mathit{z}\mathrm{\right)}$ is the MOAO controller transfer function defined as: $\begin{array}{c}\mathrm{˜}\\ {\mathit{h}}_{\mathrm{ol}}\end{array}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{=}\mathit{g}\mathrm{\times }\frac{{\mathrm{\Delta }}_{\mathit{t}}\mathrm{+}\mathrm{(}\mathrm{1}\mathrm{-}{\mathrm{\Delta }}_{\mathit{t}}\mathrm{)}\mathit{z}}{\mathit{z}\mathrm{(}\mathit{z}\mathrm{-}\mathrm{1}\mathrm{+}\mathit{g}\mathrm{)}}\mathrm{·}$(B.5)The MOAO correction transfer function of the atmospheric parallel modes, $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathit{h}}_{\mathrm{cor}}\end{array}$ is given by: $\begin{array}{c}\mathrm{˜}\\ {\mathit{h}}_{\mathrm{cor}}\end{array}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{=}\frac{\begin{array}{c}{˜}\\ {{\epsilon }}_{{\parallel }}\end{array}}{\begin{array}{c}{˜}\\ {{a}}_{{\parallel }}\end{array}}\mathrm{=}\mathrm{1}\mathrm{-}\begin{array}{c}\mathrm{˜}\\ {\mathit{h}}_{\mathrm{ol}}\end{array}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathit{.}$(B.6)

Appendix B.2: Disengaged TS measurements estimation

We note ${{S}}_{\mathrm{on}}^{\mathrm{Eng}}\mathrm{\left(}\mathit{t}\mathrm{\right)}$ the time series acquired by the TS in engaged loop. Using the voltages vector V, the interaction matrix M_int and the calibrated the fractional delay Δ_t, the disengaged TS slopes is estimated by: ${\begin{array}{c}\mathrm{􏽢}\\ {S}\end{array}}_{\mathrm{on}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}{{S}}_{\mathrm{on}}^{\mathrm{Eng}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{-}{{M}}_{\mathrm{int}}\left({\mathrm{\Delta }}_{\mathit{t}}{V}\mathrm{(}\mathit{t}\mathrm{-}\mathrm{2}\mathrm{)}\mathrm{+}\mathrm{(}\mathrm{1}\mathrm{-}{\mathrm{\Delta }}_{\mathit{t}}\mathrm{)}{V}\mathrm{(}\mathit{t}\mathrm{-}\mathrm{1}\mathrm{)}\right)\mathit{.}$(B.7)We reproduce the filtering operated by the RTC, described in Eq. (B.3), on-sky to estimate what the TS would measure if the DM stayed flat.

All Tables

All Figures