© R. J. L. Fétick et al. 2019

1. Introduction

Optical systems suffer from aberrations and diffraction effects that limit their imaging performance. For ground-based observations, the point spread function (PSF) is dramatically altered by the atmospheric turbulence that distorts the incoming wavefront (Roddier 1981). The resolution under typical conditions of a seeing-limited telescope does not exceed the diffraction limit of an ∼12 cm aperture. Modern and future large telescopes thus include adaptive optics (AO) systems (Roddier 1999) that compensate for the atmospheric turbulence thanks to wavefront sensors and deformable mirrors. The aberrated wavefront is partially corrected and telescopes may operate near their diffraction limited regime. Nevertheless the AO correction is limited by technical issues such as sensor noise, limited number of actuators, or loop delay (Martin et al. 2017; Rigaut et al. 1998). This results in a peculiar shape of the PSF made of a sharp peak due to the partial AO correction, and a wide halo caused by the residual turbulence above the AO cutoff frequency.

The PSF thus provides critical information about an optical system regarding its preliminary design, calibrations, testings, or diagnostic (Ascenso et al. 2015; Ragland et al. 2018). Image post-processing, such as deconvolution (Mugnier et al. 2004), also requires knowledge of the PSF. Deconvolution of long-exposure images using parametric PSFs has already been demonstrated in Drummond (1998) and Fétick et al. (2019). A fine model of the PSF is necessary. The substantial advantage of parametric PSFs is to compress all the important information of the physical PSF into a small number of parameters. The numerical values of these parameters might then be used for comparisons, correlations, or any statistical analysis. Moreover if the PSF parameters are correlated to physical values (e.g. turbulence strength, wind speed, AO residual phase variance), it is possible to better constrain these parameters or better understand the AO response to given observing conditions. We state that an efficient AO PSF model should fulfil the following requirements:

– Accuracy. The model must represent accurately the shape of the AO-corrected PSF, especially the two areas corresponding to its central peak and to its wide turbulent halo. The requested accuracy depends on the application of the PSF (e.g. fitting, deconvolution, turbulence monitoring).

– Versatility and robustness. The model must be used on different AO system, with different AO correction levels, for different turbulent strengths.

– Simplicity. The model must have as few parameters as possible without damaging its versatility or accuracy.

– Physical parameters. Such parameters have a physical meaning related to the observing conditions. These parameters have physical units.

The literature already provides some models of AO-corrected PSF (Drummond 1998; Zieleniewski & Thatte 2013) often based on Gaussian, Lorentzian, and/or Moffat (1969) models. A trade-off is always drawn between a simple model with few parameters but imprecise, or a more precise but also more complex model. The difficulty often comes from the description of the turbulent halo with only a few parameters. Moreover, to the best of our knowledge, these PSF models rely only on mathematical parameters without direct physical meaning or units.

We propose a long-exposure PSF model for AO-corrected telescopes that describes accurately the shape of the PSF; this model is made of a small number of parameters with physical meaning whenever possible. Our method does not parameterize the PSF directly in the focal plane, but rather from the phase power spectral density (PSD). Indeed Goodman (1968) and Roddier (1981) have shown that the phase PSD contains all the necessary information to describe the long-exposure atmospheric PSF. Working in the PSD domain allows us to include physical parameters. Then Fourier transforms give the resulting PSF in the focal plane. Our PSF model also includes pupil diffraction effects or any of the system static aberrations, provided they have been previously characterized.

In Sect. 2 we first recall the expression of the Moffat function and show its limits for AO PSF description. Then we develop our PSF model, partially based on this Moffat function. Section 3 validates the model by fitting PSFs from numerical simulations and from observations made on two Very Large Telescope (VLT) instruments. Finally Sect. 4 concludes our work and discusses direct and future applications for our PSF model.

2. Description of the PSF model

In the whole paper, we define (x_R, y_R) the reference coordinates that are respectively the detector horizontal and vertical coordinates. We also define (x, y) the proper PSF coordinates (e.g. along the major and minor PSF elongation axis), rotated by an angle θ_R with respect to the reference frame. The reference frame to PSF frame transformation can be written as

$$\begin{array}{c}\hfill \left(\begin{array}{c}x\\ y\end{array}\right)=\left(\begin{array}{cc}cos{\theta }_R& sin{\theta }_R\\ -sin{\theta }_R& cos{\theta }_R\end{array}\right)\left(\begin{array}{c}{x}_R\\ y_R\end{array}\right).\end{array}$$(1)

For the sake of simplicity we will use mainly the rotated coordinates, but it is important to keep in mind that θ_R is a crucial PSF parameter. We also simplify the notations x = x(x_R, y_R, θ_R) and y = y(x_R, y_R, θ_R).

In this section we first recall the usual Moffat PSF model, since it encompasses and generalizes Lorentzian and Gaussian models. We demonstrate the advantages of the Moffat model, but also its limitations. This motivates our search for a better PSF model. However, the mathematical expression of the Moffat function will still be used inside our more complete PSF model.

2.1. Review of the usual Moffat PSF model

The AO-corrected PSF exhibits a sharp corrected peak, with wide wings extension. The Moffat (1969) model is often used due to its good approximation of the AO PSF sharp peak (Andersen et al. 2006; Müller Sánchez et al. 2006; Davies & Kasper 2012; Orban de Xivry et al. 2015; Rusu et al. 2016). The Moffat function, of amplitude A, is written as

$$\begin{array}{c}\hfill M_{\mathrm{A}}(x,y)=\frac{A}{{(1+x^2/{\alpha }_x^2+y^2/{\alpha }_y^2)}^{\beta }},\end{array}$$(2)

with α_x, α_y, and β strictly positive real numbers. Moreover the condition β >  1 is imposed to ensure a finite integral of the function on the plane. This model encompasses the two-dimensional Lorentzian function for β = 1 and the two-dimensional Gaussian function for β → +∞. The variable β parameter thus makes the Moffat function a generalization of Lorentzian and Gaussian ones. Since the PSF has a unit energy, demonstration in Appendix A shows that the Moffat multiplicative constant is

$$\begin{array}{c}\hfill A=\frac{\beta -1}{\pi {\alpha }_x{\alpha }_y},\end{array}$$(3)

so the PSF, called h, is made of only four free parameters α_x, α_y, β, and θ_R. The Moffat PSF model thus is re-written as

$$\begin{array}{c}\hfill h(x,y)=\frac{\beta -1}{\pi {\alpha }_x{\alpha }_y}{M}_1(x,y),\end{array}$$(4)

where the notation M₁ must be understood as the Moffat M_A with a multiplicative factor A = 1.

