Joint estimation of atmospheric and instrumental defects using a parsimonious point spread function model.On-sky validation using state of the art worldwide adaptive-optics assisted instruments

Modeling the optical point spread function (PSF) is particularly challenging for adaptive optics (AO)-assisted observations owing to the its complex shape and spatial variations. We aim to (i) exhaustively demonstrate the accuracy of a recent analytical model from comparison with a large sample of imaged PSFs, (ii) assess the conditions for which the model is optimal, and (iii) unleash the strength of this framework to enable the joint estimation of atmospheric parameters, AO performance and static aberrations. We gathered 4812 on-sky PSFs obtained from seven AO systems and used the same fitting algorithm to test the model on various AO PSFs and diagnose AO performance from the model outputs. Finally, we highlight how this framework enables the characterization of the so-called low wind effect on the SPHERE instrument and piston cophasing errors on the Keck II telescope. Over 4812 PSFs, the model reaches down to 4% of error on both the Strehl-ratio (SR) and full width at half maximum (FWHM). We particularly illustrate that the estimation of the Fried parameter, which is one of the model parameters, is consistent with known seeing statistics and follows expected trends in wavelength using the MUSE instrument ($\lambda^{6/5}$) and field (no variations) from GSAOI images with a standard deviation of 0.4cm. Finally, we show that we can retrieve a combination of differential piston, tip, and tilt modes introduced by the LWE that compares to ZELDA measurements, as well as segment piston errors from the Keck II telescope and particularly the stair mode that has already been revealed from previous studies. This model matches all types of AO PSFs at the level of 4% error and can be used for AO diagnosis, post-processing, and wavefront sensing purposes.


Introduction
Adaptive optics (AO, Roddier 1999) is a game changer in the quest for high-angular resolution, especially for groundbased astronomical observations that face the presence of wavefront aberrations introduced by the atmosphere (Roddier 1981). Thanks to AO, the point spread function (PSF) delivered by an optical instrument is much narrower, by a factor up to 50 on the full width at half maximum (FWHM), than the seeinglimited scenario (Roddier 1981). Still, some correction residuals persist and render the AO PSF shape complex to model. Consequently, standard parametric models that reliably reproduce the seeing-limited PSFs, such as a Moffat function (Trujillo et al. 2001;Moffat 1969), become inefficient at describing the AO PSF. Moreover, contrary to seeing-limited observations, AO-corrected images suffer from the anisoplanatism effect (Fried 1982) that strengthens the spatial variations of the PSF on top of instrument defects. Determining the AO PSF is necessary for two major reasons. Firstly, understanding and accurately modeling the PSF morphology is key to diagnosing AO performance. From the PSF, we can identify the major contributors to the AO residual error (Beltramo- Martin et al. 2019a;Ferreira et al. 2018;Martin et al. 2017). Secondly, the image delivered by an optical instrument depends on both the science object we want to characterize and the PSF. In order to estimate the interesting astrophysical quantities, one can either use a deconvolution technique (Fétick et al. 2020(Fétick et al. , 2019aBenfenati et al. 2016;Flicker & Rigaut 2005;Mugnier et al. 2004;Fusco et al. 2003;Drummond 1998) or include the PSF as part of a model, as is performed in PSF-fitting astrometry/photometry retrieval techniques (Beltramo-Martin et al. 2019a;Witzel et al. 2016;Schreiber et al. 2012;Diolaiti et al. 2000;Bertin & Arnouts 1996;Stetson 1987) and galaxy kinematics estimation (Puech et al. 2018;Bouché et al. 2015;Epinat et al. 2010) for instance.
Nevertheless, the PSF is not systematically straightforward to determine from the focal-plane image, particularly for observations made with a small field of view (FOV), or of extended objects, crowded populations or with a low signal-to-noise ratio (S/N). Alternative methods exist, such as PSF reconstruction (Beltramo-Martin et al. 2019b,a;Wagner et al. 2018;Gilles et al. 2018;Ragland et al. 2018;Martin et al. 2016a;Ragland et al. 2016;Jolissaint et al. 2015;Exposito 2014;Clénet et al. 2008;Gendron et al. 2006;Veran et al. 1997), which is potentially very accurate (1% error level). However, this process lacks science verification and requires years to be fully implemented and op-Article number, page 1 of 16 arXiv:2009.00994v1 [astro-ph.IM] 2 Sep 2020 A&A proofs: manuscript no. aanda erational as a stand-alone pipeline. The inherent complexity has inspired several efforts to produce simpler but still efficient reconstruction techniques Fétick et al. 2019b) that have been fully validated and are now in operation for Multi Unit Spectroscopic Explorer (MUSE) (Bacon et al. 2010) at the Very Large Telescope (VLT). The analytic AO PSF model described by Fétick et al. (2019b) has shown spectacular accuracy in reproducing the MUSE Narrow Field Mode (NFM) PSF as well as the PSF of the Zurich Imaging Polarimeter (ZIMPOL) (Schmid et al. 2018) that equips the Spectro-Polarimetic High contrast imager for Exoplanets REsearch (SPHERE) instrument (Beuzit et al. 2019) PSF and is an excellent candidate to rethink the way we perform PSF reconstruction, with a single flexible algorithm that complies with all kinds of AO correction. However, several questions arise as to its capabilities. Does this model really match any type of AO PSF, regardless of the observing conditions and even at very high Strehl-ratio (SR) ? How much accuracy can we expect in a statistical sense and what are the fundamental limits ? Can we improve the description of the PSF in the presence of static aberrations ? The goal of this paper is to address these questions by (i) providing an exhaustive demonstration that this model complies with any type of AO system, telescope, and in optical and near infrared (NIR) wavelengths, (ii) assessing conditions for which the model is optimal for representing an AO PSF, and (iii) showing that this framework is robust enough to enable the joint estimation of AO residual and instrumental aberrations. In Sect. 2, we present the theoretical background and the model description as well as some upgrades we have implemented so as to treat new types of data and include static-aberration fitting. In Sect. 3, we present a statistical analysis of the model accuracy over 4812 PSFs with a discussion on the parameter estimates. Finally, in Sect. 4, we illustrate how the present model allows joint estimation of atmospheric conditions, AO performance and static aberrations, caused by the low wind effect (LWE) on SPHERE or segment cophasing error at Keck.

