Reconstruction of the ground-layer adaptive-optics point spread function for MUSE Wide Field Mode observations

Here we describe a simple, efficient, and most importantly fully operational point-spread-function(PSF)-reconstruction approach for laser-assisted ground layer adaptive optics (GLAO) in the frame of the Multi Unit Spectroscopic Explorer (MUSE) Wide Field Mode. Based on clear astrophysical requirements derived by the MUSE team and using the functionality of the current ESO Adaptive Optics Facility we aim to develop an operational PSF-reconstruction (PSFR) algorithm and test it both in simulations and using on-sky data. The PSFR approach is based on a Fourier description of the GLAO correction to which the specific instrumental effects of MUSE Wide Field Mode (pixel size, internal aberrations, etc.) have been added. It was first thoroughly validated with full end-to-end simulations. Sensitivity to the main atmospheric and AO system parameters was analysed and the code was re-optimised to account for the sensitivity found. Finally, the optimised algorithm was tested and commissioned using more than one year of on-sky MUSE data. We demonstrate with an on-sky data analysis that our algorithm meets all the requirements imposed by the MUSE scientists, namely an accuracy better than a few percent on the critical PSF parameters including full width at half maximum and global PSF shape through the kurtosis parameter of a Moffat function. The PSFR algorithm is publicly available and is used routinely to assess the MUSE image quality for each observation. It can be included in any post-processing activity which requires knowledge of the PSF.