The full fitting method will be presented in Sect. 3, but we show here a preliminary result using the Moffat function to motivate our search for better functions. Indeed, as shown in Fig. 1, the Moffat function accurately fits the central peak of an actual PSF, but poorly describes the turbulent halo. Since this halo may contain an important proportion of the PSF energy, depending on the quality of the AO correction, it is necessary to model it accurately. Adding a constant background to the model artificially improves the fitting (lower residuals). However, this method is not suitable since it poorly describes the halo and mistaking the halo for a background will yield the non-physical result of a PSF with an infinite integral on an unlimited field of view. The modulation transfer functions (MTF, bottom plot in Fig. 1), which is the modulus of the PSF Fourier transform, also shows that the Moffat does not match well the very low frequencies (halo) and does not model the telescope cutoff frequency. Similarly, none of the static aberrations of the telescope are taken into account. A more physical PSF model than a Moffat is thus required.

Fig. 1. Fitting of a SPHERE/ZIMPOL PSF (blue) using a Moffat model with background (green) and without background (dashed green). Top: PSF, bottom: MTF. The inset plot is a zoom on the low spatial frequencies.

2.2. Image formation theory

Our PSF model is based on equations of image formation from the phase PSD to the focal plane. Indeed Roddier (1981) has shown that the Fourier transform of the PSF, the optical transfer function (OTF), can be written as the product of the telescope aberrations OTF and the atmospheric turbulent OTF,

$$\begin{array}{c}\hfill \stackrel{\sim }{h}(\mathit{\rho }/\lambda )={\stackrel{\sim }{h}}_{\mathrm{T}}(\mathit{\rho }/\lambda )·{\stackrel{\sim }{h}}_{\mathrm{A}}(\mathit{\rho }/\lambda ),\end{array}$$(5)

where λ is the observation wavelength, $\stackrel{\sim }{h}$ the total OTF, ${\stackrel{\sim }{h}}_{\mathrm{T}}$ the telescope OTF, and ${\stackrel{\sim }{h}}_{\mathrm{A}}$ the atmospheric OTF. This OTF splitting equation is valid under the hypothesis of a spatially stationary phase. This is the case for a purely turbulent phase, and a good approximation for an AO-corrected phase (Conan 1994). To establish this result, Roddier also used the fact that the phase distribution follows a Gaussian process, as the sum of a large number of independent turbulent layers. The telescope OTF is simply given by the autocorrelation of the pupil transmission function, whereas the atmospheric OTF is written as

$$\begin{array}{c}\hfill {\stackrel{\sim }{h}}_{\mathrm{A}}(\mathit{\rho }/\lambda )={\mathrm{e}}^{-B_{\varphi }(\mathbf{0})}{\mathrm{e}}^{B_{\varphi }(\mathit{\rho })},\end{array}$$(6)

with B_ϕ the phase autocorrelation function defined as

$$\begin{array}{c}\hfill B_{\varphi }(\mathit{\rho })={⟨{⟨\varphi (\mathit{r},t)\varphi (\mathit{r}+\mathit{\rho },t)⟩}_t⟩}_r.\end{array}$$(7)

The Wiener–Khintchine theorem states that the PSD is the Fourier transform of the autocorrelation as

$$\begin{array}{c}\hfill W_{\varphi }(\mathit{f})=\mathcal{F}\{B_{\varphi }(\mathit{\rho })\},\end{array}$$(8)

where W_ϕ denotes the phase PSD and ℱ the Fourier operator; f and ρ are the Fourier conjugated variables. If we call u the angular variable conjugated to ρ/λ, the PSF is written as

$$\begin{array}{c}\hfill h(\mathit{u})={\mathcal{F}}^{-1}\left\{{\stackrel{\sim }{h}}_{\mathrm{T}}(\mathit{\rho }/\lambda ){\mathrm{e}}^{-B_{\varphi }(\mathbf{0})}{\mathrm{e}}^{{\mathcal{F}}^{-1}\{W_{\varphi }(\mathit{f})\}}\right\}.\end{array}$$(9)

We note that B_ϕ(0) is the residual phase variance and is equal to the integral of W_ϕ on the whole frequency plane. This equation shows that only knowledge of the pupil and the static aberrations of the term ${\stackrel{\sim }{h}}_{\mathrm{T}}$, and the phase PSD W_ϕ are necessary for the description of the long-exposure PSF. Diffraction effects – such as finite aperture, central obstruction, and spiders – only depend on the pupil geometry and are known. Static aberrations are second-order effects that can be either neglected (as we show in Sects. 3.3 and 3.4), or measured (N’Diaye et al. 2013) and then included in the ${\stackrel{\sim }{h}}_{\mathrm{T}}$ term for more accuracy. The term ${\stackrel{\sim }{h}}_{\mathrm{T}}$ being fully determined, now we only have to parameterize the residual phase PSD to model the PSF.

2.3. Parameterization of the phase PSD

Actuators controlling the deformable mirror are separated by a pitch that sets the maximal spatial frequency of the phase that can be corrected by the AO system. This is called the AO spatial cutoff frequency, defined by f_AO ≃ N_act/2D, where N_act is the linear number of actuators and D the pupil diameter. This technical limitation induces the peculiar shape of the AO residual phase PSD (and in fine a peculiar shape of the PSF). The residual phase PSD is thus separated into two distinct areas:

• AO-corrected frequencies f ≤ f_AO,

• AO-uncorrected frequencies f >  f_AO.

The uncorrected area is not affected by the AO system, and the phase PSD consequently follows the Kolmogorov law,

$$\begin{array}{c}\hfill W_{\varphi ,\text{Kolmo}}(\mathit{f})=0.023r_0^{-5/3}{f}^{-11/3}\text{, for}f>f_{}\text{AO},\end{array}$$(10)

where r₀ is the Fried parameter scaling the strength of the turbulence. The halo is thus set by the knowledge of only this r₀ parameter.

Regarding the AO-corrected area, it is difficult to parameterize the phase residual PSD since it depends on the turbulence, the magnitude of the object, the AO loop delay, and the wavefront reconstruction algorithm. Our objective is not to build a full reconstruction of the phase PSD, but to only get a model that can match it. Racine et al. (1999) and Jolissaint & Veran (2002) have shown that in extreme AO correction (small residual phase) the shape of the PSF is exactly the shape of the PSD. For partial AO correction, the shapes of the PSF and PSD are not exactly identical, but are still similar (Fétick et al. 2018). Moreover Rigaut et al. (1998) have shown that the AO residual PSD is the sum of decreasing power laws of the spatial frequency. A Moffat function used in the PSD domain would already describe two regimes due to its shape, one regime for f ≤ α and one regime for α <  f <  f_AO. Adding a constant under the Moffat allows us to describe a third regime near the AO cutoff frequency at f ≲ f_AO that is roughly similar to the shape of the aliasing PSD discussed by Rigaut et al. (1998). All the above pieces of information suggest the possibility of using the Moffat function for a parsimonious parameterization of the AO-corrected PSD, rather than using it to directly parameterize the PSF in the focal plane. The full PSD model is written as

