Issue 
A&A
Volume 635, March 2020



Article Number  A208  
Number of page(s)  14  
Section  Numerical methods and codes  
DOI  https://doi.org/10.1051/00046361/202037595  
Published online  03 April 2020 
Reconstruction of the groundlayer adaptiveoptics point spread function for MUSE wide field mode observations
^{1}
DOTA, ONERA, Université Paris Saclay, 91123 Palaiseau, France
email: thierry.fusco@onera.fr
^{2}
Aix Marseille Univ, CNRS, LAM, Laboratoire d’Astrophysique de Marseille, Marseille, France
^{3}
Univ. Lyon, Univ. Lyon1, ENS de Lyon, CNRS, Centre de Recherche Astrophysique de Lyon UMR5574, 69230 SaintGenisLaval, France
^{4}
Astrophysics Research Institute, Liverpool John Moores University, 146 Brownlow Hill, Liverpool L3 5RF, UK
^{5}
Gemini Observatory/NSF’s OIR Lab, Casilla 603, La Serena, Chile
^{6}
European Southern Observatory, KarlSchwarzschildStr. 2, 85748 Garching bei Muenchen, Germany
Received:
28
January
2020
Accepted:
24
February
2020
Context. Here we describe a simple, efficient, and most importantly fully operational pointspreadfunction (PSF)reconstruction approach for laserassisted ground layer adaptive optics (GLAO) in the frame of the Multi Unit Spectroscopic Explorer (MUSE) wide field mode.
Aims. 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 PSFreconstruction (PSFR) algorithm and test it both in simulations and using onsky data.
Methods. 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 endtoend simulations. Sensitivity to the main atmospheric and AO system parameters was analysed and the code was reoptimised to account for the sensitivity found. Finally, the optimised algorithm was tested and commissioned using more than one year of onsky MUSE data.
Results. We demonstrate with an onsky 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.
Conclusions. 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 postprocessing activity which requires knowledge of the PSF.
Key words: instrumentation: adaptive optics / methods: data analysis / methods: observational / techniques: image processing / techniques: high angular resolution / telescopes
© T. Fusco et al. 2020
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1. 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 groundbased 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 (BeltramoMartin et al. 2019; Fétick et al. 2019b).
In natural seeing observations, the PSF full width at half maximum (FWHM) is often used to quantify the achieved image quality. Most modern groundbased telescopes are equipped with a seeing monitor which provides a realtime 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 outerscale 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 multiGaussian function (Bendinelli et al. 1987; Trujillo et al. 2001; InfanteSainz et al. 2019). However, in some cases there are no bright point sources in the field of view (FoV) 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 deepfield observations like the Hubble Ultra Deep Field (UDF, Bacon et al. 2017).
The main challenges for extragalactic observations are twofold. Firstly, extragalactic 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. (2015), using stateoftheart morphodynamical 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 signaltonoise 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. 2018a,b). Even when a bright and isolated star is present in the FoV, the method described above assumes a uniform PSF over the FoV which is not always the case. We note that for seeinglimited observations with a limited FoV, the atmospheric and telescope PSF can indeed be considered as uniform with the FoV, but this is not generally the case for the instrument. With the generalisation of advanced adaptive optics (AO) systems, modern groundbased 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 wavefront sensing telemetry information and good knowledge of the system, it is theoretically possible to predict the PSF, even without a point source within the FoV (Véran et al. 1997).
Although PSFreconstruction algorithms have been in existence for a long time (Véran et al. 1997; Gendron et al. 2006; Gilles et al. 2012; Ragland et al. 2016; BeltramoMartin 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 PSFreconstruction 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 HAWKI (Pirard et al. 2004) and MUSE (Bacon et al. 2004, 2014) instruments.
With the regular use of MUSE ground layer AO (GLAO) mode, the number of nonAO 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 groundlayer fraction must be taken into account.
Here, we present a PSFreconstruction 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 WFM, we describe the AOF module that allows users to correct for the groundlayer 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 MUSEWFM and AOFGLAO correction. A sensitivity analysis using EndtoEnd simulations is provided and some algorithm parameters are then adjusted accordingly. After demonstrating the algorithm performance on simulated (and thus wellmastered) data, it is applied to a real MUSE onsky 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 PSFreconstruction (PSFR) approach. The final implementation in the MUSE pipeline is detailed and a first astrophysical application to the MUSE UDF 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).
2. MUSE wide field mode
MUSE is the ESO VLT secondgeneration widefield integral field spectrograph operating in the visible (Bacon et al. 2014), covering a simultaneous spectral range of 480–930 nm with a spectral resolution of ∼3000. Its WFM offers a FoV of 1 arcmin^{2}, sampled at 0$\stackrel{\u2033}{.}$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.
2.1. GALACSIGLAO 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 20Watt laser guide stars, and two wavefront sensing units (GALACSI and GRAAL) at each Nasmyth Platform, feeding MUSE and HawkI, 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 1 arcmin^{2} FoV of MUSE (Kolb et al. 2017). The system is robust enough to now be the “normal” mode of operation of MUSE WFM. The MUSE WFM imagequality requirements, mainly driven by the instrument sampling and FoV, have led to a very specific correction stage which is focused on groundlayer 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 wavefront 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 nonAO 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 loworder 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 GLAOcorrected 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 nonsymmetrical 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.
$$\begin{array}{c}\hfill FWH{M}_{\mathrm{Final}}=\sqrt{FWH{M}_{\mathrm{Tel}}^{2}+FWH{M}_{\mathrm{Atm},\mathrm{GLAO}}^{2}+FWH{M}_{\mathrm{MUSE}}^{2}},\end{array}$$(1)
FWHM_{Tel} stands for the telescope diffraction and any defect related to its aberrations that will not be corrected by the AO stage (hightemporalfrequency 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_{Tel} 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 ${\sigma}_{\text{Atm},\text{GLAO}}^{2}$ (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:
$${\sigma}_{\text{Atm},\text{GLAO}}^{2}={\sigma}_{\text{High Order modes}}^{2}+{\sigma}_{\text{Tip tilt}}^{2},$$(2)
$$\begin{array}{c}\text{with}\phantom{\rule{0.166667em}{0ex}}{\sigma}_{\text{High Order modes}}^{2}={\sigma}_{\text{fitting}}^{2}+{\sigma}_{\text{aliasing}}^{2}+{\sigma}_{\text{High Layers contrib}}^{2}\\ +{\sigma}_{\text{noise}}^{2}+{\sigma}_{\text{tempo}}^{2},\end{array}$$(3)
$$\begin{array}{c}\text{and}\phantom{\rule{0.166667em}{0ex}}{\sigma}_{\text{Tiptilt}}^{2}={\sigma}_{\text{TT},\text{aliasing}}^{2}+{\sigma}_{\text{TT},\text{anisoplanatism}}^{2}\\ +{\sigma}_{\text{TT},\text{noise}}^{2}+{\sigma}_{\text{TT},\text{tempo}}^{2}.\end{array}$$(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 cutoff 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 (${\sigma}_{\text{noise}}^{2}\simeq 0$). As shown in Fig. 1, ${\sigma}_{\text{High Layers contrib}}^{2}$ is typically 100 times larger than the other error terms (aliasing and temporal effects).
Fig. 1. Evolution of ${\sigma}_{\text{High Layers contrib}}^{2}$ and other AO error terms (${\sigma}_{\text{noise}}^{2}$, ${\sigma}_{\text{aliasing}}^{2}$ and ${\sigma}_{\text{tempo}}^{2}$) as a function of seeing. The grey and blue areas correspond to the possible variations of atmospheric parameters (wind speed, ${C}_{n}^{2}$ 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 onsky 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 highlayer 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).
Fig. 2. Tip tilt anisoplanatism for various atmospheric conditions obtained using data gathered on the MUSE RTC after one year of observations. The solid line corresponds to an average profile, and the grey area corresponds to the scattering of more than 400 data points gathered during more than one year by the GALACSI RTC. 
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 groundlayer 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.
2.2. EndtoEnd simulation of the MUSE groundlayer adaptiveoptics system
EndtoEnd (E2E) simulations are carried out with Object–Oriented Matlab Adaptive Optics (OOMAO; Conan et al. 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, wavefront sensor (WFS), DM, and an imager of an AO system. It can simulate NGS and LGS singleconjugate 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.
Fig. 3. MUSE WFM and AOF geometry in the FoV. 
System parameters used for the E2E simulation and the PSFR validation process.
Turbulence parameters used for the E2E simulation and the PSFR validation process.
2.3. 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 pixelwise basis. Although by definition, working on a pixelwise basis removes the need for a modelbased 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 multiwavelength 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 postprocessing 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 GLAOcorrected 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.
$$\begin{array}{c}\hfill M(x,y)={M}_{0}{({\left(\frac{x{m}_{x}}{{\alpha}_{x}}\right)}^{2}+{\left(\frac{y{m}_{y}}{{\alpha}_{y}}\right)}^{2}+1)}^{\beta},\end{array}$$(5)
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,
$$\begin{array}{c}\hfill FWH{M}_{x,y}=2{\alpha}_{x,y}\sqrt{{2}^{\frac{1}{\beta}}1},\end{array}$$(6)
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:
$$\begin{array}{c}\hfill {\mathrm{Err}}_{\mathrm{ref},\mathrm{psf}},s=\sqrt{\frac{\int {\int}_{s/2}^{s/2}{\mathrm{PSF}(x,y)\mathrm{REF}(x,y)}^{2}\mathrm{d}x\mathrm{d}y}{\int {\int}_{s/2}^{s/2}{\mathrm{REF}(x,y)}^{2}\mathrm{d}x\mathrm{d}y}}\ast 100,\end{array}$$(7)
Fig. 4. Comparison of a GLAO PSF profile computed in the MUSE WFM. The black dots show 49 regularly spaced positions in 1 × 1 arcmin^{2}. The solid black line represents the average PSF for the 49 positions and the grey area shows the full dispersion of the PSF over the entire FoV. The red line shows the Moffat fit of the mean PSF. 
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_{⟨PSF⟩,Moffat, 2 * FWHM} (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 1 × 1 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 (Fétick et al. 2019a) recently proposed a new PSF model for AOcorrected applications. This latter model relies on nine parameters and allows the user to fit an AOcorrected PSF both in its corrected area and its wings extremely accurately (see Table 3).
Evolution of PSF key parameters (FWHM and β) in the MUSE WFM FoV.
3. Pointspreadfunction 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 Fig. 5.
Fig. 5. Generic description of the PSFR algorithm. 
3.1. The groundlayer adaptiveoptics 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 shiftinvariant. 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, secondorder statistics of the residual phase and longexposure 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.
$$\begin{array}{cc}& \mathrm{PSF}(x,y)={\mathrm{PSF}}_{\mathrm{Tel}}\star {\mathrm{PSF}}_{\mathrm{GLAO}}\star {\mathrm{PSF}}_{\mathrm{MUSE}}\hfill \\ \hfill & {\mathrm{PSF}}_{\mathrm{GLAO}}(x,y)=F{T}^{1}\{exp(\frac{1}{2}FT\left\{{\mathrm{PSD}}_{\varphi}({f}_{x},{f}_{y},\lambda )\right\})\},\hfill \end{array}$$(8)
where PSD_{ϕ} is the residual phase power spectral density after GLAO correction, PSF_{Tel} is the telescope PSF defined by the telescope pupil, and PSF_{MUSE} includes pixel effects and is defined as a centred Gaussian function with a FWHM of 0$\stackrel{\u2033}{.}$2.
The main limitation of the Fourier approach is that apertureedge effects and boundary conditions that cannot be represented by shiftinvariant 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 subaperture 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 highaltitudelayer contributions).
Using Eq. (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).
3.2. 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 Sect. 2.1, the dominant term from the performance point of view is the contribution of the uncorrected highaltitude layers (${\sigma}_{\text{High Layers contrib}}^{2}$). This term depends on the three atmospheric parameters only: r_{0} at 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} at 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 1GLF (fraction of highlayer 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 16 m 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. Figure 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 onetoone 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:
$$\begin{array}{c}\hfill {\mathrm{seeing}}_{\mathit{HL}}={(1\mathrm{GLF})}^{3/5}\ast \mathrm{seeing}.\end{array}$$(9)
Fig. 6. Relative error on FWHM and β (in %) as a function of the bias on seeing estimation. 3 cases of real seeing inputs are considered: good case (0.4″), median case (0.8″) and bad case (1.2″). The GLF value is equal to 70% (which roughly corresponds to the median values measured on 1 years of MUSEWFM+GLAO operations) and L_{0} is equal to 25 m. 
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 16 m 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 diffractionlimited 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 MUSEWFM; a pessimistic case where a large amount (50%) of turbulence remains uncorrected (in the highaltitude layers).
Figure 7 shows the impact of an estimation error of the GLF (and thus of the highlayer 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 Fig. 8).
Fig. 7. Relative error on FWHM and β (in %) as a function of the bias on GLF estimation. 3 cases of real GLF inputs are considered: 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 MUSEWFM+GLAO operations) and L_{0} is equal to 25 m. 
Fig. 8. Absolute error on FWHM (in arcsec) and β (a.u.) as a function of the bias on GLF estimation. 3 cases of real GLF inputs are considered: 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 MUSEWFM+GLAO operations) and L_{0} is equal to 25 m. 
Errors on FWHM smaller than 50 mas are found for a GLF misestimation 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 32 m) are considered in the simulation in order to span the wide possible range of outerscale fluctuations. Figure 9 presents the main results obtained for the various real outerscale values as a function of outerscale 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 16 m (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.
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 32 m. 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 MUSEWFM+GLAO operations). 
4. Onsky 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 Sect. 6) and has been tested and validated on real MUSEWFM 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. 2004a,b; Fedrigo et al. 2006). For each LGSWFS signal atmospheric parameters (r_{0}, L_{0}, GLF wind speed) as well as WFS information (e.g. WFS noise) are saved every 30 s. An example of a statistical analysis obtained from these data is plotted in Fig. 10 and described in Sect. 4.1.
Fig. 10. Real Time Computer DATA: median seeing = 0.83″, median L_{0} = 16 m, median GLF = 72%. 
– The MUSE 3D images (of Global Clusters) themselves. MUSE data will be postprocessed and Moffat functions will be fitted on them. This fit will produce the “reference values” for our onsky tests. A detailed description of the MUSE data and their processing is provided in Sect. 4.2.
4.1. Adaptive optics telemetry and environmental data
The SPARTA RTC system does continuously (every 30 s) provide information on the main atmospheric parameters derived from the LGS WFS and DSM command recorded in closed loop after a pseudo openloop 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 Fig. 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 probability 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 ±3 m. 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 bimodal 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 Sect. 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.
4.2. MUSE wide field mode 3D data
The data were obtained within the MUSE globular cluster survey (Kamann et al. 2018b), 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 of 90° are applied. For this work, we consider all the data taken with the WFM AO system between October 1, 2017, and August31, 2018, that is, in observing periods P100 and P101. In total, 413 individual exposures were analysed.
4.3. Pointspreadfunction fit on MUSE data cube
We performed a fit of the PSF on the MUSE data cube using PAMPELMUSE (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, PAMPELMUSE 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).
4.4. 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 WFMAO 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 anticorrelation 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.
Fig. 11. Average relative residuals of the PSF fits in the central 2″ as a function of β (left) and the FWHM (right) of the Moffat profile used to fit the MUSE WFMAO PSF. 
5. Onsky 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 onsky data.
5.1. Pointspreadfunction 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\stackrel{\u2033}{.}2$ sampling. Measurements performed during the “Preliminary Acceptance in Europe” indicate an image quality (FWHM) of between 0$\stackrel{\u2033}{.}$20 and 0$\stackrel{\u2033}{.}$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$\stackrel{\u2033}{.}$2 FWHM. Even though a very good correlation (defined as a classical Pearson correlation coefficient) is visible in Fig. 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.
Fig. 12. Comparison of FWHM (top) and β parameter (bottom) for measured data on MUSE images (extract from Globular Custer data) and PSFR data (computed from atmospheric data obtained using AO telemetry) convolved by a 0.2″ FWHM Gaussian function (to account for MUSE pixel size). The various colors stand for various seeing values domain (black: 0.2–0.4″ – red: 0.4–0.6″ – orange: 0.6–0.8″ – green: 0.8–1.0″ – blue: 1.0–1.5″ 
5.2. Nonatmospheric 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 onsky data. To do so, the following multistep 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.
– On this data subset, for each MUSE wavelength, the GLAO PSFR computed with our algorithm were convolved by a MUSE internal PSF (PSF_{MUSE} also modelled by a Moffat function characterised by its FWHM_{MUSE}(λ) and β_{MUSE}(λ) parameters). Using a classical least square error metric we adjust FWHM_{MUSE}(λ) and β_{MUSE}(λ) in order to globally minimise the quadratic distance between the subset of MUSE onsky data and the associated PSFR convolved with PSF_{MUSE}, that is, for each wavelength λ and each data point of the subset:
$$\begin{array}{cc}& \mathrm{Min}\left\{{\left{\mathrm{PSF}}_{\mathrm{meas},i}{\mathrm{PSFR}}_{i,\lambda}^{FWHM,\beta}\right}^{2}\right\}\hfill \\ \hfill & \mathrm{with}\phantom{\rule{0.166667em}{0ex}}\mathrm{respect}\phantom{\rule{0.166667em}{0ex}}\mathrm{to}\phantom{\rule{0.166667em}{0ex}}FWHM\phantom{\rule{0.166667em}{0ex}}\mathrm{and}\phantom{\rule{0.166667em}{0ex}}\beta ,\hfill \end{array}$$(10)
with
$${\text{PSFR}}_{i,\lambda}^{FWHM,\beta}=$$(11)
$${\text{PSF}}_{\text{Tel}}\star {\text{PSF}}_{\text{GLAO}}({\lambda}_{i})\star M(FWH{M}_{\text{MUSE}}({\lambda}_{i}),{\beta}_{\text{MUSE}}({\lambda}_{i})).$$(12)
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 Fig. 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″.
Fig. 13. Full width at half maximum and β parameter estimated for the MUSE internal PSF. A set of FWHM and β is estimated for each wavelength bin between 480 and 940 nm. The blind area around 589 nm is determined by the notch filter which blocks the laser guide star light in the instrument. 
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 aberration 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 recalibration process).
5.3. Final performance
Let us now use the MUSE internal PSF in the full PSFR process and reprocess all the data (the 355 available) with the final version of the algorithm. The results are plotted in Fig. 14. Results can be compared to Fig. 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:
$$\begin{array}{c}\hfill {\mathrm{err}}_{p}={p}_{\mathrm{measured}}{p}_{\mathrm{predicted}}.\end{array}$$(13)
Figure 15 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 onsky performance of our PSFR algorithm:
Fig. 15. Probablity density function of the error on Moffat paremeters (FWHM and β) for several wavelength ranges. 
– 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 MUSEWFM pipeline and the use of PSFR for scientific observations and final astrophysical data processing.
6. 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 outerscale (L_{0}). These values can be provided directly as commandline arguments:
$ musepsfr nocolor values 1,0.7,25 MUSEPSFR version 1.0rc2 Computing PSF Reconstruction from Sparta data Processing SPARTA table with 1 values, njobs=1 Compute PSF with seeing=1.00 GL=0.70 L0=25.00  LBDA 5000 7000 9000 FWHM 0.85 0.73 0.62 BETA 2.73 2.55 2.23 
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 table contains the values for each laser, with one row every two minutes.
$ musepsfr MUSE.20180813T07:14:11.128.fits.fz MUSEPSFR version 0.31 OB MXDF0100A 20180813T07:39:21.835 Airmass 1.491.35 Computing PSF Reconstruction from Sparta data Processing SPARTA table with 13 values, njobs=1 4/13 : Using only 3 values out of 4 ... 4/13 : seeing=0.57 GL=0.75 L0=18.32 Using three lasers mode 1/13 : Using only 3 values out of 4 ... 1/13 : seeing=0.71 GL=0.68 L0=13.60 Using three lasers mode 6/13 : Using only 3 values out of 4 ... 6/13 : seeing=0.60 GL=0.75 L0=16.47 Using three lasers mode .... OB MXDF0100A 20180813T07:39:21.835 Airmass 1.491.35  LBDA 5000 7000 9000 FWHM 0.57 0.46 0.35 BETA 2.36 1.91 1.64 
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 musepsfr 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 490 nm and 930 nm (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.
7. 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 UDF (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 h 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 groundlayer conditions. For each exposure, we compute with musepsfr a polynomial approximation of the Moffat parameters FWHM(λ) and β(λ).
Fig. 16. Histogram of atmospheric parameters measured by the SPARTA realtime AO controller during the MXDF observations. The solid line displays the median value. 
We also used the imphot method to estimate the PSF. The method uses the highresolution HST broadband images with the corresponding broadband MUSE reconstructed images to derive the convolution Moffat kernel which minimises the difference between the two images. The method is described in detail in Bacon et al. (2017). We use two HST broadband images in the F606W and F775W filters which cover the MUSE wavelength range. The corresponding musepsfr 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 musepsfr 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.
Fig. 17. Comparison between Reconstructed PSF Moffat FWHM (PSFRec) and values derived from comparison with highresolution HST broadband images (IMPHOT) for two filters (F606W and F775W). 
8. 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 EndtoEnd 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 outerscale 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 MUSEWFM 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 onethird 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.
Acknowledgments
This work has been partially supported by the ANRAPPLY program (ANR19CE310011), the European Research Council consolidator grant (ERCCoG646928MultiPop) and by OPTICON, in the Horizon 2020 Framework Program of the European Commission’s (Grant number 730890 – OPTICON).
References
 Andersen, D. R., Stoesz, J., Morris, S., et al. 2006, PASP, 118, 1574 [NASA ADS] [CrossRef] [Google Scholar]
 Arsenault, R., Madec, P. Y., Hubin, N., et al. 2008, in Adaptive Optics Systems, eds. N. Hubin, C. E. Max, & P. L. Wizinowich (SPIE), Int. Soc. Opt. Photonics, 7015, 577 [Google Scholar]
 Bacon, R., Bauer, S. M., Bower, R., et al. 2004, in The SecondGeneration VLT Instrument MUSE: Science Drivers and Instrument Design, eds. A. F. M. Moorwood, & M. Iye, SPIE Conf. Ser., 5492, 1145 [Google Scholar]
 Bacon, R., Vernet, J., Borisova, E., et al. 2014, The Messenger, 157, 13 [NASA ADS] [Google Scholar]
 Bacon, R., Conseil, S., Mary, D., et al. 2017, A&A, 608, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Beckwith, S. V. W., Stiavelli, M., Koekemoer, A. M., et al. 2006, AJ, 132, 1729 [NASA ADS] [CrossRef] [Google Scholar]
 BeltramoMartin, O., Correia, C. M., Ragland, S., et al. 2019, MNRAS, 487, 5450 [NASA ADS] [CrossRef] [Google Scholar]
 Bendinelli, O., Zavatti, F., & Parmeggiani, G. 1987, JApA, 8, 343 [NASA ADS] [Google Scholar]
 Bouché, N., Carfantan, H., Schroetter, I., MichelDansac, L., & Contini, T. 2015, AJ, 150, 92 [NASA ADS] [CrossRef] [Google Scholar]
 Conan, R., & Correia, C. 2014, in Adaptive Optics Systems IV, eds. E. Marchetti, L. M. Close, & J. P. Véran (SPIE), Int. Soc. Opt. Photonics, 9148, 2066 [Google Scholar]
 Damjanov, I., Abraham, R. G., Glazebrook, K., et al. 2011, PASP, 123, 348 [NASA ADS] [CrossRef] [Google Scholar]
 Epinat, B., Amram, P., Balkowski, C., & Marcelin, M. 2010, MNRAS, 401, 2113 [NASA ADS] [CrossRef] [Google Scholar]
 Fedrigo, E., Donaldson, R., Soenke, C., et al. 2006, Proc. SPIE [Google Scholar]
 Fétick, R. J. L., Fusco, T., Neichel, B., et al. 2019a, A&A, 628, A99 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Fétick, R. J. L., Jorda, L., Vernazza, P., et al. 2019b, A&A, 623, A6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Fusco, T., Rousset, G., Rabaud, D., et al. 2004a, Journal of Optics A: Pure Appl. Opt., 6, 585 [Google Scholar]
 Fusco, T., Ageorges, N., Rousset, G., et al. 2004b, in Advancements in Adaptive Optics, eds. D. B. Calia, B. L. Ellerbroek, & R. Ragazzoni (SPIE), Int. Soc. Opt. Photonics, 5490, 118 [Google Scholar]
 Gendron, E., Clénet, Y., Fusco, T., & Rousset, G. 2006, A&A, 457, 359 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gilles, L., Correia, C., Véran, J.P., Wang, L., & Ellerbroek, B. 2012, Appl. Opt., 51, 7443 [NASA ADS] [CrossRef] [Google Scholar]
 InfanteSainz, R., Trujillo, I., & Román, J. 2019, VizieR Online Data Catalog: J/MNRAS/491/5317 [Google Scholar]
 Kamann, S. 2018, Astrophysics Source Code Library [record ascl:1805.021] [Google Scholar]
 Kamann, S., Wisotzki, L., & Roth, M. M. 2013, A&A, 549, A71 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kamann, S., Husser, T.O., Dreizler, S., et al. 2018a, MNRAS, 473, 5591 [NASA ADS] [CrossRef] [Google Scholar]
 Kamann, S., Bastian, N., Husser, T.O., et al. 2018b, MNRAS, 480, 1689 [NASA ADS] [CrossRef] [Google Scholar]
 Kolb, J., Madec, P. Y., Arsenault, R., et al. 2017, https://doi.org/10.26698/AO4ELT5.0038 [Google Scholar]
 Madec, P. Y., Arsenault, R., Kuntschner, H., et al. 2018, in Adaptive Optics Systems VI, eds. L. M. Close, L. Schreiber, & D. Schmidt (SPIE), Int. Soc. Opt. Photonics, 10703, 1 [Google Scholar]
 Moffat, A. F. J. 1969, A&A, 3, 455 [NASA ADS] [Google Scholar]
 Müller Sánchez, F., Davies, R. I., Eisenhauer, F., et al. 2006, A&A, 454, 481 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Neichel, B., Fusco, T., & Conan, J.M. 2009, J. Opt. Soc. Am. A, 26, 219 [NASA ADS] [CrossRef] [Google Scholar]
 Oberti, S., Kolb, J., Madec, P. Y., et al. 2018, in Proc. SPIE, SPIE Conf. Ser., 10703, 107031G [Google Scholar]
 Pirard, J. F., KisslerPatig, M., Moorwood, A., et al. 2004, in HAWKI: A new widefield 1 to 2.5μm imager for the VLT, eds. A. F. M. Moorwood, & I. Iye, SPIE Conf. Ser., 5492, 1763 [Google Scholar]
 Ragland, S., Jolissaint, L., Wizinowich, P., et al. 2016, in Adaptive Optics Systems V, eds. E. Marchetti, L. M. Close, & J. P. Véran (SPIE), Int. Soc. Opt. Photonics, 9909, 573 [Google Scholar]
 Trujillo, I., Aguerri, J. A. L., Cepa, J., & Gutiérrez, C. M. 2001, MNRAS, 328, 977 [NASA ADS] [CrossRef] [Google Scholar]
 Véran, J.P., Rigaut, F., Maître, H., & Rouan, D. 1997, J. Opt. Soc. Am. A, 14, 3057 [NASA ADS] [CrossRef] [Google Scholar]
 Weilbacher, P. M., Streicher, O., Urrutia, T., et al. 2014, in The MUSE Data Reduction Pipeline: Status after Preliminary Acceptance Europe, eds. N. Manset, & P. Forshay, ASP Conf. Ser., 485, 451 [Google Scholar]
 Wright, S. A., Larkin, J. E., Law, D. R., et al. 2009, ApJ, 699, 421 [NASA ADS] [CrossRef] [Google Scholar]
All Tables
Turbulence parameters used for the E2E simulation and the PSFR validation process.
All Figures
Fig. 1. Evolution of ${\sigma}_{\text{High Layers contrib}}^{2}$ and other AO error terms (${\sigma}_{\text{noise}}^{2}$, ${\sigma}_{\text{aliasing}}^{2}$ and ${\sigma}_{\text{tempo}}^{2}$) as a function of seeing. The grey and blue areas correspond to the possible variations of atmospheric parameters (wind speed, ${C}_{n}^{2}$ profiles, outer scales) for the various error items. 

In the text 
Fig. 2. Tip tilt anisoplanatism for various atmospheric conditions obtained using data gathered on the MUSE RTC after one year of observations. The solid line corresponds to an average profile, and the grey area corresponds to the scattering of more than 400 data points gathered during more than one year by the GALACSI RTC. 

In the text 
Fig. 3. MUSE WFM and AOF geometry in the FoV. 

In the text 
Fig. 4. Comparison of a GLAO PSF profile computed in the MUSE WFM. The black dots show 49 regularly spaced positions in 1 × 1 arcmin^{2}. The solid black line represents the average PSF for the 49 positions and the grey area shows the full dispersion of the PSF over the entire FoV. The red line shows the Moffat fit of the mean PSF. 

In the text 
Fig. 5. Generic description of the PSFR algorithm. 

In the text 
Fig. 6. Relative error on FWHM and β (in %) as a function of the bias on seeing estimation. 3 cases of real seeing inputs are considered: good case (0.4″), median case (0.8″) and bad case (1.2″). The GLF value is equal to 70% (which roughly corresponds to the median values measured on 1 years of MUSEWFM+GLAO operations) and L_{0} is equal to 25 m. 

In the text 
Fig. 7. Relative error on FWHM and β (in %) as a function of the bias on GLF estimation. 3 cases of real GLF inputs are considered: 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 MUSEWFM+GLAO operations) and L_{0} is equal to 25 m. 

In the text 
Fig. 8. Absolute error on FWHM (in arcsec) and β (a.u.) as a function of the bias on GLF estimation. 3 cases of real GLF inputs are considered: 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 MUSEWFM+GLAO operations) and L_{0} is equal to 25 m. 