Introduction
Achieved image quality is usually the primary parameter of successful observations, especially those performed at groundbased telescopes where atmospheric turbulence produces highly changeable conditions. For many scientific applications, precise knowledge of the achieved image quality is an absolute prerequisite. An obvious example is the comparison with higher spatial resolution space observations, like those obtained with the Hubble Space Telescope, which achieve a ten times higher resolution than classical ground-based observations in median seeing conditions. Source optimal extraction, source deblending, and image deconvolution are other examples where accurate knowledge of the point spread function (PSF) is needed (Beltramo-Martin et al. 2019;Fétick, R. JL. et al. 2019).
In natural seeing observations, the PSF full width at half maximum (FWHM) is often used to quantify the achieved image quality. Most modern ground-based telescopes are equipped with a seeing monitor which provides a real-time estimate of the FWHM. This is very convenient to get a rough estimate of the PSF, but it is usually not accurate enough in a number of science applications. Firstly, the PSF is obtained at zenith and at a given wavelength, and these parameters are usually different from the airmass and wavelength of the observation. Secondly, the measurement is taken with a small telescope and does not take into account the relative outer-scale size of the turbulence with respect to the size of the telescope, or the image quality of the telescope plus instrument system.
The easiest and best method to obtain a good estimate of the PSF is to take an a posteriori measurement of a bright, unresolved, and isolated source on the final image or data cube. This is an advantage as it takes into account all the possible effects that can alter the PSF: the atmospheric turbulence but also the instrument finite resolution, the detector sampling, and even some inaccuracy of the data reduction chain. For natural seeing observations, several more or less elaborated models exist, such as for example the Moffat function or the multi-Gaussian function (Bendinelli et al. 1987;Trujillo et al. 2001;Infante-Sainz et al. 2019). However, in some cases there are no bright point sources in the field of view because it is too small and/or is located at high galactic latitude where Galactic stars are rare. This is for example the case of deep-field observations like the Hubble Ultra Deep Field (UDF, Bacon et al. 2017). A&A proofs: manuscript no. aa_2020_37595 The main challenges for extra-galactic observations are twofold. Firstly, extra-galactic observations require long exposures over the course of different nights. Variation of the PSF over the whole observation campaign can potentially be significant: for instance, in its Wide Field Mode (WFM), the Multi Unit Spectroscopic Explorer (MUSE) PSF can change by more than 100% over several nights. At the same time, cosmological fields are usually devoid of point sources that can be used to monitor this variability (Damjanov et al. 2011). Variation of the PSF then becomes a major limitation in kinematic or morphological analyses of distant galaxies. As such, Bouché et al. (Bouché et al. 2015), using state-of-the-art morpho-dynamical 3D algorithms, show that the PSF FWHM must be known to better than 20% so as not to degrade the velocity parameters (maximum velocity, dispersion) by more than 10%. These latter authors also highlight that the shape of the PSF (ellipticity) is critical for morphological parameters such as the inclination. There is a known degeneracy between rotational velocity and inclination of the system (Wright et al. 2009;Epinat et al. 2010), and as such, the PSF ellipticity must be known to 10% or better. In this case, analytical PSF models are sufficient and better knowledge of the PSF is not critical as the accuracy of the analysis is limited by the morphokinematical model and/or signal-to-noise ratio.
In some other cases where the density of stars is too high, like in the centres of Globular Clusters, there are no isolated stars and more sophisticated techniques are needed to infer the PSF (Kamann et al. 2018b,a). Even when a bright and isolated star is present in the field of view, the method described above assumes a uniform PSF over the field of view which is not always the case. We note that for seeing-limited observations with a limited field of view, the atmospheric and telescope PSF can indeed be considered as uniform with the field of view, but this is not generally the case for the instrument. With the generalisation of advanced adaptive optics (AO) systems, modern ground-based telescopes now offer improved image quality with respect to the seeing characteristics of the telescope site. However, the AO PSF is no longer a simple function of a single atmospheric parameter (seeing) but it is a complex function of the atmospheric turbulence (e.g. the profile of the atmospheric turbulence, the coherence time, and the anisoplanetic angle), the number of actuators, the accuracy and speed of the deformable mirror and the wavefront measurements, the brightness and location of the natural and laser guide stars, and the wavelength of observations. Despite this apparent complexity, AO systems offer a unique advantage in that they measure the atmospheric turbulence in real time at the exact location of the observation and through the same system used for the scientific observation. Thanks to the wave-front sensing telemetry information and good knowledge of the system, it is theoretically possible to predict the PSF, even without a point source within the field of view (Véran et al. 1997).
Although PSF-reconstruction algorithms have been in existence for a long time (Véran et al. 1997;Gendron, E. et al. 2006;Gilles et al. 2012;Ragland et al. 2016;Beltramo-Martin et al. 2019), most of them are too complex to implement and are not robust enough to be used blindly in normal operations. The fact that AO techniques have also evolved rapidly in parallel giving birth to a large number of species (e.g. single conjugate AO, ground layer AO, laser tomography AO, multi conjugate AO) has also prevented the development of robust and stable PSF-reconstruction algorithms and their validation on sky. However, today the situation has changed with the advent of the ESO Adaptive Optics Facility (AOF, Arsenault et al. 2008;Oberti et al. 2018) at the Very Large Telescope (VLT) which is now in regular operation at UT4 since 2017 with the HAWK-I (Pirard et al. 2004) and MUSE (Bacon et al. 2004(Bacon et al. , 2014a instruments. With the regular use of MUSE Ground Layer AO (GLAO) mode, the number of non-AO expert users has increased significantly, and the need for an efficient and robust PSFreconstruction system is becoming more and more important. Another motivation is the need to qualify and rank the observations during service mode operation. The previous scheme based on the seeing monitor information cannot be used, as other critical information such as ground-layer fraction must be taken into account.
Here, we present a PSF-reconstruction algorithm developed specifically for the GLAO mode of MUSE. The paper is organised as follows. After a brief presentation of the MUSE instrument and its wide-field mode, we describe the AOF module that allows users to correct for the ground-layer contribution of the atmosphere, significantly improving the MUSE final images. We then present our PSF reconstruction scheme, its specificity, and its optimisation with respect to the typical performance of MUSE-WFM and AOF-GLAO correction. A sensitivity analysis using End-to-End simulations is provided and some algorithm parameters are then adjusted accordingly. After demonstrating the algorithm performance on simulated (and thus well-mastered) data, it is applied to a real MUSE on-sky observation. Thanks to more than one year of Globular Cluster data, we are able to provide a statistical analysis of our algorithm in operational conditions. Final results of this study are reported here with a clear demonstration of the efficiency of our final PSF-reconstruction (PSFR) approach. The final implementation in the MUSE pipeline is detailed and a first astrophysical application to the MUSE Ultra Deep Field observation is presented as an illustration of the importance and the power of the unique combination of PSF reconstruction and GLAO corrected wide field MUSE 3D cubes (Bacon et al. 2017).

MUSE Wide Field Mode
MUSE is the ESO VLT second-generation wide-field integral field spectrograph operating in the visible (Bacon et al. 2014b), covering a simultaneous spectral range of 480-930 nm with a spectral resolution of ∼3000. Its Wide Field Mode (WFM) offers a field of view of 1 arcmin 2 , sampled at 0.2 ′′ . . MUSE is composed of 24 identical channels, each one comprising a single Integral Field Unit (IFU) with an image slicer, a spectrograph, and a 4k×4k CCD. MUSE has been in regular operation since October 2014. It was used in natural seeing mode until October 2017 when its GLAO mode was als made available to the ESO community.

GALACSI-GLAO and its specificities
The MUSE GLAO mode is performed by the Ground Atmospheric Layer Adaptive Optics for Spectroscopic Imaging (GALACSI) and is part of the AOF, a full AO system with a deformable secondary mirror of 1170 actuators, four 20-Watt laser guide stars, and two wave-front sensing units (GALACSI and GRAAL) at each Nasmyth Platform, feeding MUSE and Hawk-I, respectively. Commissioned in 2017, GALACSI provides improved image quality (e.g. 10-50% improved FWHM) in the MUSE wavelength range (480 -930 nm) and over the full of MUSE wide field mode. The MUSE WFM image-quality requirements, mainly driven by the instrument sampling and field of view (FoV), have led to a very specific correction stage which is focused on ground-layer correction. Hence, it does not try to reach the telescope diffraction limit but rather to improve the 'equivalent' seeing (Oberti et al. 2018;Madec et al. 2018). In that respect, both the AO error budget breakdown and the PSF shape are very different from the classical AO ones. In addition, the performance criterion is no longer the Strehl Ratio but rather some parameters related to the PSF shape such as its FWHM or ensquared energy in various box sizes. Even though these parameters could be related to the residual wave-front variance, the relation is far less obvious than the for the Strehl Ratio and this new paradigm in terms of PSF shape and performance has to be analysed and taken into account in the PSF reconstruction scheme.
Analysis of the non-AO PSF on MUSE sky data has shown that a circular MOFFAT gives an accurate description of the PSF core and wing (Moffat 1969;Andersen et al. 2006;Müller Sánchez et al. 2006). The smooth evolution of the Moffat shape parameters (FWHM, β) with wavelength can also be fitted with a low-order polynomial. This model has been extensively used with success for science analysis since MUSE began operation and its relevance has been fully demonstrated. From this basis and because the GLAO mode only provides a very partial correction and the resulting PSF is far from being diffraction limited. GLAO correction can be seen as a seeing improvement and in that respect a simple yet efficient way to describe a GLAO-corrected PSF is still to consider a Moffat function that can be fully described by its FWHM and its kurtosis (β coefficient which characterises the wing shape of the PSF). The FWHM could be non-symmetrical if we need to account for residual anisoplanatism effects. More details and a justification of the Moffat choice for the GLAO PSF parameters is given in Sect. 2.3.
For the sake of simplicity and clarity, let us now focus on the FWHM and derive an error budget for a typical MUSE WFM PSF.
FWHM Tel stands for the telescope diffraction and any defect related to its aberrations that will not be corrected by the AO stage (high-temporal-frequency wind shake and/or vibrations, field aberrations, etc.); FWHM Atm, GLAO stands for the FWHM extension due to the uncorrected part of the atmospheric residual phase and FWHM MUSE stands for the FWHM extension due to the MUSE configuration: its coarse sampling and its own internal optical aberrations (not corrected by the AO loop). In the following we assume that FWHM T el and FWHM MUSE are constant whatever the observation, whereas FWHM Atm, GLAO is a timedependant contribution (this assumption is discussed in Sect. 5.2. There is no simple, analytical way to link FWHM Atm, GLAO and σ 2 Atm, GLAO (the residual variance after AO correction). Nevertheless, it is straightforward to say that they follow the same monotonic behaviour. In other words, being able to identify the dominant terms of the residual variance will give us critical items for developing an efficient model of the GLAO part of the MUSE PSF. In that respect, a full GLAO error budget can be developed as follows: Atm, GLAO = σ 2 High Order modes + σ 2 Tip tilt , with σ 2 High Order modes = σ 2 fitting + σ 2 aliasing + σ 2 High Layers contrib and σ 2 Tip tilt = σ 2 TT, aliasing + σ 2 TT, anisoplanatism + σ 2 TT, noise + σ 2 TT, tempo (4) In the above error budget list, which gathers all the known error terms for such an instrument (due to spatial sampling, measurement noises, temporal error, anisoplanatism, etc.), there are several points worth highlighting: The fitting error mainly acts on high spatial frequencies, that is, those higher than the AO cut-off frequency defined as 1/(2d), with d being the deformable mirror (DM) spacing on the PSF wings for example. The laser guide stars (LGS) are bright enough (20 Watts emitted on sky) to neglect the measurement noise on LGS WFS (σ 2 noise ≃ 0). As shown in Fig. 1, σ 2 High Layers contrib is typically 100 times larger than the other error terms (aliasing and temporal effects).  High Layers contrib and other AO error terms (σ 2 noise , σ 2 aliasing and σ 2 tempo ) as a function of seeing. The grey and blue areas correspond to the possible variations of atmospheric parameters (wind speed, C 2 n profiles, outer scales) for the various error items.
The tip tilt (TT) contribution can be decomposed into several terms among which noise and anisoplanatism are by far the dominant ones. The choice of natural guide star (NGS) in the technical FoV is mainly driven by its limit magnitude. The WFS characteristics could also be adapted to accommodate low flux NGS by changing the integration time. These two combined aspects lead to observing configurations where the noise term is never dominant in the error budget. The noise is neglected in the following; we note that it would have been straightforward to take it into account in our algorithm by adding a classical noise measurement and noise propagation term. A combination of simulations, AOF design, and AOF on-sky data shows that such an addition is of no real benefit in the MUSE WFM case.
By design the NGS is always further than 1 arcmin from the optical axis. In that case, a good approximation assuming a twolayer model for the turbulence is to consider full decorrelation of the high-layer contribution and a full TT correction of the ground layer. This is confirmed by Fig. 2 where the TT anisoplanatism contribution is plotted for various atmospheric conditions. It is shown that the decorrelation hypothesis is very well validated but also that the final TT contribution due to the anisoplanatism effects remains very small (typically smaller than 50 mas) with respect to the MUSE pixel size (200 mas).
A&A proofs: manuscript no. aa_2020_37595 Impact of TT anisoplanatism for strong seeing condition and low ground Layer In all cases, the full decorrelation of TT anisoplanatism combined with the small contribution of TT noise leads to a nonelongated PSF (this has been experimentally confirmed on all the MUSE WFM images since the beginning of the instrument operation more than three years ago). The TT star is only here to ensure that ground-layer contributions of the atmosphere, the telescope pointing, and wobble aspects are corrected for.
The analysis of the various error terms clearly shows that in the specific case of GLAO, the PSF is mainly impacted by three terms: the 'fitting' and 'high layer (HL) contribution' terms for the LGS (i.e. the High Order mode correction), and the TT anisoplantism for the NGS contribution. The latter three terms are the only ones considered in the following in our GLAO PSFR algorithm.

End-to-End simulation of the MUSE ground-layer adaptive-optics system
End-to-End (E2E) simulations are carried out with Object-Oriented Matlab Adaptive Optics -OOMAO (Conan & Correia 2014), which is a Matlab communitydriven toolbox dedicated to AO systems. The OOMAO toolbox is based on a set of classes representing the source, atmosphere, telescope, wave-front sensor(WFS), Deformable Mirror (DM), and an imager of an AO system. It can simulate NGS and LGS single-conjugate AO (SCAO) and tomography AO systems on monolithic and/or segmented telescopes up to the size of the Extremely Large Telescope (ELT).
We used OOMAO for simulating the full AOF system in order to validate our PSFR algorithm and to provide a comprehensive sensitivity analysis of the PSFR performance. The simulation parameters (from the system and the environment view points) are listed below. Figure 3 presents the problem geometry and more precisely the positions of LGS, NGS, and the directions of interest in the FoV. Tables 1 and 2 present the system and atmospheric parameters used in the simulation and sensitivity analysis.

Choice of PSF model
One of the critical issues for any PSF reconstruction algorithm is the choice of PSF model. The most natural basis for describing the PSF is a pixel-wise basis. Although by definition, working on a pixel-wise basis removes the need for a model-based approximation, it is often not very well adapted to operational constraints because it requires a lot of memory and storage capacity. This is especially true for multi-wavelength instruments where one (or several if there are field variations to account for) PSF has to be computed and stored per wavelength bin. The choice of PSF model also strongly depends on the type of AO system (and thus AO correction) considered as well as the observation and post-processing requirements related to the astrophysical science cases. For MUSE WFM, the GLAO system only provides a very partial AO correction, and therefore the two critical parameters that have been identified to cover most of the science case requirements in terms of PSF knowledge are the PSF FWHM, which provides information on the data quality and final image resolution; and the level of the PSF wings, which provides information on the energy spread by the PSF in the FoV (spaxel contamination).
Considering these two aspects and the typical shape of a very partially GLAO-corrected PSF, a Moffat model is particularly well adapted for the description of the MUSE WFM PSF. The Moffat PSF can be mathematically described as follows.
where α x (resp. y) and β are directly related to the function FWHM, M 0 stands for the global amplitude factor, and m x , y for the absolute focal plane positions. Furthermore, where β is a very good marker of the shape of the PSF wings; the poorer the correction, the larger the β. It has been shown that a Moffat model with a β value of greater than 4 is very well adapted for describing a purely turbulent PSF (Trujillo et al. 2001). Figure 4 shows a comparison of a simulated GLAO PSF with OOMAO (and the nominal parameters defined above) with a Moffat fit of this PSF. We define a criterion on the PSF profile with respect to a given reference for a given focal plane area (s)  as follows: (7) where REF(x, y) stands for the reference PSF (considered as the true one). This error parameter is used to evaluate the accuracy with which a Moffat can actually fit GLAO PSFs. Let us first focus on the Moffat description of the PSF. In that case, Err <PS F>,Mo f f at,2 * FWH M (as defined in Eq 7) is equal to 1.0, 1.1, and 2% for imaging wavelengths of 500, 700, and 900 nm, respectively. This description of a GLAO PSF is therefore extremely accurate. The redder the wavelength, the better the correction, and therefore the more structured the PSF. This means that the model errors will be greater for the larger wavelengths than for smaller ones. Nevertheless an extensive analysis of the GLAO PSF in various atmospheric conditions and for the whole MUSE WFM spectral range shows that a Moffat fit always gives better results than a few percent which fully validates our choice of a Moffat description for the GLAO PSF. Let us now focus on the FoV evolution of the PSF. As mentioned above, we simulated 47 regularly spaced PSFs in a 1x1 arcmin 2 FoV with our OOMAO simulator. For each PSF, we fitted a Moffat function and we can therefore analyse the evolution of the Moffat parameters (FWHM and β). The FWHM rms error in the FoV is smaller than 20 mas and the value of β is smaller than 0.1. From the previous simulation results we can consider a single PSF and apply it to the whole FoV. We note that, for further development, a more complex PSF model could be investigated. For example, R. Fetick ) recently proposed a new PSF model for AO-corrected applications. This latter model relies on nine parameters and allows the user to fit an AO-corrected PSF both in its corrected area and its wings extremely accurately.

Point-spread-function reconstruction for MUSE Wide Field Mode
This section is dedicated to the presentation of the PSFR algorithm and its performance analysis on simulated data. Using the output from the previous section, we now focus on the three Moffat parameters (α x , α y and β) for each wavelength. One single PSF (resulting from the average of nine PSFs evenly distributed across the FoV) is estimated per wavelength bin. The whole PSFR process is summarised in Figure 5.

The ground-layer adaptive-optics PSFR algorithm
The starting point of the PSFR algorithm is to consider a Fourier basis to describe the whole problematic. Here, it is assumed that everything (phase propagation, WFS measurements, DM commands) is linear and spatially shift-invariant. Hence, all the usual operators are diagonal with respect to spatial frequencies and simply act as spatial filters in the Fourier domain. It follows that each equation can be written frequency by frequency Neichel et al. (2009). The main advantage of the Fourier basis is its diagonal aspect in the frequency domain. It follows that any reconstruction algorithm may be derived and evaluated one Fourier component at a time. In addition, second-order statistics of the residual phase and long-exposure PSF can be evaluated directly without the need for iterations. By avoiding the convergence problem, simulation times are cut down by orders of magnitude at VLT scales.
where PS D φ is the residual phase power spectral density after GLAO correction, PS F T el is the telescope PSF defined by the telescope pupil, and PS F MUS E includes pixel effects and is defined as a centred Gaussian function with a FWHM of 0.2 ′′ . The main limitation of the Fourier approach is that apertureedge effects and boundary conditions that cannot be represented by shift-invariant spatial filters are neglected. Hence, the Fourier modelling only applies to the idealised case of an infinite aperture system, and all effects of incomplete beam overlap in the upper atmospheric layers are neglected. However, in the frame of MUSE WFM, the size of the telescope aperture is large enough (with respect to the sub-aperture diameter and to r 0 ) to satisfy this assumption. Moreover, the GLAO system and its simple averaging process is very well adapted to the Fourier representation: it allows the user to simply and directly focus on the dominant error terms in the error budget (fitting and high-altitudelayer contributions).
Using Equation 8, PSFs are computed for each wavelength at nine positions in the FoV and then averaged. From the averaged PSF, a 2D Moffat fit is performed using a classical Least Square algorithm and the three main Moffat parameters α x (λ i ),α y (λ i ), and β(λ i ) are stored for each λ i bin (for MUSE WFM, the number of bins is equal to 3000).

Sensitivity analysis
A comprehensive analysis of the PSFR algorithm has been provided in close interaction with ESO and MUSE teams during the development process. Here we present a very small subset of the full analysis in order to illustrate the main conclusions. From the GLAO error budget presented in Section 2.1, the dominant term from the performance point of view is the contribution of the uncorrected high-altitude layers (σ 2 High Layers contrib ). This term depends on the three atmospheric parameters only: r 0 @0.5µm (or the seeing value s = 0.1/r 0 in arcsec); ground layer fraction (GLF) -it is worth noting that the combination of GLF and r 0 @0.5µm gives the contribution of the uncorrected high turbulent layers -and L 0 (in m) which is the outer scale of the turbulence.
The combination of r 0 and 1-GLF (fraction of high-layer turbulence) gives the contribution of the uncorrected layers. In the following we study the impact of incorrect values for these three parameters on the PSFR performance. As proposed above, the system PSF is described by a Moffat function and we study the impact of the incorrect parameter values on both the Moffat FWHM and β.
Let us first focus on the seeing. Three cases are considered, assuming a 70% GLF seeing and a 16m L0 (both extracted from typical/median values observed at Paranal): an optimistic case of 0.4" seeing condition; a typical case of 0.8" seeing condition; and a pessimistic case of 1.2" seeing. Figures 6 shows the evolution of the error on FWHM and on β as a function of an error on the seeing estimation. Firstly, we can clearly see that the error on the FWHM estimation is strongly correlated with the error on the seeing with almost a one-to-one relationship. This can be easily understood in the case of partial correction ('typical' and 'pessimistic' cases). In that case the high uncorrected turbulence is responsible for the broadening of the PSF and the PSF FWHM will be directly linked to the high altitude seeing which is given by the following relationship: seeing HL = (1 − GLF) 3/5 * seeing.
As shown previously, the GLAO correction does not only affect the FWHM of the PSF but also its shape (especially far from the optical axis). This shape is represented by the β parameter in the case of a Moffat. Looking at β we can see that this parameter is significantly less affected by an error on the seeing parameter except for the optimistic case (when most of the turbulence is located near the ground). In that case, the correction becomes quite efficient and the shape of the PSF starts to change significantly and therefore the β parameter starts to play a greater role in the overall PSF description. Let us now consider an error on the GLF estimation. Here, again, three cases are considered (assuming a 0.8" seeing and a 16m L 0 ): an optimistic case where 90% of the turbulence is located near the ground and is therefore corrected by the GLAO system. In that case the AO performance becomes quite important and a diffraction-limited core appears; a typical case where 70% of the turbulence is located near the ground. In that case   GLAO provides a significant reduction of the PSF FWHM without achieving the diffraction limit. This case is meant to represent the typical performance expected with the GLAO system for MUSE-WFM; a pessimistic case where a large amount (50%) of turbulence remains uncorrected (in the high-altitude layers). Figure 7 shows the impact of an estimation error of the GLF (and thus of the high-layer uncorrected fraction of the turbulence) on the PSFR accuracy (looking at the Moffat parameters). Here again a linear behaviour between the error on the GLF and the estimated FWHM is found. It is worth noting that even though the relative error on FWHM and β is relatively high for a small estimation error in the high GLF fraction case, the absolute values remain reasonable (see Figure 8).
Errors on FWHM smaller than 50 mas are found for a GLF mis-estimation of typically ±10%. The "good case" (90 % of the turbulence near the ground) is worth to be analysed. In that case, the GLAO system provides a very good correction and PSF are close to be diffraction limited. In that regime, the PSF shape becomes more complex and the impact of inaccurate atmospheric parameters has a more significant impact on the PSFR accuracy. Finally let us focus on the last important atmospheric parameter (especially for a large aperture telescope), the outer scale L 0 . In this case, both seeing and GLF are fixed to their median values (0.8" and 70%). Four L 0 values (8, 16, 24 and 32m) are considered in the simulation in order to span the wide possible range of outer-scale fluctuations. Figure 9 presents the main results obtained for the various real outer-scale values as a function of outer-scale input in the PSFR algorithm. Here, we see that the outer scale estimation is clearly not critical as soon as the real atmospheric outer scale is larger than typically 16m (i.e. two times the telescope diameter). Below this limit, an exact measurement of L 0 becomes important. Fortunately, small L 0 are rare (from Paranal measurements). More importantly, when the outer scale is small, its signature on the WFS signal becomes relatively strong. Therefore, its estimation from RTC data should be accurate enough assuming that the atmospheric parameter measurement from the RTC data process is properly calibrated and validated. The details of these measurements and calibration processes are reported below. Good case (90%), median case (70%) and bad case (50%). the seeing value is equal to 0.8" (which roughly corresponds to the median values measured on 1 years of MUSE-WFM+GLAO operations) and L 0 is equal to 25m. Fig. 9. Relative error on FWHM and β (in %) as a function of the bias on L 0 estimation. 4 cases of real L 0 inputs are considered: 8,16, 24 and 32m. the seeing value is equal to 0.8" and GLF is equal to 70% (which roughly corresponds to the median values measured on 1 years of MUSE-WFM+GLAO operations).

On-sky data
The previous sections give an overview of the PSFR algorithm and of its performance measured on simulated data. This algorithm has been implemented in Python (see Section 6) and has been tested and validated on real MUSE-WFM data acquired during commissioning, science verification, and the early operation periods. The available data can be split into two main categories: -The RTC (also known as SPARTA) data. The instantaneous WFS measurement and DSM command are used to compute statistics on which turbulence models are fitted (Fusco et al. 2004b,a;Fedrigo et al. 2006). For each LGS-WFS signal atmospheric parameters (r 0 , L 0 , GLF wind speed) as well as WFS information (e.g. WFS noise) are saved every 30 seconds. An example of a statistical analysis obtained from these data is plotted in Figure 10 and described in Section 4.1. -The MUSE 3D images (of Global Clusters) themselves. MUSE data will be post-processed and Moffat functions will be fitted on them. This fit will produce the 'reference values' for our on-sky tests. A detailed description of the MUSE data and their processing is provided in Section 4.2.

Adaptive optics telemetry and environmental data
The SPARTA RTC system does continuously (every 30 seconds) provide information on the main atmospheric parameters derived from the LGS WFS and DSM command recorded in closed loop after a pseudo open-loop reconstruction. In the following, we use 392 data sets associated to observations of Globular Clusters performed by MUSE during the periods of October 1, 2017, to August 31, 2018; we kept 355 of them. We discarded data with (i) computation issues (aberrant values) and (ii) very large seeing (> 1.5"). From the remaining RTC data, turbulence parameters are estimated (in the LOS) and a statistical analysis is provided (see Figure 10). The median values of seeing (0.83"), L 0 (16 m) and GLF (72%) are fully compatible with the common Paranal values now recorded for more than 20 years. The seeing proba- bility density function (PDF) has a typical Poisson distribution shape. The L 0 PDF is almost symmetric around its median value with a tight FWHM of only ±3m. Only a very small percentage of the measurements exhibit an outer scale smaller than 8m and the validity of those data points remains questionable. Finally, the GLF PDF is more structured, with some kind of bi-modal shape, and a significant part of the occurrences are found in the 80-90% domain. This aspect could be important in the following. The combination of this large GLF with small seeing values should lead to very good GLAO performance and produce neardiffraction limit images (at least at the AO system focal plan output before entering the MUSE spectrograph). This has two main implications: (i) a higher sensitivity to MUSE internal defects (see Section 5.2) and (ii) a more complex final PSF shape than that coming from the Moffat assumption. This latter effect will probably be one of the main limitations of the current method.

MUSE Wide Field Mode 3D data
The data were obtained within the MUSE globular cluster survey (Kamann et al. 2018a), which is carried out as part of the MUSE guaranteed time observations. The survey targets the central regions of Galactic globular clusters with a series of relatively short exposures. In order to detect variable stars, the observations of each cluster are split into different epochs, with time lags of hours to months between them. Each individual epoch is split into three exposures, in between which derotator offsets

Point-spread-function fit on MUSE data cube
We performed a fit of the PSF on the MUSE data cube using Pam-pelMuse (Kamann 2018;Kamann et al. 2013), a software package designed for the analysis of integral field data of crowded stellar fields such as globular clusters. PampelMuse uses a reference catalogue containing the world coordinates and magnitudes of the sources in the observed field as input and first identifies the subset of available sources that can be resolved from the integral field data. In a subsequent step, PampelMuse determines the coordinate transformation from the input catalogue to the integral field data as well as the PSF as a function of wavelength. This information is finally used to optimally extract the spectra of the resolved sources. The MUSE data considered in this study were analysed using an analytical Moffat profile as PSF. Both the width of the Moffat (parametrised by the FWHM) and its kurtosis (parametrised by the parameter β) were optimised during the analysis and were allowed to change with wavelength. As mentioned above, a standard observation of a globular cluster consists of three exposures with derotator offsets. By default, the exposures are combined before the analysis, in order to homogenise the image quality across the FOV. However, for this project we analysed the individual exposures, which allows for a more direct comparison with the atmospheric parameters gathered during the observations. As the resampling into a data cube can produce artefacts if only a single exposure is used, Pampel-Muse has been modified to work on pixel tables, an intermediate data format used by the MUSE pipeline that does not require resampling (see Weilbacher et al. 2014).

Results
In order to investigate the quality of the PSF fits, we proceeded as follows. When analysing an exposure, PampelMuse selects a number of bright and reasonably isolated stars that are used to optimise the PSF model. The optimisation is done iteratively. The contributions of nearby stars that could potentially disturb the fits are subtracted using an initial PSF model, after which the model is refined by fitting single PSF profiles to the selected stars. The refined model is then used to improve the subtraction of the nearby stars. The steps are repeated until the fluxes of the PSF stars have converged.
After convergence, we extracted radial profiles of the PSF stars and compared them to the radial profiles of the models. We note that before extracting the profiles from the integral field data, we again subtracted the contributions of nearby stars. By subtracting the model profile from the measured one and dividing the result by the measured profile, we determined the relative residuals for each PSF fit. Those were averaged for the 50% brightest PSF stars. Finally, we measured the RMS deviation from zero of the mean relative residuals within the central 2 ′′ . This value, which is shown as a function of β and the FWHM of the fitted Moffat PSF in Fig. 11, serves as our criterion for the agreement between our PSF model and the actual MUSE WFM-AO PSF. It can be understood as the typical residual flux in a pixel after subtraction of a star relative to its recorded flux. The results depicted in Fig. 11 show that the residuals of the PSF fits are typically < 10%, although some cases exist where the  Moffat profiles seem to provide a less accurate fit to the actual PSF. While no obvious trend with the fitted values of β is visible, there is an anti-correlation between the strength of the residuals and the value of the fitted FWHM. While the fit residuals are typically < 5% for observations with FWHM > 0.6 ′′ , stronger residuals are observed for smaller FWHM values. We attribute this behaviour to the PSF becoming critically sampled. The spatial sampling of MUSE in the WFM is 0.2 ′′ , meaning that a PSF with a FWHM of 0.4 ′′ will be approximately Nyquist sampled. Hence, as the width of the PSF approaches this limit, it becomes increasingly difficult to recover its true shape. A direct consequence of this observation is that one has to expect larger PSF residuals for data obtained under better conditions.

On-sky performance of the PSF reconstruction
Using both the PSFR estimate computed with the RTC data and the associated results of the PSF fit on the Globular Cluster images, we can now test and assess the performance of our algorithm on real on-sky data.

Point-spread-function reconstruction and first comparison with MUSE data
The Fourier algorithm provides us with a GLAO PSF but does not account for any instrumental defects. The MUSE image quality is mostly dominated by its 0 ′′ . 2 sampling. Measurements performed during the 'Preliminary Acceptance in Europe' indicate an image quality (FWHM) of between 0.20 ′′ . and 0.27 ′′ . , depending of the channel, and over the full wavelength range. Without any additional available information, the most straightforward way to include that instrumental defects is to convolve the GLAO PSF with a Gaussian function of 0.2 ′′ . FWHM. Even though a very good correlation (defined as a classical Pearson correlation coefficient) is visible in Figure 12 (more than 90% for FWHM and more than 70% for β), a bias characterised by an underestimation of the FWHM is clearly visible with a typical value of 0.25 and 0.3", respectively.

Non-atmospheric part to the PSF
The bias observed in the previous section is mainly attributed to the MUSE internal PSF. The very nature of the instrument makes a precise measurement of the internal PSF, both in the FoV and for each wavelength channel, very challenging (if not impossible). No such measurement was available with a sufficient spatial and spectral resolution. To deal with this specific issue we decided to measure the overall MUSE internal PSF (including both the optics and the detector) using on-sky data. To do so, the following multi-step process was applied: -Identification of a subset of data points (among the 355 available). We chose the best data points of the data set (those for which, at the reddest wavelength, the PSF FWHM computed on the MUSE images is better than 0.4 arcsec). This corresponds to 95 of the 392 data points, that is typically 27% of the data.
This process allows us to obtain an associated MUSE internal PSF characterised by its FWHM and β parameter for each MUSE wavelength channel. This PSF includes both detectors and optical defects. The final results are plotted in Figure 13. The results are fully compatible with the MUSE original specification and with the very few and incomplete MUSE PSF internal measurements made in the laboratory during the Assembly Integration and Test (AIT) period. The MUSE internal PSF (detector included) goes from 0.35" to 0.3" (in the reddest part of the instrument spectrum). Assuming a 0.2" detector FWHM, this corresponds to a full optical error budget of between 0.25 and 0.2". We note that we assumed here that the MUSE internal aberrations were, are, and will be fully static temporally speaking. By design of the instrument and because the instrument lays on the VLT Nasmyth platform we strongly believe that this assumption will remain correct at the level of accuracy required for the MUSE PSF reconstruction algorithm (i.e. that any temporally variable internal aberrations will not affect the internal PSF by more than a few tens of mas in terms of FWHM). Considering that level of GLAO correction and the level of required accuracy on the reconstructed PSF, an error on the static aber- ration of a few tens (up to a few hundreds) of nanometers will be completely negligible compared to the atmospheric contribution. This is far larger than any expected temporal evolution of the instrumental aberrations. The comparison between our measurements and some partial (for only a very limited number of the MUSE channels) internal data acquired during the AIT stage of the instrument seems to fully confirm this hypothesis. In any case, a follow up of the internal aberrations with time could be organised using the same procedure proposed here in order to fully validate the hypothesis.
The MUSE PSF is now included in the complete PSFR algorithm in order to obtain the final performance. Although it has been computed using only the best available data, it is now applied to all the data, assuming that this MUSE internal PSF remains constant during the lifetime of the project (or at least between two re-calibration process).

Final performance
Let us now use the MUSE internal PSF in the full PSFR process and re-process all the data (the 355 available) with the final version of the algorithm. The results are plotted in Figure 14. Results can be compared to Figure 12. The correlation of both FWHM and β remains identical but the bias has completely disappeared for ALL the processed data, showing the pertinence of the MUSE internal PSF for the whole set of data. In order to obtain more quantitative results, we propose in the following to compute an error metric between the final computed parameters using the full PSFR process (including the MUSE internal PSF) and the measured parameters obtained on the Globular Cluster images: (12) Figure 14 shows the PDF of the err p for both FWHM and β for several wavelength ranges. It also shows the PDF cumulative function in each case. From these various plots, we can extract the final on-sky performance of our PSFR algorithm: -For FWHM, the bias is 10mas and the 1 and 2 σ dispersion are respectively 60 mas and 160 mas (which means that the error is smaller than 60 mas in 68% of the cases and smaller than 160 mas in 95.4% of the cases). -For β, the bias is -0.1 and the 1 and 2 σ dispersions are respectively 0.26 and 0.6 (which means that the error is smaller than 0.26 in 68% of the cases and smaller than 0.6 in 95.4% of the cases).
These results, combined with the simulation analysis, fully demonstrate (with respect to our initial scientific requirements) the accuracy and reliability of the PSFR algorithm. This study validates the proposed strategy and allow us to pass to the next level of the project: the final implementation in the MUSE-WFM pipeline and the use of PSFR for scientific observations and final astrophysical data processing.

Implementation
The PSF reconstruction algorithm is implemented as a Python package (muse_psfr), and its source code is available on GitHub 1 .
The algorithm requires three values provided by SPARTA (the AOF Real Time Computer Fedrigo et al. 2006): the seeing, the ground layer fraction (GLF), and the outer-scale (L 0 ). These values can be provided directly as command-line arguments: It is also possible to provide a raw MUSE file. Since the GLAO commissioning, the MUSE raw files contain a FITS table (SPARTA_ATM_DATA) containing the atmospheric turbulence profile estimated by SPARTA. This  The last option is to use the Python API directly, which gives access to more parameters: -Number of reconstructed wavelengths: To reduce computation time, the muse-psfr command reconstructs the PSF at three wavelengths: 500, 700, and 900 nm. But it is possible to reconstruct the PSF at any wavelength, with the compute_psf_from_sparta function. This function reconstructs by default for 35 wavelengths between 490nm and 930nm (which can be specified with the lmin, lmax, and nl parameters) -Number of reconstructed directions: Since the spatial variation is negligible over the MUSE FOV, the reconstruction is done by default only at the centre of the FOV. This can be changed in compute_psf_from_sparta with the npsflin parameter.
The documentation 2 gives more information about the Python API and the various parameters.

Example of application: MUSE deep field
We used the algorithm to estimate the PSF for the MUSE eXtreme Deep fields (MXDF) obervations performed with MUSE in the area of the Hubble Ultra Deep Field (Beckwith et al. 2006). The aim of the project (Bacon et al in prep) is to perform the deepest ever spectroscopic deep field by accumulating more than 100 hours of integration in a single MUSE field. The observations were performed between August 2018 and February 2019 with the GLAO mode. The field location was selected to be in the deepest region of the UDF and to have a tip/tilt star bright enough to ensure a good GLAO correction, plus a fainter star in the outskirts of the MUSE FOV for the slow guiding system. These two requirements allow us to achieve the best possible spatial resolution for the given atmospheric conditions. However, given the poor star density at this location, it was not possible to simultaneously have a PSF star in the FOV and thus an alternative way to estimate the PSF was required.
A total of 377 exposures with 25 mn integration time was obtained. As shown in Fig. 16, exposures were taken in a variety of seeing and ground-layer conditions. For each exposure, we compute with muse-psfr a polynomial approximation of the Moffat parameters FWHM(λ) and β(λ).
We also used the imphot method to estimate the PSF. The method uses the high-resolution HST broad-band images with the corresponding broad-band MUSE reconstructed images to derive the convolution Moffat kernel which minimises the difference between the two images. The method is described in de-  tail in Bacon et al. (2017). We use two HST broad-band images in the F606W and F775W filters which cover the MUSE wavelength range. The corresponding muse-psfr FWHM value at the filter central wavelength is shown for comparison in Fig. 17. The two methods are in good agreement with a scatter of 0.06 arcsec rms for both filters. We note that a systematic mean offset of 0.06 arcsec is measured between the two methods, the imphot method giving higher FWHM than muse-psfr when the PSF is small. This bias is most likely due to the way the sampling is taken into account but it is difficult to come to any conclusions on this matter without an independent ground truth measurement.

Conclusions
We present a simple, efficient, and fully operational (from the astrophysical image analysis of MUSE WFM images) PSFreconstruction algorithm based on a Fourier analysis of the GLAO residual phase statistics completed by dedicated information and measurements concerning the instrument itself. A detailed analysis of the GLAO error budget has allowed us to both simplify and optimise the algorithm. It has been thoroughly and successfully tested with respect to complete End-to-End simulations. A sensitivity analysis allowed us to determine the required accuracy in terms of input parameters. It is shown that precise knowledge (typically a few percent accuracy) on seeing, ground layer, and outer-scale values is required to fulfill the astronomer requirements. The code was then tested with respect to real onsky data obtained during the commissioning and the science verification of the coupling of MUSE-WFM and the AOF. Gathering almost 400 independent observations of globular clusters, PSF reconstruction was performed for each set of data and a statistical analysis of the results was performed. Using a subsample of the data (only those obtained under the best observation conditions) it has been possible to estimate the MUSE internal PSF (at each wavelength). After integration of the MUSE internal PSF into the algorithm, we demonstrated that it is now capable of reconstructing the critical parameters of a PSF (represented by the FWHM and the kurtosis parameters of the Moffat function) with the accuracy required by astronomers in 90% of the observing cases. More precisely, we obtain an error on the PSF FWHM of smaller than 60 mas (less than one-third of a pixel) for 68% of the cases and 160 mas (smaller than the pixel) for 95.4% of the cases. Similarly, for the β parameter, an error smaller than 0.26 is obtained for 68% of the cases and smaller than 0.6 for 95.4%. After this successful validation, the algorithm was implemented as a python package and can be now used routinely with the MUSE 3D data. A first example of application is presented here for the MUSE deep field observations. It is now a fully operational tool available for all users of MUSE.