$$\begin{array}{c}\hfill W_{\varphi }(\mathit{f})={[\frac{\beta -1}{\pi {\alpha }_x{\alpha }_y}\frac{M_{\mathrm{A}}(f_x,f_y)}{1-{(1+\frac{f_{}{\text{AO}}^2}{{\alpha }_x{\alpha }_y})}^{1-\beta }}+C]}_{f\le f_{}\text{AO}}+{\left[W_{\varphi ,\text{Kolmo}}(\mathit{f})\right]}_{f>f_{}\text{AO}},\end{array}$$(11)

where the Moffat normalization factor ensures a unit integral of the Moffat on the area f ≤ f_AO (see Appendix A). Constant C is an AO-corrected phase PSD background. It is useful to model the residual AO PSD near the AO cutoff, where the Moffat function is close to zero. Thus the AO residual phase variance on the circular domain below the AO cutoff frequency is directly

$$\begin{array}{c}\hfill {\sigma }_{}{\text{AO}}^2=A+C\pi f_{}{\text{AO}}^2.\end{array}$$(12)

Parameter A, added to the C background contribution, has the physical meaning of being the residual variance. An example of our PSD model is given Fig. 2. We do not impose continuity at the AO cutoff frequency, so the PSD might be locally discontinuous. Indeed the transition area f ≃ f_AO between corrected and uncorrected frequencies can lead to strong PSD gradients, which are modelled by an eventual PSD discontinuity.

Fig. 2. Three components of the PSD model: the Moffat (blue) and the constant contribution (orange) below the AO cutoff frequency, and the Kolmogorov spectrum (green) above the AO cutoff frequency. Discontinuity has been exaggerated by reducing C to show this degree of freedom in our model. Plotting is in logarithmic–logarithmic scale.

Our PSF model based on the PSD is made of the following set of seven parameters: 𝒮 = {α_x, α_y, β, θ_R, C, r₀, A} in the asymmetric case. This reduces to five parameters in the symmetric case (setting α_y = α_x and θ_R = 0). Even though symmetric PSFs are sufficient in the majority of cases, asymmetries make it possible to consider PSFs elongated due to strong wind effects or anisoplanetism. Once the parameters 𝒮 are set, the PSF is then computed from the AO PSD and static aberrations using Eq. (9).

For the reader interested in deriving the Strehl ratio from our model, we have to compute the integral of the Kolmogorov spectrum above the AO cutoff frequency,

$$\begin{array}{cc}\hfill {\sigma }_{}^2\text{halo}& =0.023r_0^{-5/3}2\pi {\displaystyle {\int }_{f_{}\text{AO}}^{+\infty }{f}^{-11/3}f\mathrm{d}f}\hfill \\ \hfill & =0.023\frac{6\pi }{5}{(r_0·f_{}\text{AO})}^{-5/3}.\hfill \end{array}$$(13)

The Strehl ratio consequently is written as

$$\begin{array}{cc}\hfill {\text{S}}_{}\text{R}& =exp[-({\sigma }_{}^2\text{AO}+{\sigma }_{}^2\text{halo})]\hfill \\ \hfill & =exp[-A-C\pi f_{}^2\text{AO}-0.023\frac{6\pi }{5}{(r_0·f_{}\text{AO})}^{-5/3}].\hfill \end{array}$$(14)

3. Validation

3.1. PSF fitting method

In this section we deal with images of PSFs (the data) that may come from numerical simulations or observations of stars on VLT instruments. The fitting method consists in finding the PSF parameters so that the model PSF minimizes the square distance to the data PSF:

$$\begin{array}{c}\hfill \mathcal{L}(\mathcal{S},\gamma ,\zeta ,{\delta }_x,{\delta }_y)={\displaystyle \sum _{i,j}{w}_{i,j}{[\gamma ·h_{i,j}(\mathcal{S},{\delta }_x,{\delta }_y)+\zeta -d_{i,j}]}^2,}\end{array}$$(15)

where h_i, j is the discretized model of PSF on the pixels (i, j), 𝒮 its set of parameters, and d_i, j is the data PSF. Since the PSF model (given by Eqs. (9) and 11) has a unit flux, it is scaled by γ to match the flux of the data PSF, and ζ accounts for a possible background. The shifts δ_x and δ_y centre the PSF with sub-pixel precision on the data (by multiplication of the OTF with the correct phasor). The weighting factor w_i, j is the inverse of the noise variance, which takes into account the photon noise and the detector read-out noise. As noted by Mugnier et al. (2004), for high fluxes (typically greater than ten photons per pixel), the Poisson photon noise becomes nearly Gaussian and the weighting factor is written as

$$\begin{array}{c}\hfill w_{i,j}=\frac{1}{max\{d_{i,j},0\}+{\sigma }_{}{\text{RON}}^2}·\end{array}$$(16)

In this case, our approach can be seen from a statistical point of view as maximizing the likelihood of the data d_i, j corrupted by photon and read-out noise. We thus minimize ℒ, which is the neg-logarithm of the likelihood for a Gaussian process.

Let us now note that the minimum of ℒ has an analytic solution for γ and ζ (see Appendix B for a full demonstration). We actually do not need to numerically minimize over these two parameters, the least-square criterion only relies on the PSF intrinsic parameters 𝒮 and the position parameters (δ_x, δ_y) as

$$\begin{array}{c}\hfill \mathcal{L}{}^{\prime }(\mathcal{S},{\delta }_x,{\delta }_y)=\mathcal{L}(\mathcal{S},\widehat{\gamma },\widehat{\zeta },{\delta }_x,{\delta }_y),\end{array}$$(17)

where $\widehat{\gamma }$ and $\widehat{\zeta }$ are the analytic solutions for the flux and the background, respectively. At each iteration of the minimization process, the minimizer evaluates our ℒ′ criterion with a new set of parameters (𝒮, δ_x, δ_y). The current PSF estimate h(𝒮, δ_x, δ_y) is computed, then the analytic solutions $\widehat{\gamma }$ and $\widehat{\zeta }$ are computed. The quantity $\widehat{\gamma }·h(\mathcal{S},{\delta }_x,{\delta }_{\mathit{y}})+\widehat{\zeta }$ is used to compute the residuals with the data d. Residuals are then provided to the minimizer to estimate a new set of parameters (𝒮, δ_x, δ_y).

3.2. OOMAO end-to-end simulations

The Object-Oriented Matlab Adaptive Optics (OOMAO) toolbox, presented by Conan & Correia (2014), provides end-to-end simulations. For each time step, OOMAO generates a turbulent wavefront with a Von-Kármán spectrum defined as

$$\begin{array}{c}\hfill W_{\varphi ,\text{VK}}(\mathit{f})=0.023r_0^{-5/3}{[{\left(\frac{1}{L_0}\right)}^2+f^2]}^{-11/6},\end{array}$$(18)

