Physics-based model of the adaptive-optics corrected point-spread-function

Context. Adaptive optics (AO) systems greatly increase the resolution of large telescopes, but produce complex point spread function (PSF) shapes, varying in time and across the field of view. This PSF must be accurately known since it provides crucial information about optical systems for design, characterisation, diagnostics and image post processing. Aims. We develop here a model of the AO long exposure PSF, adapted to various seeing conditions and any AO system. This model is made to match accurately both the core of the PSF and its turbulent halo. Methods. The PSF model we develop is based on a parsimonious parameterization of the phase power spectral density with only five parameters to describe circularly symmetric PSFs and seven parameters for asymmetrical ones. Moreover, one of the parameters is directly the Fried parameter r0 of the turbulence s strength. This physical parameter is an asset in the PSF model since it can be correlated with external measurements of the r0, such as phase slopes from the AO real time computer (RTC) or site seeing monitoring. Results. We fit our model against endtoend simulated PSFs using OOMAO tool, and against on sky PSFs from the SPHERE ZIMPOL imager and the MUSE integral field spectrometer working in AO narrowfield mode. Our model matches the shape of the AO PSF both in the core and the halo, with a sub 1 percent relative error for simulated and experimental data. We also show that we retrieve the r0 parameter with subcentimeter precision on simulated data. For ZIMPOL data, we show a correlation of 97 percent between our r0 estimation and the RTC estimation. Finally, MUSE allows us to test the spectral dependency of the fitted r0 parameter. It follows the theoretical lambda 6 fifth evolution with a standard deviation of 0.3 cm. Evolution of other PSF parameters, such as residual phase variance or aliasing, is also discussed.


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

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

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 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 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 where the notation M 1 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.

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 where λ is the observation wavelength, h the total OTF, hT the telescope OTF, and hA 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 with B φ the phase autocorrelation function defined as The Wiener-Khintchine theorem states that the PSD is the Fourier transform of the autocorrelation as where W φ denotes the phase PSD and F the Fourier operator; f and ρ are the Fourier conjugated variables.If we call u the angular variable conjugated to ρ/λ, the PSF is written as 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 hT , 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 hT term for more accuracy.The term hT being fully determined, now we only have to parameterize the residual phase PSD to model the PSF.

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: The uncorrected area is not affected by the AO system, and the phase PSD consequently follows the Kolmogorov law, where r 0 is the Fried parameter scaling the strength of the turbulence.The halo is thus set by the knowledge of only this r 0 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 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 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.Our PSF model based on the PSD is made of the following set of seven parameters: S = {α x , α y , β, θ R , C, r 0 , 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 S 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, The Strehl ratio consequently is written as

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: where h i, j is the discretized model of PSF on the pixels (i, j), S 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 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 L, which is the neg-logarithm of the likelihood for a Gaussian process.
Let us now note that the minimum of L 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 S and the position parameters (δ x , δ y ) as where γ and ζ are the analytic solutions for the flux and the background, respectively.At each iteration of the minimization process, the minimizer evaluates our L criterion with a new set of parameters (S, δ x , δ y ).The current PSF estimate h(S, δ x , δ y ) is computed, then the analytic solutions γ and ζ are computed.The quantity γ • h(S, δ x , δ y ) + ζ is used to compute the residuals with the data d.Residuals are then provided to the minimizer to estimate a new set of parameters (S, δ x , δ y ).

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 where we have chosen the outer scale L 0 = 30 m.Since 1/L 0 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 0 from 7.5 cm to 25.0 cm.All r 0 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 = (λ 1 /λ 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 seconds.For 0.1 and 1 second, the random OOMAO phase screen did not converge towards a stable state, leading to a strong bias in the r 0 estimation.For a 10 second exposure PSF, the random fluctuations of the phase are correctly averaged.This is confirmed by the 100 second exposure PSF, which gives the same r 0 estimation as the 10 second case.Since a 100 second exposure is computationally demanding and does not significantly improve the results, we performed all our tests on the 10 second exposure time.

Parameter
Table 2. Typical range of PSF parameters for OOMAO simulations and SPHERE/ZIMPOL instrument, lower bounds and values used as initial guess for the minimizer.Typical ranges are indicative and may vary according to the considered instrument.The value 'eps' denotes the machine precision.Parameters do not have any upper bound.
Each PSF is then fitted using the L (S, δ 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 0 = 18 cm.Using the same initial parameters {S, δ 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 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 −3 with a standard deviation of 2.2 × 10 −3 .For comparison, fitting with a Moffat model gives an average relative error of h = 3.4 × 10 −2 with a standard deviation of 1.6 × 10 −2 .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 γ analytically estimated by our model fitting method, and the OOMAO flux that is directly the sum of the data on all the pixels: 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%].

Fried parameter r 0 estimation
As shown in Fig. 4, our r 0 estimation is consistent with the OOMAO value of r 0 .We find the best linear fit to be r 0,FIT = 1.038 r 0,OOMAO − 0.132, where values are given in centimetres.The Pearson correlation coefficient is C Pearson = Cov(r 0,FIT , r 0,OOMAO )/ Var(r 0,FIT ) • Var(r 0,OOMAO ) = 0.99992.This result fully confirms our r 0 estimation with respect to OOMAO simulations with sub-centimetre precision.

AO estimation
Theoretically our model should also be able to retrieve the residual variance σ 2 AO on the corrected area and follow a λ −2 power  law. Figure 5 shows the fitting estimation of this variance versus the wavelength.This data is then fitted with curves of equation aλ −2 .Except the two outliers for minimal r 0 = 7.5cm at low wavelength (λ 500nm), the λ −2 power law is a good estimation of the σ 2 AO 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 σ 2 AO parameter for PSF estimation.The λ −2 spectral dependence is also an asset to shift the PSF from one wavelength to another.

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 −5/3 0 (Rigaut et al. 1998), we look for similar r 0 dependencies for the PSF constant C. Figure 6 shows a clear decrease of C with r 0 .Fitting the estimated C with a r −5/3 0 power law shows a good match, with small residuals for nearly all r 0 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 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.

High performance imager ZIMPOL
The Spectro Polarimetric High-contrast Exoplanet REsearch (SPHERE) instrument (Beuzit et al. 2008;Beuzit et al. 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 0 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 0 given by our fitting and r 0 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 0 of the sample, the median r 0 , and the largest r 0 , 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 −3 , with a standard deviation of 1.2 × 10 −3 .We only consider diffraction due to the telescope 8 m aperture and its central obstruction in the static OTF hT (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 hT 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 hT might bias the atmospheric term hA during fitting procedure and slightly offset the estimation of the S 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.
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 0 estimated by fitting and SPARTA is r 0,SPARTA = 3.41 r 0,FIT − 16.82. ( The Pearson correlation coefficient between the two series r 0,FIT and r 0,SPARTA is C Pearson = Cov(r 0,FIT , r 0,SPARTA )/ Var(r 0,FIT ) • Var(r 0,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 0 estimations and telescope r 0 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 minutes (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 0 with respect to the MASS/DIMM r 0 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 0 ), 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 0 , 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.

MUSE integral field spectrograph
The Multi-Unit Spectroscopic Explorer MUSE (Bacon et al. 2006(Bacon et al. , 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 5nm 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 −3 .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.
The evolution of r 0 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 0 gives a spectral averaged estimation of r 0 = 13.3 cm at 500 nm.Our fitted r 0 matches well the theory, with a standard deviation of 0.3 cm.(orange line).The best match between data and theory is achieved for r 0 = 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 0 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 0 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 C Pearson = 0.96, which is similar to the SPHERE/ZIMPOL case.However, we get the linear relationship r 0,SPARTA = 1.96 r 0,FIT − 5.31, 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 0 is set to 25 m, whereas for MUSE the L 0 is estimated jointly with r 0 .This SPARTA double trend is corroborated by MASS/DIMM information (Fig. 8, middle graph).On the other hand, r 0 estimation using our PSF model gives similar results on both SPHERE/ZIMPOL and MUSE instruments.This confirms the robustness of our PSF fitting method.

Conclusions
In this article we developed a parametric model of longexposure 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 0 and the residual AO variance σ 2 AO .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 0 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).where A is the 2 × 2 matrix of the system defined as w i, j 1 h i, j h i, j h 2 i, j . (B.4) In order to invert A, we need to make sure that det(A) 0. The determinant of A is This states that the determinant of A 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(A) 0 and the analytic solution for γ and ζ is written as w i, j d i, j 1 h i, j .(B.8)

Fig. 1 .
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 insert plot is a zoom on the low spatial frequencies.

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

Fig. 3 .
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 r 0 .Solid curves: OOMAO.Dashed: fitting.Dotted: residuals.All PSFs, for all wavelengths, are sampled at Shannon-Nyquist.

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

Fig. 5 .
Fig. 5. Estimation of the σ 2 AO from PSF fitting versus the wavelength (dots).Colours correspond to the seven different values of OOMAO r 0 .Curves of parametric equations σ 2 = aλ −2 are fitted on the data.

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

Fig. 7 .
Fig. 7. Three ZIMPOL PSFs (top), model fittings (middle), residuals (bottom).Left: Minimal r 0 of the sample.Middle: Median r 0 .Right: Maximal r 0 .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.

Fig. 8 .
Fig. 8. Zenital r 0 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: r 0,SPARTA versus r 0,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: r 0,SPARTA and r 0,FIT versus r 0,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).

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

Fig. 10 .
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 r 0 = 13.3 cm.The grey area corresponds to missing wavelengths due to the sodium notch filter.
Consequently, in order to get a Moffat of unit integral over the whole plane, one must choose the Moffat amplitude factor as The minimum of L (given in Eq. 15) has an analytic solution for γ and ζ since nulling the partial derivative of L towards these parameters gives∂L ∂γ = 0 ⇐⇒ i, j w i, j h i, j γ • h i, j + ζ − d i, j = 0 j γ • h i, j + ζ − d i, j = 0 ⇐⇒ γ i, j w i, j h i, j + ζ i, j w i, j = i, j w i, j d i, j .(B.2)These two equations are linear in γ and ζ.They can be written within the matrix formalism

Table 1 .
OOMAO parameters summary for our PSF simulations.