In the text 
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 32 m. 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 MUSEWFM+GLAO operations). 

In the text 
Fig. 10. Real Time Computer DATA: median seeing = 0.83″, median L_{0} = 16 m, median GLF = 72%. 

In the text 
Fig. 11. Average relative residuals of the PSF fits in the central 2″ as a function of β (left) and the FWHM (right) of the Moffat profile used to fit the MUSE WFMAO PSF. 

In the text 
Fig. 12. Comparison of FWHM (top) and β parameter (bottom) for measured data on MUSE images (extract from Globular Custer data) and PSFR data (computed from atmospheric data obtained using AO telemetry) convolved by a 0.2″ FWHM Gaussian function (to account for MUSE pixel size). The various colors stand for various seeing values domain (black: 0.2–0.4″ – red: 0.4–0.6″ – orange: 0.6–0.8″ – green: 0.8–1.0″ – blue: 1.0–1.5″ 

In the text 
Fig. 13. Full width at half maximum and β parameter estimated for the MUSE internal PSF. A set of FWHM and β is estimated for each wavelength bin between 480 and 940 nm. The blind area around 589 nm is determined by the notch filter which blocks the laser guide star light in the instrument. 

In the text 
Fig. 14. As in Fig. 12, but the pixel PSF has been replaced by a full internal MUSE PSF. 

In the text 
Fig. 15. Probablity density function of the error on Moffat paremeters (FWHM and β) for several wavelength ranges. 

In the text 
Fig. 16. Histogram of atmospheric parameters measured by the SPARTA realtime AO controller during the MXDF observations. The solid line displays the median value. 

In the text 
Fig. 17. Comparison between Reconstructed PSF Moffat FWHM (PSFRec) and values derived from comparison with highresolution HST broadband images (IMPHOT) for two filters (F606W and F775W). 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.