where we have chosen the outer scale L₀  =  30 m. Since 1/L₀ ≪ f_AO, the Von-Kármán spectrum is consistent with our PSF model using the Kolmogorov spectrum above the AO cutoff frequency. OOMAO then propagates the wavefront through the telescope, simulates the wavefront sensor measurement, performs the wavefront reconstruction, and simulates the mirror wavefront deformation. For each time step, we get a short exposure PSF. Integration over time allows us to retrieve the long-exposure PSF. It is important to notice that the method to compute the PSF is then very different from our model, which directly uses the residual phase PSD in Eq. (9). We used OOMAO to generate a set of PSFs, corresponding to different wavelengths from 0.5 um to 2.18 um, and seven different values of r₀ from 7.5 cm to 25.0 cm. All r₀ are given at 500 nm. We translate them from the observation wavelength to the reference wavelength of 500 nm using the theoretical spectral dependency r_0, λ1/r_0, λ2  =  (λ₁/λ₂)^6/5. For all the PSFs, the telescope parameters are kept unchanged D = 8 m, N_act = 32, sampling at Shannon–Nyquist for all wavelengths. The phase screen consists in one frozen flow (Taylor’s hypothesis) turbulent layer translating at v = 10 m s⁻¹. Using these parameters, we generated PSFs corresponding to exposure times of 0.1, 1, 10, and 100 s. For 0.1 and 1 s, the random OOMAO phase screen did not converge towards a stable state, leading to a strong bias in the r₀ estimation. For a 10 s exposure PSF, the random fluctuations of the phase are correctly averaged. This is confirmed by the 100 s exposure PSF, which gives the same r₀ estimation as the 10 s case. Since a 100 s exposure is computationally demanding and does not significantly improve the results, we performed all our tests on the 10 s exposure time (Table 1).

Each PSF is then fitted using the ℒ′(𝒮, δ_x, δ_y) criterion given to an optimizer (e.g. Levenberg–Marquardt, trust region, or Markov chain Monte Carlo). We used the Trust Region Reflective algorithm, called “trf”, from the Python/SciPy (Jones et al. 2001) library. This algorithm is gradient based, said to be robust, and allows bounds on the parameters. The robustness of this algorithm was verified for our applications of it to PSF fitting, even though the convexity of the problem is not demonstrated. So far, we have not encountered any local minimum and residuals are always small. For all PSFs, the same initial conditions are provided to the fitting algorithm (see Table 2), in particular we used the same value of r₀ = 18 cm. Using the same initial parameters {𝒮, δ_x, δ_y}_init for all fits ensures that our model is suited for minimization procedures and that convergence is ensured even if starting far from the true values. Fitting results are presented on Fig. 3. Our model fits well the OOMAO-generated PSF on both the corrected and the uncorrected area, residuals being on average one to two decades below the PSF. Let us define the relative error between fitted PSF and data PSF as

$$\begin{array}{c}\hfill {ϵ}_h=\frac{\sqrt{\sum _{i,j}{[\gamma ·h_{i,j}(\mathcal{S},{\delta }_x,{\delta }_y)+\zeta -d_{i,j}]}^2}}{{\sum }_{i,j}{d}_{i,j}}·\end{array}$$(19)

Fig. 3. OOMAO PSF fitting with our model. Left: circular average for PSFs (given in photons). The vertical grey line corresponds to the AO cutoff radius. Right: corresponding circular average for OTFs (normalized to unity at the null frequency). From top to bottom, three wavelengths are scanned from 500 nm to 1220 nm. Colours correspond to three values of the OOMAO required r0. Solid curves: OOMAO. Dashed: fitting. Dotted: residuals. All PSFs, for all wavelengths, are sampled at Shannon–Nyquist.

This error is the L2 norm of the differences between fitting and data, relative to the flux. Considering all the OOMAO fitting, we find an average relative error ϵ_h = 6.4 × 10⁻³ with a standard deviation of 2.2 × 10⁻³. For comparison, fitting with a Moffat model gives an average relative error of ϵ_h = 3.4 × 10⁻² with a standard deviation of 1.6 × 10⁻². Using our model thus increases the fitting accuracy by a factor of approximately 5 with respect to the former Moffat model of Sect. 2.1. Regarding the OTFs, the fit is also accurate on the whole frequency range, from low frequencies (mainly the halo) to high frequencies (PSF peak and telescope cutoff). Our model slightly over-estimates frequencies just below the telescope cutoff frequency but this has never been an issue in our applications, such as deconvolution, and is still much better than the Moffat OTF (which has no telescope cutoff frequency and give a poor estimation of the low frequencies).

Regarding the flux, we consider the relative error between the flux $\widehat{\gamma }$ analytically estimated by our model fitting method, and the OOMAO flux that is directly the sum of the data on all the pixels:

$$\begin{array}{c}\hfill {ϵ}_{\gamma }=\frac{\widehat{\gamma }-{\sum }_{i,j}{d}_{i,j}}{{\sum }_{i,j}{d}_{i,j}}·\end{array}$$(20)

On all our OOMAO simulations, we find an average relative error of −1.96%, indicating a small underestimation of the flux with our fitting method. The standard deviation of this relative error is 1.11%, and the range of variation is [ − 3.43%,2.77%].

3.2.1. Fried parameter r₀ estimation

As shown in Fig. 4, our r₀ estimation is consistent with the OOMAO value of r₀. We find the best linear fit to be

$$\begin{array}{c}\hfill r_{0,\text{FIT}}=1.038r_{0,\text{OOMAO}}-0.132,\end{array}$$(21)

Fig. 4. Fried parameter r0 estimated by fitting versus the r0 used in OOMAO to generate the PSF. All r0 are given at 500 nm. Here are shown results on 84 different PSFs, corresponding to seven values of r0 and 12 different wavelengths. The line is the linear fit between our r0 estimation and OOMAO r0. Crosses show residuals |r0, FIT − r0, OOMAO|. A log–log scale is used to show on the same graph both data and residuals.

where values are given in centimetres. The Pearson correlation coefficient is 𝒞_Pearson = Cov(r_0, FIT, r_0, OOMAO)/$\sqrt{\text{Var}(r_{0,\text{FIT}})·\text{Var}(r_{0,\text{OOMAO}})}=0.99992$. This result fully confirms our r₀ estimation with respect to OOMAO simulations with sub-centimetre precision.

3.2.2. AO residual variance ${\sigma }_{\text{AO}}^2$ estimation

Theoretically our model should also be able to retrieve the residual variance ${\sigma }_{\text{AO}}^2$ on the corrected area and follow a λ⁻² power law. Figure 5 shows the fitting estimation of this variance versus the wavelength. This data is then fitted with curves of equation aλ⁻². Except the two outliers for minimal r₀ = 7.5 cm at low wavelength (λ ≃ 500 nm), the λ⁻² power law is a good estimation of the ${\sigma }_{\text{AO}}^2$ evolution. This result gives confidence in the estimated parameter. Data from the real time computer (RTC) could be used in the future to provide the ${\sigma }_{\text{AO}}^2$ parameter for PSF estimation. The λ⁻² spectral dependence is also an asset to shift the PSF from one wavelength to another.

Fig. 5. Estimation of the ${\sigma }_{\text{AO}}^2$ from PSF fitting versus the wavelength (dots). Colours correspond to the seven different values of OOMAO r0. Curves of parametric equations σ2 = aλ−2 are fitted on the data.