PSF model
This model was originally introduced by Fétick et al. (2019b) and validated on SPHERE/ZIMPOL and MUSE NFM data. In order to increase the range of applicability of such a model, we upgraded it by including (i) the model of the Apodized Lyot Coronagraph (APLC, Soummer (2005)) used on SPHERE (focal-plane mask not considered) (ii) additional degrees of freedom to adjust static aberrations over a specific modal basis, which can be Zernike modes described on the pupil, tip-tilt and piston modes for each pupil area delimited by the spiders so as to describe the so-called LWE or piston modes for each pupil segment (Laginja et al. 2019;Leboulleux et al. 2018;Ragland 2018).
First of all, the model assumes the stationarity of the phase of the electric field in the pupil plane (Roddier 1981), which allows the system optical transfer function (OTF) to be split into a static parth Static (telescope + internal aberrations) and an AO residual spatial filterk AO as follows, where ρ is the separation vector within the pupil, λ is the imaging wavelength, and ρ/λ is the angular frequencies vector. The static OTF derives from the instrument pupil mask P(r) and the optical path difference (OPD) map ∆ Static in the pupil plane, which gives for a monochromatic beam at wavelength λ h Static (ρ/λ) = P P(r)P * (r + ρ) exp (2iπ/λ(∆ Static (r) − ∆ Static (r + ρ))) dr. (2) The SPHERE/IRDIS data we have treated (see Sect. 3.1 for a more complete description) were obtained during PSF calibration procedure using an off-axis stars. Therefore, the incoming beam was not getting through the focal-plane mask and the pupil mask model results from the multiplication of the apodizer (amplitude) function A and the Lyot stop L (Soummer 2005) as follows P(r) = A(r).L(r). ( The 2D functions A(r) and L(r) are illustrated in Fig. 1.
In this paper, we split the static aberration contribution in two terms: the Non Common Path Aberrations (NCPA) calibration (non-coronagraphic) noted ∆ NCPA (Jia et al. 2020;Vigan et al. 2019;Lamb et al. 2018;Vassallo et al. 2018;N'Diaye et al. 2016;Sitarski et al. 2014;Jolissaint et al. 2012;Sauvage et al. 2011;Robert et al. 2008;Mugnier et al. 2008;Blanc et al. 2003;Fusco et al. 2003), to which is added another contribution decomposed over a modal basis M to be specified. This latter could be simply Zernike modes over the whole pupil or per pupil area delimited by the spiders, so as to include potential LWE that may introduce severe asymmetries in the PSF structure. Analyses of SPHERE images (Milli et al. 2018;Sauvage et al. 2016Sauvage et al. , 2015 showed that the LWE mainly introduces piston, tip, and tilt differential aberrations between pupil areas separated by the spiders, which represents 12 parameters to be adjusted. Also, for segmented pupil telescopes like Keck telescopes, this basis can be segment piston or petal modes (Ragland 2018). We analyze the capacity of this model to identify such aberrations in Sect. 4. The static OPD is calculated following where µ Stat (k) and M k are respectively the coefficient and 2D shape of the k th mode over n m modes. To include segment piston errors, M k must be defined as the pattern formed by the k th segment (36 in total for the Keck pupil) and a k is the piston value in meters. Eventually, one may estimate both the atmospheric parameters and static coefficients and we illustrate such a joint estimation on Keck data in Sect. 4.1. Finally, the description of the AO residual spatial filterk AO in Eq. 1 relies on two assumptions, which are, (i) the exposure is infinitely long, and (ii) the residual phase is a Gaussian statistical process. Thus,k AO is fully described by the residual phase power spectrum density (PSD) as follows (Roddier 1981): where D AO , B AO , and W AO are the residual phase structure function, autocovariance function, and PSD respectively, k is the spatial frequencies vector and F [x] is the 2D Fourier transform of x. In Fétick et al. (2019b), the PSD is described as a split function to separate the AO-corrected spatial frequencies from the uncorrected frequencies that follow a Kolmogorov's −11/3 power law (Kolmogorov 1941) where k = |k|, r 0 is the Fried parameter (Fried 1966), k AO the equivalent AO cut-off frequency beyond which the AO correction no longer occurd, C a constant value and M an asymmetric Moffat function (Moffat 1969) that depends on five shape parameters A, α x , α y , θ, and β and on the spatial frequency vector k = (k x , k y ) as follows, where ψ is a normalization factor that ensures that the integral of the Moffat function over the AO-corrected area is given by A: (Fétick et al. 2019b) Contrary to the classical use of a Moffat function in astronomy (Trujillo et al. 2001;Moffat 1969), we stress that the Moffat function used in the present model serves in the description of the PSD and not the focal-plane PSF. This latter is deduced by the inverse Fourier transform of the OTF given in Eq. 1. Consequently, the PSF FWHM is actually driven by parameters A and r 0 rather than α x , α y , and β. In order to really capture the physical meaning of this model, below we present a description of each of the seven parameters that the model relies on and their impact on the PSF: r 0 constrains the uncorrected PSD that corresponds to the PSF wings, that is, the part of the PSF that remains untouched by any AO correction. For a system with n act × n act actuators to compensate for the incoming wavefront, this breaking occurs in the focal plane at approximately n act × λ/(2D), with D the primary mirror diameter and for a Nyquist-sampled PSF (Roddier 1999). -C is a constant value that changes the gap between the AOcorrected part and the uncorrected high-spatial frequencies.
This parameter plays a major role near the AO cutoff in both the PSD and PSF planes, and compares to the aliasing error that contaminates wavefront measurements due to the WFS discrete sampling (Bond et al. 2018a;Correia & Teixeira 2014;Jolissaint 2010;Rigaut et al. 1998).
-A is the total energy in nm 2 contained in the Moffat PSD model (constant C not included) and thus characterizes the AO residual wavefront error that connects to the PSF SR, that is, the intensity of the PSF peak compared to the diffractionlimit scenario.
α x , α y , and θ govern the elongation and skewness of the PSD as well as the direction of the elongation thanks to the angle θ. The PSD FWHM is proportional to α x , α y parameters; for the same amount of energy, a PSD with a larger FWHM indicates that the PSD flattens and the correction homogenizes across spatial frequencies.
β represents the asymptotic slopes of the PSD at large spatial frequencies within the AO correction radius. From Fourier analysis of AO performance (Bond et al. 2018a;Correia & Teixeira 2014;Jolissaint 2010;Rigaut et al. 1998), we know that the PSD pattern introduced by the wavefront measurement noise follows a k −2 power law for WFSs sensitive to the first order derivative of the wavefront such as the Shack-Hartmann . Asymptotically, we would have β = 1 for a Shack-Hartmann based AO system if the measurement noise is the only error in the wavefront reconstruction process. Moreover, in the situation of an AO system correcting all atmospheric aberrations with the same relative level (which does not occur as the sensitivity of the wavefront sensor is aberration-dependent (Fauvarque et al. 2016)), we would obtain β = 11/6 1.83 according to the von Kármánn expression of the atmospheric PSD (von Karman 1948). Nevertheless, larger values of β can be observed. Indeed, for very large values of β, the Moffat distribution converges towards a Gaussian shape (Trujillo et al. 2001), as does the PSF consequently. We expect this behavior with poor tip-tilt correction for instance or in the presence of strong telescope wind-shake or telescope/dome vibrations. As a summary, β should usually range between 1 and 1.83 for nominal AO correction but can reach higher values in the presence of nonatmospheric aberrations.
We illustrate in Fig. 2 azimuthal profiles of PSDs and corresponding PSFs for various sets of model parameters. This figure shows clearly how r 0 drives the PSF halo and how the parameter A increases the area underneath the PSD curve, which corresponds to the residual wavefront error. We also observe that a larger value of β sharpens the PSD. Knowing that the total energy remains constant, this situation corresponds to less efficient low-order modes compensation. Finally, a larger value of α increases the PSD FWHM and distributes more energy over the highest AO-correct modes. The shape of the PSF is mainly conditioned by parameters A and r 0 , which are the two most significant parameters of this model. Parameters C, α, θ, and β are secondary parameters that shape the PSD in order to more accurately reproduce the PSF structure in various observing conditions and especially in the presence of sub-optimal AO control or non-atmospheric aberrations. Henceforth, this model is a combination of a PSF determination tool and an AO diagnosis facility gathered up together into a single and parsimonious analytical framework.

Overview of data
We aim to provide the most exhaustive on-sky validation of the PSF analytical and parsimonious model proposed by Fétick et al. (2019b). To achieve this goal, we collected 4812 PSFs obtained from different observatories, AO systems, AO modes settings for the model parameters, assuming that the PSD is symmetric, i.e. α x = α y = α and θ = 0. No static aberrations were added here and the parameter C was set to C = 10 −4 rad 2 m 2 . The PSF maximum is normalized to the SR value. and spectral bands covering optical and NIR wavelengths. The GALACSI/MUSE and GeMS/GSAOI instruments deliver simultaneous PSFs at different wavelengths and field positions respectively. If we consider a single PSF for each observation they produce, we obtain a total of 1880 PSFs. However, although the realizations of atmospheric residual wavefront are not independent from a PSF to another within the same observation, these data allow to test the model in the presence of chromatic and field-dependent instrumental aberrations. In order to conserve the same fitting process for all data, we have fitted each PSF with the same starting point and regardless the results obtained on another PSF acquired during the same observation. A description of each data set is given below • SPHERE/ZIMPOL@VLT.SPHERE is a facility at the VLT that was installed in 2015 and relies on an extreme AO (EXAO) system (Sphere Ao for eXoplanets Observations; Fusco et al. 2016;Sauvage et al. 2015) to deliver a very efficient atmospheric correction level to three science instruments, which are ZIMPOL (Schmid et al. 2018 to extract the intensity image, subtract a bias frame, and correct for the flat-field. These data are particularly useful for testing the model in optical wavelengths and under strong atmospheric residual regime. • SPHERE/IRDIS@VLT. A total of 237 PSFs were obtained from the SPHERE Data Centre client (Delorme et al. 2017) using the Keyword Frame type set to IRD_SCIENCE_PSF_MASTER_CUBE. These PSFs were acquired over the last five years with the N_ALC_YJH_S APLC (Vigan et al. 2010) during PSF calibration and with a pixel scale of 12.5 mas/pixel. These data were collected using the dual band filters DB_H23, and DB_K12, and they are useful for testing the model under very high SR regime and for validating the LWE retrieval presented in Sect. 4.2. • Keck AO/NIRC2@Keck II. We obtained 355 PSFs using the narrow field mode of NIRC2 with a pixel scale of 9.94 mas/pixel and using the Fe II and K cont filters. The Keck AO system on the Keck II telescope was operated in single conjugated AO (SCAO) mode using a natural (Wizinowich et al. 2000) or an on-axis laser (Wizinowich et al. 2006) guide star. These PSFs were obtained during PSF reconstruction engineering nights in 2013 and 2017 (Ragland et al. , 2016 and such data are especially useful to validate the model under the influence of remaining piston cophasing errors, as we discuss in Sect. 4.1. • SOUL/LUCI@LBT These data were delivered in 2020 from the to two LUCI NIR spectro-imagers assisted with the pyramid SCAO (PSCAO) SOUL AO system (Pinna et al. 2016) driven by a pyramid WFS (Ragazzoni & Farinato 1999). These PSFs were acquired using the H filter (1.653 µm) and with a sampling of 15 mas/pixel. Although few (11 for our purpose) data sets have been obtained so far in comparison to others systems, these PSFs pave the way to a demonstration of PSF determination strategies in the presence of a pyramid WFS, which will be the baseline for the SCAO mode of HARMONI (Thatte 2017) • CANARY/CAMICAZ@WHT CANARY Myers et al. 2008) was designed as a pathfinder for demonstrating the reliability and robustness of the multiobject AO (MOAO) concept proposed for assisting very large field (>1') multi-object spectrographs, such as MO-SAIC for the ELT (Hammer et al. 2016 . We focused our analysis on the narrow field mode (NFM) of MUSE, which delivers a laser tomography AO (LTAO)-corrected field covering a 7.5" × 7.5" FOV, providing near-diffraction-limited images with a sampling of 25 mas/pixel. These PSFs were obtained during the commissioning phase in 2018 and are particularly useful for analyzing the model outputs with respect to the wavelength and demonstrating that the model also complies with under-sampled PSFs. Also, we performed a spectral binning to reach a spectral width of 5 nm (91 PSFs per cube if we remove the notch filter wavelengths), resulting in a total of 1986 PSFs.
• GEMS/GSAOI@GEMINI The Gemini South Adaptive Optics Imager (GSAOI) is a NIRd camera that benefits the correction provided by the Gemini Multi-conjugate Adaptive Optics ( In total, we created a dictionary of 4812 PSFs covering several observatories, AO correction types, optical and NIR wavelengths, as summarized in Table 1. Having this diversity of data is important for spanning the full range of possible AO correction levels and assessing which conditions must be met to achieve an accurate PSF representation.

PSF fitting
In order to fit the model over the image and retrieve the associated parameters, we followed the same strict process for each of the 4812 PSFs in our dictionary, which is described as follows -Define a model of the image including the PSD degrees of freedom (no static aberrations retrieval) µ AO as well as four additional scalar factors δ x , δ y , γ, and ν in order to finely adjust the PSF position, flux, and constant background level, The image model is calculated over a given number of pixels that is 10% larger than the on-sky images so as to mitigate aliasing effects due to the Fast Fourier Transform algorithm. The size of the on-sky image support is instrumentdependent and truncated in order to mitigate the noise contamination. When possible, we crop the on-sky image to conserve a FOV up to twice the AO cutoff. -Define a criterion to minimize where d i j is the (i, j) pixel of the 2D image and W i j the weight matrix defined by The weight matrix accounts for the noise variance, that is, both photon noise and read-out noise, and allows us to maximize the robustness of the fitting process (Mugnier et al. 2004). -Perform the minimization using a nonlinear and iterative recipe based on the trust-reflective-region algorithm (Conn et al. 2000). We did not use specific regularization techniques on top of the weighting matrix as we have taken care of selecting good S/N images.

PSF model accuracy
In Fig. 3, we present a visual comparison of on-sky PSFs, bestmatch model and residual map for the seven instruments considered in this analysis. The model adapts to any kind of AO correction; the maps of residuals are visually similar to first order among all systems. On the Keck AO/NIRC2 image, we see some structures that are static speckles, probably introduced by a remaining cophasing error, and residual NCPA as illustrated in Sect. 4.1. On SOUL/LUCI and CANARY/CAMICAZ images, we also see a persistent pattern that can be explained by the presence of static aberrations not included in the model and the exposure time that was not sufficiently long (a few seconds of exposure) to average the atmospheric speckles.
In Fig. 4 Fig. 4, there is no specific bias, except for SPHERE/IRDIS for reasons mentioned above, as well as for SOUL/LUCI owing to the small amount of data we have access to so far. Overall, we conclude that (i) the model becomes biased for very high SR observations (SR > 80%), calling for the introduction of instrumental defects to improve the model accuracy, (ii) the SR is estimated with a 1-σ precision of 1 %, and (iii) the FWHM is estimated with a 1-σ precision of 3 mas, which correspond to approximately to one-fifth down to one-tenth the pixel scale depending on the instrument.

Exploitation of the model outcomes for diagnosing observing conditions and AO performance
As discussed above, fitting the shape of the residual phase PSD allows to retrieve atmospheric parameters and AO performance. Therefore, the goals in this section are to (i) confirm that the output parameters r 0 and A follow expected trends and give confidence in the retrieval process, and (ii) provide statistics on α, β parameters to asses which values they should reach for a nominal AO correction and thus discriminate situations of sub-optimal AO correction. We have excluded the parameter θ from this analysis as it corresponds to a PSF orientation only.   Table 2: Individual statistics per system and imaging wavelength on SR and FWHM estimates. The columns "Median" give median values estimated on images, while columns "bias", "std" and "Pearson" give the median of SR/FWHM error (in percent) and the Pearson correlation coefficient respectively.

Primary parameters estimates
The seeing is estimated from the PSF wings fitting relying on a Kolmogorov expression of the PSD. This measurement technique has proven to be robust (Fétick et al. 2019b) as it consists in determining a single parameter from a significant number of pixels. Thanks to the large redundancy across the pixels of the spatial signature that the algorithms is seeking out, this approach still gives meaningful results with moderate S/N in comparison to external profilers (Fétick et al. 2019b). However, contrary to these latter, this focal-plane-based seeing determination includes all turbulence effects that contribute to impact the PSF wings, such as the dome seeing (Lai et al. 2019;Conan et al. 2019), which the external profilers are not sensitive to as they are apart from the dome. Consequently, estimating the r 0 from the focalplane image allows to diagnose more accurately the AO performance in comparison to an external profiler.
The target here is to verify that the retrieved seeing values are consistent with what we know about the observing sites. The first verification we made is to analyze statistics on the seeing estimates presented in Table 3. Given that the Keck data were acquired in February, August, and September seasons, the median seeing at Mauna Kea is consistent with the literature (Ono et al. 2016;Miyashita et al. 2004). Similar observations can be drawn for La Palma by comparing to either CANARY telemetry data (Martin et al. 2016b)

or the RoboDIMM (O'Mahony 2003).
We also find consistency with analysis by Masciadri et al. (2014) Article number, page 7 of 16  Table 3: Median and 1-σ standard-deviation of seeing values retrieved from the PSF-fitting process. Seeing values are given at zenith and 500 nm.
In addition, we highlight in Fig. 5 the r 0 estimates at zenith and 500 nm obtained on a single image of GeMS/GSAOI and a single cube of GALACSI/MUSE. For the GEMS/GSAOI case, we have r 0 measurements obtained from several PSFs distributed over the field. Given that r 0 is assessed from the PSF wings that contain the non-AO-corrected spatial frequencies, we expect r 0 estimates to be independent of the PSF position in the FOV. Similarly, we multiplied r 0 estimates obtained on GALACSI/MUSE images by a factor (500/λ c ) 6/5 , where λ c is the central wavelength of the considered image in the data cube, in order to disable the theoretical dependency of r 0 with respect to λ.
Using our model, we retrieve r 0 as 14.2 cm±0.4 and 10.4 cm±0.2 cm for GEMS/GSAOI and GALACI/MUSE, respectively. The achieved 3% 1-σ error of our estimates shows that this method is robust and precise compared to telemetrybased techniques that typically reach 10% precision .
For completeness, Fig. 6 shows histograms for the A parameter (the wavefront error) and the function A =f(r 0 ) for the specific case of CANARY working in SCAO (i.e., the largest SCAO PSF dictionary we have). We clearly see that the wavefront standard-deviation error follows a r −5/6 0 law as we expect from the von Kármán PSD (von Karman 1948), which is proportional to r −5/3 0 for a SCAO system whose on-axis PSF is not influenced by the atmospheric turbulence profile. This illustrates the agreement between the different retrieved parameters.
Moreover, the histogram of A estimates given in Fig. 6 shows that the retrieved wavefront errors vary within a meaningful range: the SPHERE AO system performs better than others, as expected from such an extreme AO system, with a median of 140 nm, which includes all the aberrations that blur the AO-corrected part of the PSF. The present SPHERE results gather PSFs obtained with SPHERE/ZIMPOL, acquired with a AO loop running at low frequency (300 Hz, m V 11), and SPHERE/IRDIS fir which 25% were contaminated by a strong LWE. Therefore, this wavefront error of 140 nm seems consistent with the literature . The Keck AO system reaches 230 nm, which also complies with Wizinowich (2012). From GALACSI/MUSE NFM images, we retrieve 285 nm of wavefront error, which agrees with the Oberti et al. (2018) analysis Moreover, the GALACSI system was not yet fully optimized to operate in LTAO and a recent acquisition in 2019 already showed improvements on SR. For CA-NARY/CAMICAZ, we obtain a wavefront errors of 320 nm, 396 nm and 417 nm in SCAO, MOAO, and GLAO mode respectively, which compares well with Martin et al. (2017) and Vidal et al. (2014). Finally, GeMS/GSAOI data unveil a median wavefront error of 570 nm from PSFs extracted from the field at 20" up to 60" off-axis, which corresponds to a SR of 7.8% at the edges of the field and complies with analysis of the performance of GeMS Rigaut et al. 2012).
We stress that this value of wavefront error is not deduced from the image SR but directly from the integral of the PSD, which can be determined from the model parameters (Fétick et al. 2019b). This corresponds to the real mathematical definition of the wavefront error and is not influenced by the Maréchal approximation (Parenti & Sasiela 1994) that is biased at low SR. Consequently, this model is also a robust focal plane-based wavefront error estimator.

Secondary parameters estimates
We have emphasized in Sect. 2 that parameters α and β govern the PSD shape and assist in the AO correction diagnosis. Figure 7 presents histograms for each system; the β parameter has a relatively narrow histogram with a peak located at β = 1.82 and a 1-σ standard deviation of 0.6, while the α histograms seem particularly instrument-dependent with values from 0.001 rad.m for efficient AO correction to 1 rad.m for marginal AO correction (e.g., in the bluest visible wavelengths of MUSE). Our first conclusion here is that an optimized AO system should provide a PSD with a β parameter comprised between 1 and 1.9, as discussed in Sect. 2. For larger values of β, there is an excess of residual error into low-order spatial frequencies. For instance, we see that on Keck AO/NIRC2 images the median β increases from 1.7 up to 2.1 in NGS and LGS modes, respectively, indicating the presence of additional loworder modes introduced here mainly from the focal anisoplanatism (Wizinowich 2012;van Dam et al. 2006). In addition, we have observed cases with β > 2.5 on Keck data due to a strong wind-shake effect that was enlarging the PSF FWHM by a factor three compared to nominal performance. Both α and β parameters grow up in the presence of strong wind shake and are wavelength dependent as we see on GALACSI/MUSE histograms. However, from GEMS/GSAOI data analysis, we do not find clear trend with respect to the field position, which may be owing to the uniform correction across the field provided by the MCAO mode of GEMS.
We illustrate here that those two parameters carry additional and relevant information on the AO correction. However, the exact connection with the AO status is not straightforward to identify. To enable this identification, we are currently developing a convolutionnal neural network (Herbel et al. 2018) that we train to estimate the model parameters from the AO control loop data, such as wavefront measurements. As we can collect a very large amount of data for the purpose of estimating a few tens of parameters, solving this problem is definitely achievable with dedicated simulations and on-sky data that we will continue to be collected regularly among observatories.

Discussion on the influence of exposure time and spectral bandwidth
One of the major assumptions in the model proposed by Fétick et al. (2019b) concerns the infinitely long exposure hypothesis. This model is therefore not capable of reproducing atmospheric speckles that average when taking a sufficiently long exposure.
As the method relies on the second-order statistical moment of the residual phase, the time necessary to achieve a convergence of the PSD shape is highly dependent on atmospheric parameters but may be achieved in few seconds (Martin et al. 2012). Thanks to SOUL/LUCI data, we analyzed the PSF accuracy when fitted  on short exposure images with integration times from 0.157 s to 60 s. We find that the model matches the PSF down to 1-2 s of exposure and for a PSF acquired at 1.6 µm, below which the parameters estimation begins to degenerate as illustrated in Fig.  8 (left) through the average of absolute parameters variations.
We have also noticed that the PSF shape remains stable from few seconds exposure, which explains the stability of retrieved parameters.
Moreover, using GALACSI/MUSE data, we tested the model on a polychromatic image, from 2.5 nm up to 470 nm of spectral width by binning monochromatic PSF observed simultaneously with MUSE. When compensating for the beam dispersion (recentering PSFs and then stack), the parameter estimation does not deviate by more than 4% over the whole spectral width span as presented in Fig. 8 (right). Consequently, the presence of chromatic static aberrations and atmosphere chromaticity do not prevent the model from very accurate characterization of the PSF. When including the beam dispersion (no recentering before stacking), which is 350 mas (14 pixels) for the largest spectral width, the parameter estimation degrades up to 10% and goes beyond the 4% threshold after 300 nm of spectral width. As a result, the model remains robust and reliable, even in the presence of such a strong beam dispersion.

Joint atmospheric, AO performance and static aberrations retrieval
The previous section illustrates that this model accurately and reliably reproduces the AO residual component of the PSF. Considering that the model is parsimonious (7 AO parameters and 4 photometric and astrometric parameters) and the large amount of pixels (10,000 for a 100×100 image) that the model-fitting may rely on, one can attempt to retrieve more degrees of freedom, such as static aberrations. In order to estimate the static coefficients, we replace in the criterion given in Eq. 10h Static byh Static (µ Static ), which becomes a parametric OTF term. Besides, determining the wavefront from a single long-exposure image suffer from strong degeneracies by essence and this approach does not overcome this problem by using a non-linear minimization algorithm to estimate the problem solution. We distinguish two sorts of degeneracy -Sign ambiguity : this occurs when trying to fit pair modes for which the PSF is not sensitive to the sign . Therefore, we concentrate the static aberrations retrieval on piston, tip and tilt modes only. Moreover, we bound the solution domain between -λ/2 to λ/2 in order to overcome phase wrapping issues with piston modes. -Local minima : two different wavefront patterns may produce a similar but slightly different PSF and make the PSFfitting algorithm converge towards a local minimum of the criterion. To assess how much this problem affect our algorithm, we have compared the PSF model with various conditions of piston, tip and tilt levels in order to highlight possible combinations (10,000) of modes that could produce a similar PSF. From a vector of aberrations, we have tested different permutations of the elements of this vector and in 99% of cases, the relative PSF variation is larger than 1%, which is large enough to change the structure of the PSF and retrieve the correct wavefront map, as long as the S/N is sufficient (>50). We also rely on an analysis from Gerwe et al. (2008) that shows that for a fit of 35 Zernike polynomials on a segmented pupil, there are no local minima as long as the initial guess for the static coefficient remains within ± 0.2λ = 330 nm rms in H-band from the optimal solution. We are in a situation where the static aberrations we attempt to retrieve are already mitigated from the Keck telescope active control (Chanan et al. 1988) and the VLT spiders coating (Milli et al. 2018). Consequently, the aberrations level we must retrieve remains sufficiently weak to mitigate the presence of local minima. At a level greater than 300 nm (Keck AO residual is 280 nm) would impact drastically the PSF and degrade the telescope science exploitation so much that this aberration would be necessarily visible and mitigated as much as possible.
Henceforth, we are in a good situation to estimate piston, tip and tilt static modes on Keck and VLT images obtained on bright star.

Keck cophasing error retrieval
The analysis presented in this section focuses on a smaller sample of 129 Keck AO/NIRC2 data acquired in NGS mode and with a high SR (>40%) in 2013 (Ragland et al. 2016). The residual NCPA map was calibrated at the beginning of the night. We compared model-fitting performance in three situations: (i) a fit of PSD and photometric and astrometric parameters (11 values), (ii) a fit of these parameters when including, following Eq. 2, the static aberration map calibrated at the beginning of the night and, (iii) a joint adjustment of the PSD, photometry/astrometry and the 36 piston values corresponding to each Keck pupil segment (47 parameters in total), including the calibration static aberration map. Figure 9 shows the PSFs and the residual maps produced from the three model-fitting strategies considered in this section The envelope shows the ± 4% of variations around the median value and the spectral bandwidth is systematically centered around 700 nm. Error bars are assessed by averaging results over the GALACSI/MUSE data sets. and obtained for three cases acquired in February, August, and September, 2013. We firstly conclude that accounting for the residual NCPA static aberrations improves the model fitting as revealed by the residual map. This is also confirmed by the results shown in Tab. 4 which gives the median and 1σ standard deviation of SR/FWHM estimates over the 129 data samples as well as the mean square error (MSE) determined from the residual map. We firstly observe that biases and dispersion values on the SR are larger than in Tab. 2, owing to the data sample we consider in the present analysis and that gathers only high SR data for which the residual atmospheric aberrations are weaker. This advocates for including the calibrated static aberrations map into the model in order to mitigate this bias in the model. From Fig. 9, we also observe that these static aberrations manifest mainly as static speckles in the PSF and consequently contribute marginally to shape the PSF core, which explain why the statistics on the FWHM estimates are not influenced by the reduction of the data sample. Moreover, the results shown in both Fig.  9 and Table 4 highlight that the model accuracy is slightly improved by including the NCPA residual map in the model, with still 50% improvement on the MSE. Furthermore, the cophasing map retrieval allows to drastically enhance the model-fitting results, for example by up to a factor three on the MSE compared to the first scenario, as also illustrated in Fig. 9 and Tab. 4. The fact that the model efficiently reproduces the speckles observed on sky images suggests that these speckles are mainly introduced by cophasing errors due to the telescope.
In addition, in order to further demonstrate the strength of this framework, we present in Fig. 10 the retrieved cophasing maps for three specific cases. The corresponding NCPA map used in the model is presented in Ragland et al. (2016). We unmistakably observe the presence of a stair mode, that introduces a constant step between adjacent segment rows or columns and  Table 4: Median and 1-σ dispersion of SR and FWHM errors and mean square error obtained over the 129 H-band and high SR (>40%) Keck/NIRC2 images treated in this analysis for the three implemented strategies: (i) PSD/stellar parameters retrieval, (ii) including the residual NCPA map in the model, and (iii) also including the estimation of Keck II cophasing errors. Systematically, accounting for the NCPA map allows us to slightly reduce the biases on the estimates and diminish the MSE, but best model-fitting results are obtained thanks to cophasing error retrieval that decreases the MSE by almost a factor three.
that has been already highlighted from data acquired in 2017 using an alternative technique (Ragland 2018). Moreover, each of those three images was obtained with different elevation-azimuth telescope positions that are reported in Table 5, which shows that the wavefront standard deviation value is very close to observations reported by Ragland (2018) for a 45 • telescope elevation. Also, the stair case peak-to-valley (PV) energy increases with low telescope elevation and rotates with respect to the telescope azimuth. This observation reveals a connection between the presence of this stair mode and global flexure over the primary mirror that is controlled to compensate for the gravity effect. Deeper analyses will confirm the presence of such a residual flexure on the Keck telescoped and whether or not this framework now to allow us to correct for this stair mode to generally improve the scientific exploitation of the Keck, but also future segmented telescopes.  Table 5: Summary of peak-to-valley (PV) and 1-σ standard deviation of retrieved cophasing error map regarding the corresponding telescope elevation/azimuth.
Measuring the piston map from the focal-plane image is feasible from results obtained with the present analysis, but necessitates to have a star bright enough in the field to calibrate these aberrations. We must ensure that the shape of the primary mirror is controlled without the aid of a focal-plane-based technique in order to achieve the best image quality regardless of the observed field. We are again in a situation where we have to retrieve a few parameters (47) from a potentially large amount of data delivered by segment position sensors, temperature, pressure and humidity sensors, atmospheric parameters near the dome and even AO telemetry, which can deliver information about the piston error when using a pyramid WFS (Bond et al. 2018b;Schwartz et al. 2018). It is not easy to investigate the connection between this ensemble of data and the piston map, and the use of neural networks for this purpose will be explored.

SPHERE low wind effect retrieval
the analysis presented in this section focuses on 176 SPHERE/IRDIS data obtained with good SR conditions (>40%) so as to test two different fitting strategies, namely (i) fitting the PSD and photometry/astrometry parameters (11 parameters) and (ii) fitting these latter 11, and 12 additional parameters corresponding to piston, tip, and tilt of each VLT pupil area delimited by the spiders in order to account for the LWE. According to Sauvage et al. (2016), this description of the LWE allows the main impact of this effect to be reproduced, i.e., a strong PSF asymmetry. Figure 11 presents the results of the two strategies implemented in our study over three particular cases for which a strong LWE is observed. We observe, especially in Fig. 11, that the sole parametrization of the PSD is not rich enough to precisely reproduce the strong asymmetry, which becomes a solved problem thanks to the aberration parametrization we propose. Regarding the estimated PSF shape and values given in Table 6, we show clear evidence that (i) describing the LWE as a combination of differential piston, tip, and tilt allows to accurately characterize the PSF asymmetry, and (ii) accounting for this description in the PSF model ensures that biases and dispersion on SR and FWHM estimates are drastically mitigated. Moreover, the average MSE value over the 176 data samples is significantly diminished, that is, by a factor two. On the three cases illustrated in Fig. 11,Table 7 also reports the SR and FWHM errors as well as the MSE obtained from the PSF-fitting. These results show clearly a gain on accuracy for these particular cases for which the LWE impact on the PSF is significant comparatively to the statistical trend observed on the 176 data sets. This gain is partially mitigated on the statistical analysis owing to the fact that only 25% of the data were noticeably contaminated by the LWE. Moreover, the atmospheric part of the PSF model can mimic PSF asymmetries by modifying the ratio α x /α y in Eq. 7. Despite it leads to an inaccurate representation of the LWE impact on the PSF, it permits to mitigate the SR and FWHM estimation error comparatively to a symmetric atmospheric PSF model. As a summary, fitting the twelve extra parameters to represent the LWE allows can enhance the PSF metrics estimation by a factor up to five.
In order to confirm the robustness of the LWE retrieval, we have compared the estimated static aberration map with ZELDA measurements ) taken during SPHERE commissioning nights in 2014. The PSF-fitting was performed using the differential tip-tilt sensor (DTTS) that delivers 32×32 pixels images (Baudoz et al. 2010;Sauvage et al. 2015) acquired simultaneously with SPHERE/IRDIS using the ZELDA focal-plane mask to measure optical aberrations within the pupilplane. Moreover, in order to calibrate ZELDA measurements (that are in ADU) to reconstruct the phase, we followed the process described by Sauvage et al. (2015) and removed the mean ADU value over the pupil and adjusted a multiplicative factor to obtain the closest PSF possible from the DTTS observation as reported in Fig. 12. According to this calibration, we obtained standard deviations of 173 nm rms and 146 nm rms on the ZELDA map and the DTTS image-based map,respectively, which leads to a quadratic difference of 92 nm rms. This residual includes the internal aberrations in the IRDIS science path that do not impact the DTTS images and reach 50 nm and standard deviation (Beuzit et al. 2019). Also, this residual is likely mostly due to higher order aberrations not included in the differential piston, tip, and tilt of the static map model, which suggests that further improvements can be pursued to characterize the LWE. In conclusion, the framework presented in this paper is a powerful tool to assist in the estimation of internal aberrations. Furthermore, we are able to assess these aberrations from the DTTS image that delivers a post-AO and non-coronagraphic image regardless the coronagraphic mask employed during SPHERE/IRDIS observations. Consequently, this technique allows the joint estimation of atmospheric and instrumental defects using tip-tilt sensors measurements. Future work will address the extension of this strategy to ELT instruments that will rely on a 2x2 Shack-Hartmann WFS to measure low-order modes and that will provide PSFs from which we will be able to assess the telescope aberrations and calibrate a PSF model for science exploitation.

Conclusions
This paper revisits and improves the analytical framework proposed by (Fétick et al. 2019b), which now includes a parametrization of static aberrations for a joint retrieval of atmospheric parameters, AO performance and static aberrations.
We demonstrate in an exhaustive manner, using 4812 PSFs obtained from four different observatories and seven optical or NIR instruments, that the proposed model matches the PSF of any AO flavor within 4% error, even for high-SR observations. We also illustrate that the retrieved parameters carry relevant information about the AO performance and the atmospheric conditions, especially seeing and wavefront error, that shows agreement with the literature.
Finally, we illustrate that this model, upgraded with additional degrees of freedom to estimate static aberrations, allows the atmospheric parameters, AO performance and static aberrations, to be retrieved simultaneously. Particularly, the framework presented in this paper allows us to assess (i) the Keck pupil segment piston errors and especially the presence of a stair mode that was already pointed out (Ragland 2018) and, (ii) the LWE on SPHERE/IRDIS images as a combination of differential piston, tip and tilt over the four VLT pupil segments delimited by the spiders.
This model is a unique tool that gathers an AO diagnosis and a PSF estimation facility in the simplest and the most parsimonious way possible. However, we have handled it as a parametric model so far and the next step of this work is to enable a forward estimation of its parameters from contextual data expect the focal-plane image. We emphasized that the connection of these parameters with the observing conditions is not easily made; nevertheless, thanks to the parsimony of this model, this regression problem consists in assessing a few parameters (up to a few tens) from a large amount of data, which is provided by the AO telemetry, all the sensors within the telescope and the dome and the external meteorological profilers, which can be achieved with the use of neural networks. We are currently developing convolutional neural networks capable of directly estimating the model outputs from either an imaged PSF or a subsample of AO telemetry and we will present this work in a dedicated publication.  Table 6: Median and 1-σ dispersion of SR and FWHM errors obtained over the 176 H-band and high SR (>40%) SPHERE/IRDIS images treated in this analysis for the two implemented strategies: (i) PSD/stellar parameters retrieval only (11 parameters in total) and (ii) including piston, tip, and tilt retrieval for the four pupil segments (23 parameters in total). Systematically, accounting for the LWE model allows us to reduce the biases and the dispersion on PSF estimates.
Case LWE ∆ SR (pts) ∆ FWHM (mas) MSE (%)  Table 7: Impact of the LWE fit compared to a pure atmospheric model on SR, FWHM and MSE metrics for the three cases presented in Fig. 11.