3.2.3. Constant C estimation

The C term in Eq. 11 accounts for multiple sources of residual PSD, including wavefront aliasing and other AO residual errors. Since this constant is dominated by the Moffat PSD in the core, it becomes more important near the AO cutoff, where the aliasing dominates. Since the aliasing scales in $r_0^{-5/3}$ (Rigaut et al. 1998), we look for similar r₀ dependencies for the PSF constant C. Figure 6 shows a clear decrease of C with r₀. Fitting the estimated C with a $r_0^{-5/3}$ power law shows a good match, with small residuals for nearly all r₀ values. The power law is not exactly −5/3 ≃ −1.67 but is closer to −1.46 for this OOMAO case. However, one can still think about normalizing the constant in Eq. 11 by

$$\begin{array}{c}\hfill C=C^{\prime }{r}_0^{-5/3}\end{array}$$(22)

Fig. 6. Estimation of the PSF constant C versus the r0 given at the observed wavelength (dots). A $r_0^{-5/3}$ fitting equation (solid line) is applied on the data. Residuals between each data point and the $r_0^{-5/3}$ power law are also shown (crosses).

and perform fitting over C′ instead of C. This would reduce the variation range of this parameter. Reducing bounds or standard deviation of a parameter is an asset for constraining the model and improving minimization processes.

3.3. High performance imager ZIMPOL

The Spectro Polarimetric High-contrast Exoplanet REsearch (SPHERE) instrument (Beuzit et al. 2008, 2019) of the VLT includes the powerful SPHERE Adaptive optics for eXoplanets Observation (SAXO) system described in Fusco et al. (2014) and Sauvage et al. (2010). The AO real time computer is built on the ESO system called SPARTA (Fedrigo et al. 2006), which stands for Standard Platform for Adaptive optics Real Time Applications. In particular, for each observation, SPARTA is able to give an estimate of r₀ from the mirror voltages and wavefront sensor slopes.

The Zurich IMaging POLarimeter (ZIMPOL) instrument (Schmid et al. 2018) is mounted at the focal plane of SPHERE. ZIMPOL is also used as a very efficient imager at visible wavelengths. One of its applications in non-coronagraphic mode is the observation of asteroids (Vernazza et al. 2018; Viikinkoski et al. 2018; Fétick et al. 2019) as part of an ESO Large Program (ID 199.C-0074, PI P. Vernazza). PSFs from stars were observed with the ZIMPOL N_R filter (central wavelength 645.9 nm, width 56.7 nm) during the Large Program. When PSFs are saved together with the SPARTA telemetry, we are able to correlate r₀ given by our fitting and r₀ given by SPARTA. In our sample, 28 PSFs were saved along with the SPARTA telemetry. These PSFs were obtained during different nights, on stars of different magnitudes, with various seeing conditions. We fitted these 28 PSFs with our PSF model. Figure 7 shows three of the 28 fittings, for the smallest r₀ of the sample, the median r₀, and the largest r₀, respectively. Our fitted PSFs match the shape of the core and halo. The average of relative error defined in Eq. (19) is ϵ_h = 2.5 × 10⁻³, with a standard deviation of 1.2 × 10⁻³. We only consider diffraction due to the telescope 8 m aperture and its central obstruction in the static OTF ${\stackrel{\sim }{h}}_{\mathrm{T}}$ (see Eq. (9)). Consequently non-circularly-symmetric effects, such as the spiders or static aberrations visible on Fig. 7, are not modelled. Even if the spiders could have been included in our model, we have deliberately chosen to ignore them in the ${\stackrel{\sim }{h}}_{\mathrm{T}}$ term since they are negligible in comparison to the other dominant effects (AO residual core, turbulent halo, 8 m aperture diffraction). When performing PSF fitting, the contribution of these effects not taken into account in ${\stackrel{\sim }{h}}_{\mathrm{T}}$ might bias the atmospheric term ${\stackrel{\sim }{h}}_{\mathrm{A}}$ during fitting procedure and slightly offset the estimation of the 𝒮 parameters. Moreover these ZIMPOL images are field stabilized, meaning rotating spiders, which are harder to model. Pupil stabilized images would make the description of the spider diffraction effect easier.

Fig. 7. Three ZIMPOL PSFs (top), model fittings (middle), residuals (bottom). Left: minimal r0 of the sample. Middle: median r0. Right: maximal r0. The main differences are due to some static aberrations not taken into account in our model (only the pupil and its central obstruction are taken into account). The hyperbolic arcsine of the intensity is shown to enhance details. The same intensity scale is used per column (data, model, residuals), but differs between columns.

The top graph on Fig. 8 shows that the values estimated by SPARTA (median = 22 cm) are greater than the values given by fitting (median = 13 cm). The best linear fit between r₀ estimated by fitting and SPARTA is

$$\begin{array}{c}\hfill r_{0,\text{SPARTA}}=3.41r_{0,\text{FIT}}-16.82.\end{array}$$(23)

Fig. 8. Zenital r0 at 500 nm estimated on MUSE (filled circles, 40 data points) and SPHERE/ZIMPOL (empty circles, 27 data points) using three methods: PSF fitting, SPARTA, and MASS/DIMM. Top: r0, SPARTA versus r0, FIT. A different linear tendency is found for MUSE (plain line) and SPHERE/ZIMPOL (dashed line). Shaded areas show the standard deviation between data points and the best linear fit. Middle: r0, SPARTA and r0, FIT versus r0, MASS/DIMM. Two linear tendencies are identified for SPARTA, and only one for PSF fitting. Bottom: histograms of seeing estimated with the three methods on both instruments. Dashed vertical lines show median values of 0.46″ (SPARTA), 0.69″ (MASS/DIMM), and 0.83″ (PSF fitting).

The Pearson correlation coefficient between the two series r_0, FIT and r_0, SPARTA is 𝒞_Pearson = Cov(r_0, FIT, r_0, SPARTA)/$\sqrt{\text{Var}(r_{0,\text{FIT}})·\text{Var}(r_{0,\text{SPARTA}})}=0.97$. From this data it appears that the estimates of r_0, FIT and r_0, SPARTA are not identical, however they show a strong correlation. We further investigated the difference between the SPARTA and the fitting estimates thanks to the ESO atmospherical monitoring using the Multi-Aperture Scintillation Sensor (MASS) combined with the Differential Image Motion Monitor (DIMM). Since the MASS/DIMM instrument is located apart from the telescopes, it does not see exactly the same turbulent volume as the telescopes and does not suffer the same dome effect. There might be some uncertainties between MASS/DIMM r₀ estimations and telescope r₀ estimations (PSF fitting or RTC) due to the spatial evolution of the turbulence. Nevertheless this instrument is a valuable indicator of the Paranal atmospheric statistics. For each PSF observation we retrieved the associated MASS/DIMM seeing estimation within a delay of ±3 min (see Fig. 8, middle and bottom graphics). The median seeing estimated by the MASS/DIMM is 0.69″, to be compared with a median seeing of 0.46″ for SPARTA and 0.83″ for PSF fitting. The over-estimation of the SPARTA r₀ with respect to the MASS/DIMM r₀ has been already discussed by Milli et al. (2017). The exact origin of the difference between these three estimations has not be found. Nevertheless we note that estimations with SPARTA are based on RTC measurements of the low spatial frequencies of the phase (sensitive to the Von-Kármán outer scale L₀), whereas our fitting method is based on the PSF halo corresponding to the high spatial frequencies. Our PSF fitting method might be sensitive to telescope internal wavefront errors if they have not been previously calibrated and taken into account in the PSF model.

Even if there is still an uncertainty on the true value of r₀, the strong correlation between SPARTA and fitting estimations is sufficient for many applications. Indeed it is still possible to get r_0, SPARTA from telemetry, use Eq. (23) to translate it into r_0, FIT, and get an estimate of the PSF halo. This method constrains the model for future PSF estimations without having access to the actual image of the PSF.

3.4. MUSE integral field spectrograph

The Multi-Unit Spectroscopic Explorer MUSE (Bacon et al. 2006, 2010) is an integral field spectrograph (IFS) working mainly in the visible, from ∼465 nm to ∼930 nm. MUSE is equipped with the Ground Atmospheric Layer Adaptive Optics for Spectroscopic Imaging (GALACSI) adaptive optics system (Ströbele et al. 2012) to improve its spatial resolution in two different modes, the so-called narrow-field mode (NFM) and wide-field mode (WFM), to correct different sizes of field of view. The AO facility uses four laser guide stars (LGS; Calia et al. 2014) to perform a tomographic reconstruction of the turbulent phase. A 589 nm dichroic is present in the optical path to avoid light contamination from the sodium AO lasers, so no scientific information is available around this wavelength. Let us also note that MUSE is undersampled (sampling is 25 mas in NFM) on the whole available spectrum. Our PSF model manages the undersampling issue by oversampling the PSD and the OTF to safely perform numerical computations. The given PSF is then spatially binned to retrieve the correct sampling. The shape of the PSF can be retrieved, but this method has lower precision on the parameters’ estimation inherent to undersampled data. These differences between the MUSE instrument and SPHERE/ZIMPOL allow us to test the versatility of our PSF model. Additionally the spectral resolution is an asset to validate our model at different wavelengths and to study the spectral evolution of our PSF parameters.

During the May–June 2018 commissioning phase, MUSE observed multiple targets in narrow-field mode. Among these targets, we have access to 40 PSFs observed on different stars, during different nights and at different seeing conditions. These selected PSFs have been spectrally binned into 92 bins of 5 nm each to increase the signal to noise ratio and reduce the number of fittings. Then fitting is performed independently, spectral bin by spectral bin, without any spectral information on the targets or the atmosphere. Figure 9 shows one MUSE datacube PSF fitting at three different wavelengths. The evolution of the AO correction radius is clearly visible in both the data and the model. As for ZIMPOL, we did not take into account static PSF (except the occulted pupil diffraction), which is the main visible difference between data and model. Fainter stars visible in the field did not affect the fitting and appear clearly in the residuals. For the 40 datacubes PSF, the relative error is ϵ_h = 3.3 × 10⁻³. This result is similar to the previous ones on OOMAO and SPHERE/ZIMPOL. Secondary stars in the field also count in the residual error computation.

Fig. 9. MUSE PSFs (top), fittings (middle), and residuals (bottom) of the same star at three different wavelengths over the 92 spectral bins actually fitted. The hyperbolic arcsine of the intensity is shown to enhance details.

The evolution of r₀ with the wavelength is shown on Fig. 10 for one datacube PSF. A least square fitting between our data points and the theoretical λ^6/5 evolution of the r₀ gives a spectral averaged estimation of $\overline{r_0}=13.3$ cm at 500 nm. Our fitted r₀ matches well the theory, with a standard deviation of 0.3 cm.

Fig. 10. Estimation of the r0 by fitting of the PSF MUSE at 92 different wavelengths (blue dots), and comparison with a theoretical law in λ6/5 (orange line). The best match between data and theory is achieved for r0 = 13.3 cm. The grey area corresponds to missing wavelengths due to the sodium notch filter.

So far each spectral bin is fitted independently, however the spectral deterministic trend we recover is an asset for PSF determination. It makes possible the fitting of the whole datacube with only one r₀ parameter. The statistical contrast – ratio of the number of measurements over the number of unknowns – would be increased. It would improve fitting robustness, especially for faint stars where the halo is strongly affected by noise.

The results of fitting on the 40 datacubes give statistical information on the r₀ estimation. As for the SPHERE/ZIMPOL case, we have access to SPARTA and MASS/DIMM data to correlate with our fitting estimations (see Fig. 8). The Pearson correlation coefficient between PSF fitting and SPARTA is 𝒞_Pearson = 0.96, which is similar to the SPHERE/ZIMPOL case. However, we get the linear relationship

$$\begin{array}{c}\hfill r_{0,\text{SPARTA}}=1.96r_{0,\text{FIT}}-5.31,\end{array}$$(24)

which is different from the SPHERE/ZIMPOL (Eq. (23)). The exact origin of this different trend is unknown, nevertheless let us note that the actual implementation of the phase PSD estimation is slightly different for the SPHERE and for the MUSE instruments. For SPHERE the Von-Kármán outer scale L₀ is set to 25 m, whereas for MUSE the L₀ is estimated jointly with r₀. This SPARTA double trend is corroborated by MASS/DIMM information (Fig. 8, middle graph). On the other hand, r₀ estimation using our PSF model gives similar results on both SPHERE/ZIMPOL and MUSE instruments. This confirms the robustness of our PSF fitting method.

4. Conclusions

In this article we developed a parametric model of long-exposure AO-corrected PSF. The particularity of this model is to parameterize the phase PSD using a Moffat core and a turbulent Kolmogorov halo. This model also incorporates prior knowledge of the telescope, such as the optical cutoff frequency, the obstruction and spider shapes, and even the static aberrations if they are calibrated, for example by phase diversity (Mugnier et al. 2008). This model only requires five parameters for circularly symmetrical PSFs, and seven for asymmetrical ones. The sparsity of this PSF model makes it suitable for numerical computation, such as minimization algorithms or least-square fits. Tests on both simulated and real data validated the appropriateness of our model.

One substantial advantage of our model over focal plane models is to use physical parameters such as the Fried parameter r₀ and the residual AO variance ${\sigma }_{\text{AO}}^2$. Since these parameters are physical, their values in our PSF model can be correlated to external measurements. Tests on both OOMAO simulations and on-sky data (from the SPHERE/ZIMPOL and MUSE instruments) confirmed the physical meaning of the r₀ parameter used in our PSF. The ultimate goal would be to only use physical parameters in the PSF description.

Our model has already shown usability for different seeings and different instruments, with different AO-correction quality. This shows the robustness and versatility of the model. We also plan to use it to parameterize the PSF for the future instruments on bigger telescopes such as the Extremely Large Telescope (ELT).

Finally, the small number of parameters makes this model suited for image post-processing techniques such as deconvolution of long-exposure images. Deconvolution using parametric PSFs has already been demonstrated by Drummond (1998) and Fétick et al. (2019). We plan to develop a myopic deconvolution algorithm estimating both the observed object and the PSF parameters in a marginal approach similar to Blanco & Mugnier (2011).

Acknowledgments

This work was supported by the French Direction Générale de l’Armement (DGA) and Aix-Marseille Université (AMU). This work was supported by the Action Spécifique Haute Résolution Angulaire (ASHRA) of CNRS/INSU co-funded by CNES. This project has received funding from the European Union’s Horizon 2020 research and innovation program under grant agreement No. 730890. This material reflects only the authors views and the Commission is not liable for any use that may be made of the information contained therein. This study has been partly funded by the French Aerospace Lab (ONERA) in the frame of the VASCO Research Project. The authors are thankful to Pierre Vernazza for providing data related to his ESO Large Program ID 199.C-0074.

References

Appendix A: Integral of a truncated Moffat

A Moffat function as given in Eq. (2) shows elliptical contours (see Fig. A.1). In this appendix we compute the integral of the Moffat function inside one of these elliptical contours, called 𝔼(R_x, R_y), of semi-major axis R_x and semi-minor axis R_y = (α_y/α_x)R_x. The integral to calculate is

$$\begin{array}{c}\hfill I(R_x,R_y)={\displaystyle \int {\int }_{\mathbb{E}(R_x,R_y)}{M}_{\mathrm{A}}(x,y)\mathrm{d}x\mathrm{d}y.}\end{array}$$(A.1)

Fig. A.1. Visualization of an elongated Moffat (colour map), its elliptical level curves (white), and the ellipse 𝔼(Rx, Ry) inside which the integral is computed. Dashed circles of radius Rx and Ry help to visualize the error done when computing the integral over the ellipse instead of a circle. A very elongated Moffat is shown here (αx = 1.5αy).

Let us perform the change of variables

$$\begin{array}{c}\hfill \varphi :\{\begin{array}{cc}{\mathbb{R}}_{+}^{\ast }\times [-\pi ,\pi [& ⟶\mathbb{R}\times \mathbb{R}\backslash (0,0)\hfill \\ (r,\theta )& ⟶({\alpha }_{x}rcos\theta ,{\alpha }_{y}rsin\theta )\hfill \end{array}.\end{array}$$(A.2)

The determinant of the Jacobian is

$$\begin{array}{c}\hfill det{J}_{\varphi }=\left|\begin{array}{cc}{\alpha }_{x}cos\theta & -r{\alpha }_{x}sin\theta \\ {\alpha }_{y}sin\theta & r{\alpha }_{y}cos\theta \end{array}\right|={\alpha }_x{\alpha }_{y}r.\end{array}$$(A.3)

The integral of the Moffat in the ellipse is rewritten as

$$\begin{array}{cc}\hfill I(R_x,R_y)& ={\displaystyle {\int }_{-\pi }^{\pi }{\int }_0^{R_x/{\alpha }_x}|det{J}_{\varphi }|M_{\mathrm{A}}(\varphi (r,\theta ))\mathrm{dr}\mathrm{d}\theta }\hfill \\ \hfill & =2\pi A{\alpha }_x{\alpha }_y{\displaystyle {\int }_0^{R_x/{\alpha }_x}{[1+r^2]}^{-\beta }}\hfill \\ \hfill & =A\frac{\pi {\alpha }_x{\alpha }_y}{\beta -1}\{1-{[1+{(R_x/{\alpha }_x)}^2]}^{1-\beta }\},\hfill \end{array}$$(A.4)

where we assumed β ≠ 1.

In the nearly circular regime α_x ≃ α_y and R_x ≃ R_y so the integral over the ellipse 𝔼(R_x, R_y) is nearly equal to the energy in a disk of radius R_x. This assumption is made for our model of PSD to compute the residual variance below the AO cutoff frequency in Eq. (11).

Moreover it is possible to calculate the integral of the Moffat function over the whole plane. The assumption β >  1 is mandatory to get a bounded integral (finite energy). Under this assumption, letting R_x → +∞ and R_y → +∞, one finds

$$\begin{array}{c}\hfill I(+\infty ,+\infty )=A\frac{\pi {\alpha }_x{\alpha }_y}{\beta -1}·\end{array}$$(A.5)

Consequently, in order to get a Moffat of unit integral over the whole plane, one must choose the Moffat amplitude factor as

$$\begin{array}{c}\hfill A=\frac{\beta -1}{\pi {\alpha }_x{\alpha }_y}·\end{array}$$(A.6)

Appendix B: Analytic solution for the flux and background

The minimum of ℒ (given in Eq. (15)) has an analytic solution for γ and ζ since nulling the partial derivative of ℒ towards these parameters gives

$$\begin{array}{cc}\hfill \frac{\partial \mathcal{L}}{\partial \gamma }=0& ⟺{\displaystyle \sum _{i,j}{w}_{i,j}{h}_{i,j}[\gamma ·h_{i,j}+\zeta -d_{i,j}]=0}\hfill \\ \hfill & ⟺\gamma {\displaystyle \sum _{i,j}{w}_{i,j}{h}_{i,j}^2+\zeta \sum _{i,j}{w}_{i,j}{h}_{i,j}=\sum _{i,j}{w}_{i,j}{h}_{i,j}{d}_{i,j}}\hfill \end{array}$$(B.1)

and

$$\frac{\partial \mathcal{L}}{\partial \zeta }=0⟺{\displaystyle \sum _{i,j}{w}_{i,j}\left[\gamma \cdot h_{i,j}+\zeta -d_{i,j}\right]=0}$$(B.2)

$$⟺\gamma {\displaystyle \sum _{i,j}{w}_{i,j}{h}_{i,j}+\zeta {\displaystyle \sum _{i,j}{w}_{i,j}}={\displaystyle \sum _{i,j}{w}_{i,j}}{d}_{i,j}.}$$(B.3)

These two equations are linear in γ and ζ. They can be written within the matrix formalism

$$\begin{array}{c}\hfill \mathcal{A}·\left(\begin{array}{c}\zeta \\ \gamma \\ \end{array}\right)={\displaystyle \sum _{i,j}{w}_{i,j}{d}_{i,j}\left(\begin{array}{c}1\\ h_{i,j}\\ \end{array}\right),}\end{array}$$(B.4)

where 𝒜 is the 2 × 2 matrix of the system defined as

$$\begin{array}{c}\hfill \mathcal{A}={\displaystyle \sum _{i,j}{w}_{i,j}\left(\begin{array}{cc}1& h_{i,j}\\ h_{i,j}& h_{i,j}^2\\ \end{array}\right).}\end{array}$$(B.5)

In order to invert 𝒜, we need to make sure that det(𝒜)≠0. The determinant of 𝒜 is

$$\begin{array}{cc}\hfill \text{det}(\mathcal{A})& =\left({\displaystyle \sum _{i,j}{w}_{i,j}}\right)\left({\displaystyle \sum _{i,j}{w}_{i,j}{h}_{i,j}^2}\right)-{\left({\displaystyle \sum _{i,j}{w}_{i,j}{h}_{i,j}}\right)}^2\hfill \\ \hfill & =\left({\displaystyle \sum _{i,j}{\sqrt{w_{i,j}}}^2}\right)\left({\displaystyle \sum _{i,j}{(\sqrt{w_{i,j}}{h}_{i,j})}^2}\right)-{\left({\displaystyle \sum _{i,j}(\sqrt{w_{i,j}})(\sqrt{w_{i,j}}{h}_{i,j})}\right)}^2.\hfill \end{array}$$(B.6)

We can now define the two vectors $\sqrt{W}=\{\sqrt{{\mathit{w}}_{i,j}}\}$ and $\sqrt{W}H=\{\sqrt{{\mathit{w}}_{i,j}}{h}_{i,j}\}$. Using these notations the determinant is rewritten as

$$\begin{array}{c}\hfill \text{det}(\mathcal{A})=\Vert \sqrt{W}{\Vert }^2\Vert \sqrt{W}{H\Vert }^2-{|⟨\sqrt{W},\sqrt{W}H⟩|}^2.\end{array}$$(B.7)

We recognize the Cauchy-Schwarz inequality. It follows that det(𝒜)≥0, with equality only if $\sqrt{W}$ and $\sqrt{W}H$ are colinear vectors. The colinearity is written as

$$\begin{array}{cc}\hfill \sqrt{W}//\sqrt{W}H& ⟺\exists k\in \mathrm{R}\mid \forall (i,j),k\sqrt{w_{i,j}}=\sqrt{w_{i,j}}{h}_{i,j}\hfill \\ \hfill & ⟺\exists k\in \mathrm{R}\mid \forall (i,j)\text{where}{w}_{i,j}\ne 0,k=h_{i,j}.\hfill \end{array}$$(B.8)

This states that the determinant of 𝒜 is null only if the PSF model h is constant on each point where w_i, j ≠ 0. Since our PSF is not constant on the domain w_i, j ≠ 0, we ensure that det(𝒜)≠0 and the analytic solution for γ and ζ is written as

$$\begin{array}{c}\hfill \left(\begin{array}{c}\zeta \\ \gamma \\ \end{array}\right)={\mathcal{A}}^{-1}{\displaystyle \sum _{i,j}{w}_{i,j}{d}_{i,j}\left(\begin{array}{c}1\\ h_{i,j}\\ \end{array}\right).}\end{array}$$(B.9)

All Tables

All Figures

	Fig. 1. Fitting of a SPHERE/ZIMPOL PSF (blue) using a Moffat model with background (green) and without background (dashed green). Top: PSF, bottom: MTF. The inset plot is a zoom on the low spatial frequencies.
In the text

	Fig. 2. Three components of the PSD model: the Moffat (blue) and the constant contribution (orange) below the AO cutoff frequency, and the Kolmogorov spectrum (green) above the AO cutoff frequency. Discontinuity has been exaggerated by reducing C to show this degree of freedom in our model. Plotting is in logarithmic–logarithmic scale.
In the text

Fig. 3. OOMAO PSF fitting with our model. Left: circular average for PSFs (given in photons). The vertical grey line corresponds to the AO cutoff radius. Right: corresponding circular average for OTFs (normalized to unity at the null frequency). From top to bottom, three wavelengths are scanned from 500 nm to 1220 nm. Colours correspond to three values of the OOMAO required r0. Solid curves: OOMAO. Dashed: fitting. Dotted: residuals. All PSFs, for all wavelengths, are sampled at Shannon–Nyquist.

In the text

Fig. 4. Fried parameter r0 estimated by fitting versus the r0 used in OOMAO to generate the PSF. All r0 are given at 500 nm. Here are shown results on 84 different PSFs, corresponding to seven values of r0 and 12 different wavelengths. The line is the linear fit between our r0 estimation and OOMAO r0. Crosses show residuals |r0, FIT − r0, OOMAO|. A log–log scale is used to show on the same graph both data and residuals.

In the text

	Fig. 5. Estimation of the ${\sigma }_{\text{AO}}^2$ from PSF fitting versus the wavelength (dots). Colours correspond to the seven different values of OOMAO r0. Curves of parametric equations σ2 = aλ−2 are fitted on the data.
In the text

	Fig. 6. Estimation of the PSF constant C versus the r0 given at the observed wavelength (dots). A $r_0^{-5/3}$ fitting equation (solid line) is applied on the data. Residuals between each data point and the $r_0^{-5/3}$ power law are also shown (crosses).
In the text

Fig. 7. Three ZIMPOL PSFs (top), model fittings (middle), residuals (bottom). Left: minimal r0 of the sample. Middle: median r0. Right: maximal r0. The main differences are due to some static aberrations not taken into account in our model (only the pupil and its central obstruction are taken into account). The hyperbolic arcsine of the intensity is shown to enhance details. The same intensity scale is used per column (data, model, residuals), but differs between columns.

In the text

Fig. 8. Zenital r0 at 500 nm estimated on MUSE (filled circles, 40 data points) and SPHERE/ZIMPOL (empty circles, 27 data points) using three methods: PSF fitting, SPARTA, and MASS/DIMM. Top: r0, SPARTA versus r0, FIT. A different linear tendency is found for MUSE (plain line) and SPHERE/ZIMPOL (dashed line). Shaded areas show the standard deviation between data points and the best linear fit. Middle: r0, SPARTA and r0, FIT versus r0, MASS/DIMM. Two linear tendencies are identified for SPARTA, and only one for PSF fitting. Bottom: histograms of seeing estimated with the three methods on both instruments. Dashed vertical lines show median values of 0.46″ (SPARTA), 0.69″ (MASS/DIMM), and 0.83″ (PSF fitting).

In the text

	Fig. 9. MUSE PSFs (top), fittings (middle), and residuals (bottom) of the same star at three different wavelengths over the 92 spectral bins actually fitted. The hyperbolic arcsine of the intensity is shown to enhance details.
In the text

	Fig. 10. Estimation of the r0 by fitting of the PSF MUSE at 92 different wavelengths (blue dots), and comparison with a theoretical law in λ6/5 (orange line). The best match between data and theory is achieved for r0 = 13.3 cm. The grey area corresponds to missing wavelengths due to the sodium notch filter.
In the text

	Fig. A.1. Visualization of an elongated Moffat (colour map), its elliptical level curves (white), and the ellipse 𝔼(Rx, Ry) inside which the integral is computed. Dashed circles of radius Rx and Ry help to visualize the error done when computing the integral over the ellipse instead of a circle. A very elongated Moffat is shown here (αx = 1.5αy).
In the text