1.

Introduction

Next generation workhorse facilities for 8- to 10-m class ground-based telescopes [e.g., enhanced resolution imager and spectrograph (ERIS) and MCAO-assisted visible imager and spectrograph (MAVIS) at Very Large Telescope; System for coronagraphy with High order Adaptive optics from R to K band (SHARK)-visible and SHARK-near-infrared (NIR) at Large Binocular Telescope (LBT); Keck All-sky Precision Adaptive (KAPA)-optics and Keck Planet Finder (KPF) at W. M. Keck Observatory] will improve adaptive optics (AO) assisted imaging and spectroscopy, both in the visible and at NIR wavelengths. The great advances obtained by these AO systems will pave the road for extremely large telescopes [e.g., Giant Magellan Telescope (GMT), European Southern Observatory’s Extremely Large Telescope (ELT), and Thirty Meter Telescope (TMT)] in the future to obtain sharper and deeper images than the ones obtained with big space telescopes (e.g., Hubble Space Telescope and James Webb Space Telescope), exploiting all their ground-breaking imaging capabilities.

In many scientific applications exploiting diffraction limit AO capabilities,^1–6 a serious challenge experienced during data analysis is the precise knowledge of the two-dimensional point spread function (PSF): approaching the diffraction limit of ground-based cornerstone facilities is not an easy task. A severe drawback in this sense is that the PSF of an AO system is a very complex entity, which has a high number of degrees of freedom involved, varies spatially in the field of view (FOV), is not constant in time, and behaves differently at each wavelength. Accurate knowledge of the PSF is mandatory to perform high-quality photometric, spectroscopic, and astrometric analysis of astronomical AO data-sets. For instance, a detailed knowledge of the PSF is an unavoidable ingredient to carry out transformational science cases, e.g., resolved stellar clusters beyond the Milky Way, stellar populations beyond the Local Group up to the Virgo cluster, galaxy size evolution, stellar mass growth across cosmic time, substructures in intermediate-$z$ galaxies, detailed studies of high-redshift galaxies, and the formation of the smallest galaxies and proto-globular clusters at high $z$. In addition, a precise reconstruction of the PSF for astronomical images is fundamental, e.g., for: (i) quality control of scientific frames (especially blank extra-galactic fields); (ii) flux calibration of spectroscopic data (since the PSF depends on the wavelength); (iii) unbiased photometry in crowded fields; (iv) image restoration (e.g., Richardson–Lucy deconvolution); and (v) forward modeling. Currently, there is a growing interest in the development of tools to infer the observed PSF from AO systems. The three main approaches of PSF estimation are (1) based on focal-plane data (e.g., PSF classification, PSF fitting/interpolation, PSF marginal estimation and deconvolution, empirical PSF derivation); (2) based on pupil-plane data [e.g., PSF simulation, PSF reconstruction (PSF-R)]; and (3) hybrid techniques that combine the use of point sources in the observed scientific (co-added) images, models of atmospheric turbulence, and instrumental aberrations, hence optimizing the reconstruction (e.g., adaptive PSF-R, myopic deconvolution, calibrated PSF-R). All these methods have been extensively described in Ref. 7. Focal-plane-based methods by default include the non common path aberrations (NCPA) components, which means that imperfections in the optics of the instruments are automatically included in the reconstruction. But they are also sensitive to the image noise, and they can be considered representative of a specific positions in the field and/or for a particular wavelength. Moreover, they require the presence of at least one non-saturated point-like source in the imaged FOV. Pupil-plane based methods, instead, refer to techniques aimed at deriving the post-AO PSF only from telemetry control loop data, e.g., time series of wavefront sensor (WFS) measurements or applied deformable mirror commands and associated calibration, without any recourse to focal plane data. NCPA or segment differential piston should be taken into account separately, following dedicated independent calibrations. Testing the result of pupil-plane methods can only be carried out via comparison of PSFs in the science data.

PSF-R is a versatile and flexible approach, which can also be the initial step for hybrid or adaptive methods. The PSF-R technique is the only viable possibility when no other solutions are possible, due to several reasons: e.g., there are no suitable point sources in the science frames, the target is faint and complex, as in dense and crowded fields (i.e., the core of globular clusters, where only the PSF cores are detected for the majority of the stars), the bright stars are saturated, the integration times are relatively short, observations have been carried out in poor atmospheric conditions or faint guide stars have been used. Indeed, for both single conjugated adaptive optics (SCAO⁸) and multi conjugated adaptive optics (MCAO^9,10) observing modes, the AO reference star can be out of the instrumental field of view.

Moreover, in laser assisted¹¹ MCAO mode, only faint tip–tilt stars are usually present. PSF-R is particularly useful in spectroscopic mode, where no PSF is available in the focal plane, or in non-sidereal tracking mode, where the stars in the field of view may be blurred along one axis and so cannot be used as a PSF reference, or where the Solar System targets (used also as AO reference) can be spatially extended.

The PSF-R approach is particularly convenient for the extra-galactic targets where the field of view is typically void of bright point sources, necessary to characterize the PSF of the science data using other methods. As an example, simulations of stellar densities with TRILEGAL¹² code indicate that, in a typical extra-galactic field, one can expect on average $\sim 3$ stars per square arcmin in the range $Ks=16-25$ Vega mag. Considering the typical field of view of AO assisted instruments (of the order of $1{\mathrm{arcmin}}^2$), the statistical probability of having at least one available star is less than 35% to 45%. The PSF-R remains the unique way forward.

In this paper, we present the first results of our ongoing project aimed at obtaining a PSF-R tool, i.e., the estimation of the PSF using only AO control-loop telemetry, following the seminal work of Ref. 13. The adopted PSF-R method is described in Ref. 14 and involves several approximations to reconstruct the PSF from the residual phase power spectrum density. Such algorithm was previously tested on simulated data only. Therefore, the present work is the first application of the method to real telescope data and is meant to validate the approach. For this reason, we limited the present analysis to high signal-to-noise ratio (SNR), on-axis observations to reduce the degree of freedom of the problem and avoid degeneracies in multi-parameter space (e.g., when taking into account wider AO reference SNR range, off-axis angles, and seeing values). We investigate the performance of the selected software with archival, SCAO data from the LBT¹⁵. Specifically, the data consist of two distinct set of observations of bright, on-axis point sources, preliminary described in Ref. 16. They both have been taken as part of the single conjugated adaptive optics upgrade for LBT (SOUL¹⁷) validation operations, using LBT Utility Camera in the Infrared (LUCI¹⁸). Moreover, to quantify and compare the level of accuracy obtained, we present an idealized scientific case in the area of applicability of the PSF-R, showing the gain offered by this technique to the analysis of the photometry and morphology of a high-$z$ galaxy.

The paper is organized as follows: the PSF-R algorithm is described in Sec. 2; the data are presented in Sec. 3; in Sec. 4 the performance of the PSF-R algorithm is discussed evaluating also the impact on the measure of the morphological parameters of compact galaxies; and Sec. 5 is dedicated to the summary and conclusions. Finally, an interesting relation linking the encircled energy (EE), the full width at half maximum (FWHM), and Strehl ratio (SR) of the observed PSF is discussed in Appendix A.

2.

PSF-R Method

In this section, we describe the method used to reconstruct the on-axis PSF from AO telemetry data, the one described in Ref. 14, which is by itself an update of the classical algorithm developed by Ref. 13. As detailed in Ref. 14, the algorithm is highly flexible with respect to the WFS, and we performed minor tuning to adapt it to SOUL + LUCI. Specifically, to our knowledge this is the first successful implementation of a PSF-R algorithm to a pyramid WFS AO system. This is especially relevant also in sight of the upcoming next generation AO instruments that will make wide use of this kind of WFS. The core idea of this approach is to compute the optical transfer function (OTF), the Fourier transform of the PSF derived from the pupil mask. This is advantageous because the OTF can subsequently be factorized into several independent components (dependent on wavelength and line-of-sight), which can be either calculated from the wavefront residual phase or using tabulated values derived from appropriate calibrations.

The OTF can be defined as

Following the mathematical analysis of Refs. 13 and 14, ${\mathrm{OTF}}_{\mathrm{AO}}$ is described as

To compute $D_{\parallel }$, one needs the residual phase in the space of the mirror modes. However, we only have an estimate for it, namely, the reconstructed residual phase computed from the WFS data, which is perturbed by noise and aliasing effects. Therefore, since the assumptions from Refs. 13 and 14 of fast frame rate and a least squares reconstructor are fulfilled, we can use the same decomposition as in Ref. 13 to decompose $D_{\parallel }$ further as

Table 1

Structure functions computed for the PSF-R. A description on how they are obtained is also provided.

The last term needed for the computation of the post-AO OTF (${\mathrm{OTF}}_{\mathrm{AO}}$) is ${\mathrm{OTF}}_{\mathrm{tel}}$, which is computed from the calibrated pupil mask and then multiplied by the aforementioned filters following Eq. (4).

$K_{\mathrm{NCPA}}$ in Eq. (1) is computed from the calibrated phase map and depends on the specific instrument configuration and telescope pointing. It therefore varies in the field.

The final PSF is then retrieved by transforming the ${\mathrm{OTF}}_{\mathrm{AO}+\mathrm{NCPA}}$ via inverse Fourier transform into the ${\mathrm{PSF}}_{\mathrm{REC}}$, which is also integrated in wavelength over the pass-band of the selected filter.

We note that focusing only on on-axis point source, the anisoplanatic component can be neglected, thus reducing the degrees of freedom in the PSF-R. Moreover, considering that during the observations described in this work the loop was sampled at frequencies $>500\mathrm{Hz}$, and we expect that non-negligible physical vibrations occur at much slower variations, we can assume that the SOUL + LUCI system already accounted for all these variations. Therefore, we do not need to account for any additional term related to vibrations, as it is already present in the used data.

3.

Data

As the main goal of this work is to apply the selected algorithm for PSF-R to real AO data, we used SCAO near infrared images and associated synchronous telemetry acquired with the SOUL system mounted on the LBT and feeding the LUCI near-infrared camera. This, in turn, is equipped with three optical zoom configurations delivering the same number of distinct resolutions (see Ref. 19 for a detailed description of the instrument set-up). We use the diffraction limited mode, called N30, a $2048\times 2048\mathrm{px}$ configuration with ${30}^{\prime \prime }\times {30}^{\prime \prime }$ field of view ($14.95\mathrm{mas}/\mathrm{px}$).

We used archive data, consisting of two independent datasets. The first set refers to observations of an artificial source (i.e., an illuminated fiber—hereafter “daytime” dataset) with controlled atmospheric condition.²⁰ The second dataset consists of observations of a natural point source (i.e., a bright star—hereafter “nighttime” dataset). The log of the observations is presented in Table 2.

Table 2

Log of the SOUL+LUCI at LBT observations.

3.2.

Reconstructed PSFs

In the following, we provide additional technical details on the application of the adopted PSF-R method in the specific case of SOUL + LUCI data. We only used the telemetry data associated with the scientific observations to reconstruct the PSFs.

The telemetry data provided by the SOUL team include: AO WFS slopes data history, control matrix, interaction matrix, gain vector, and the pupil definition. In particular, we obtained pupil images taken with LUCI in the pupil image mode used for the alignment of the instrument. This imaging mode unfortunately does not include the additional pupil obstruction exclusive of the N30 camera. We derived them from the camera mechanical drawing, to obtain the pupil shape for daytime dataset. For nighttime, additional contributions have been taken into account in defining the pupil shape. In this configuration, in fact, LUCI sees also the secondary mirror mounting (including the LBT swing arm) shadow occurring in the pupil mask. While the N30 camera system is fixed with the LUCI detector, the swing arm projection depends on the orientation of the LUCI camera that is continuously updating to follow the apparent rotation of the sky. The derived pupil shape is used as a binary mask assuming negligible modulation of the amplitude (optical system transmission) over the aperture. Using the LUCI images we verify the pixelscale value, which results 14.95 mas/px for both the H and FeII filters.

In all the considered observations, NCPA between WFS optical train and the LUCI camera are compensated by the AO system making the pyramid loop closing on non-zero reference slopes. This correction, however, leaves residual depending on the pyramid optical gain²¹ and un-calibrated reciprocal flexures of the WFS with respect to LUCI. These residuals may count up to 100  nm¹⁹ of wavefront error. Lacking proper calibrations, we measured these residuals by analyzing the observed PSF: this method, applied on long exposure PSF, has low sensitivity. In the present case we observed residuals (tip–tilt free) accounting for about 40  nm in both daytime and nighttime. The measured phase maps for NCPA are shown in Fig. 1. As a consequence, the small residual values also indicate that the pyramid WFS gain measures registered in the telemetry data are accurate. Even if negligible, these estimated wavefront errors were also applied in the PSF-R algorithm.

Fig. 1

NCPA mean residuals, after AO correction, (a) daytime and (b) nighttime. The pixel values are in nm.

3.4.

Daytime

For this particular dataset, laboratory-like conditions have been used: optical turbulence disturbance has been generated and corrected²² using the adaptive secondary²³ of the LBT, with the telescope and AO system configured for daytime.²⁰ Specifically, a seeing of 1.2 arcsec was simulated with an average wind speed of $15\mathrm{m}/\mathrm{s}$.¹⁹ The simulated point source has an equivalent $R$-band magnitude of 10 mag. Three sets of 30 1 s-images have been obtained in $H$ band (1.5 to $1.8\mu \mathrm{m}$, centered at $1.65\mu \mathrm{m}$), achieving a temporal frequency for the AO correction of 500  Hz. The resulting saving rate for the associated telemetry is $500\mathrm{frames}/\mathrm{s}$. In all the science frames the point source is not saturated, within the linear regime of the detector. The resulting SNR per pixel for the observed PSF goes from 250 at peak to a mean value of 10 at a radial distance of 300 mas (20 px) from PSF centre and reaching mean SNR per pixel$=1$ around 0.9 arcsec (60 px). The observed PSF for daytime is presented in Fig. 2(a). The measured SR—defined as the ratio between the peak values of a selected PSF and of the diffraction limited one—is 0.59 and the FWHM is 2.95 px. Figure 2(b) shows the reconstructed PSF. The match between the two PSFs can be qualitatively inferred from the right panel of the same figure, where the normalized residuals are presented. No evident residual PSF structures are present. A quantitative assessment on the performance of the reconstructed PSF is presented in the following sections.

Fig. 2

Daytime PSFs. (a) observed PSF. (b) Reconstructed PSF. Both are in arbitrary logarithmic color scales. (c) Normalized residuals between observed and reconstructed PSF. The residuals are normalized pixel-per-pixel with the observed PSF. The color scale is linear and the colorbar is on the bottom. There are no evident residual structures of the two PSFs. All the images are $\sim 8{\mathrm{arcsec}}^2$ ($190\mathrm{px}\times 190\mathrm{px}$).

3.5.

Nighttime

The nighttime dataset has been obtained in the narrow-band FeII filter, centered at $1.65\mu \mathrm{m}$, with FWHM $0.02\mu \mathrm{m}$. The use of a narrow-band filter reduces the impact of atmospheric dispersion on data: LUCI is not equipped with an atmospheric dispersion corrector and the associated blur of the PSF is of the order of few pixels, so we properly take into account the effect while reconstructing the PSF. The dataset consists of on-axis observation of a bright star (HD 873, $R=8.57\mathrm{mag}$). The star has been also used as AO reference, thus obtaining the best correction from the AO system along the direction of the target. In this case, 20 images of 0.313  s have been considered. The measured seeing was 1.20 arcsec and the wind speed 3.4  m/s. Using a bright target, the frequency of AO correction could be raised to 1700  Hz, with a resulting saving rate for the associated telemetry of 850  frames/s (decimated by a factor 2). The measured SR is 0.58 and the FWHM 43.95 mas (2.93  px). Also in this case, all the images of the point source are not saturated and within the linear regime of the detector. The final SNR per pixel of the observed PSF is 50 at the peak, decreasing to a mean value of 10 at 300 mas. For this dataset, a mean SNR per pixel = 1 is reached around 1.3 arcsec (90 px) from PSF centre. The difference between the SNR radial trend of the two dataset is due to the different bandwidth of the filters used for each dataset: a broadband one for daytime and a narrowband one for nighttime. The variation in environmental conditions during the observations is also a contributing factor. The observed PSF is shown in Fig. 3(a), the reconstructed one in the central panel. In the reconstructed PSF it can be noted a feeble ring-like structure, with a radius of $\sim 0.6\mathrm{arcsec}$ ($\sim 40\mathrm{px}$). This approximately coincides with the expected location of the control radius, which delimits the region of the PSF affected by AO correction. This ring-like structure leaves a trace also in the residuals, shown in Fig. 3(b), and it is due to the method used in the adopted PSF-R algorithm to treat corrected spatial modes.

Fig. 3

Nighttime PSFs. (a) Observed PSF. (b) Reconstructed PSF. (c) Normalized residuals between observed and reconstructed PSF. Images size, color codes and scales are the same as in Fig. 2.

4.

Results and Discussion

In this section, we detail the performance of the adopted PSF-R algorithm. A qualitative analysis of Figs. 2 and 3 corroborates the good match between the observed and reconstructed PSF. In particular, the residuals shown on the right panels suggest that the core of the PSF is well reproduced. More quantitatively, a collection of the most relevant properties of the observed and reconstructed PSFs is given in Table 3. For each dataset and PSF, we report the measured SR and FWHM. The EE – defined as the amount of PSF energy that is contained within a fixed circular aperture, centered on the PSF center – in the core of the PSF, ${\mathrm{EE}}_{\mathrm{core}}$, is also listed. The radius of the aperture, in this case, has been fixed to 40.4 mas (2.7 px): the expected FWHM of the diffraction limited PSF. ${\mathrm{R}}_{50}$ is the aperture (in mas) that contains 50% of the total light. Finally, to quantify the match between the observed PSF and the reconstructed one, two parameters are listed in Table 3. The first is ${\chi }_{\mathrm{RED}}^2$, which is the reduced ${\chi }^2$ between observed and reconstructed PSFs. It has been obtained using GALFIT²⁴ on the observed PSF with the reconstructed PSF as input. The center position of a “psf” model has been set as the only free parameter for the fit and an input rms image has been created directly from the science data. The second parameter is the correlation level between the observed PSF and the reconstructed one. It is measured pixel-per-pixel over the whole PSF image ($8{\mathrm{arcsec}}^2$) by means of Pearson $R$ index. For comparison, the same parameter list is also provided in Table 4 for the diffraction limited case. We also include in the table the results associated to the best fit of the observed PSFs with an analytic profile, the sum of two Moffat functions. The choice is motivated by the fact that fitting a single Moffat function does not provide a good match on both the core and the wings of the observed PSF simultaneously (Ref. 25 and references therein).

Table 3

Properties of the observed and reconstructed PSFs for both datasets. See text for details.

Table 4

As in Table 3 but for the diffraction limited case and the best fit to the observed PSF with an analytic profile that is the sum of two Moffat functions.

Table 3 shows that the reconstructed PSFs match extremely well the observed ones, both for daytime and nighttime data. Indeed, the SR is recovered with an accuracy of 1.7% in both cases. To further discuss the other parameters, we compare the radial profile of the observed and reconstructed PSFs. They are reported in Fig. 4(a) for daytime and Fig. 4(b) for nighttime. The radial profile of the diffraction limited PSF is also provided for comparison. For completeness, the relative differences between observed and reconstructed PSF are also reported on the insets on top. In both cases, the radial profile of the reconstructed PSF nicely follows the observed one, particularly in the central regions. This is confirmed by the measured FWHM of the reconstructed PSFs that are accurate at the level of 1% and 4% for daytime and nighttime, respectively. This translates into differences of the order of 2  mas (0.13 px) in the worst case. Even the ${\mathrm{EE}}_{\mathrm{core}}$ in the core is fairly reproduced, with a precision better than 4% in both cases. Moving outward to around 50  mas (3.5 px), the light distribution of the reconstructed PSF keeps following the observed one with errors on $R_{50}$ of the order of 7% and 4% for daytime and nighttime respectively (0.25 px in the worst case). For larger apertures, small deviations (relative difference $<1$) are present in the radial profiles of the reconstructed PSFs with respect to the observed ones. In particular, we note the bump in the radial profile of the nighttime reconstructed PSF, associated to the ring-like structure visible in Fig. 3. This feature, however, does not seem to affect significantly the match between observed and reconstructed PSF, as demonstrated by the associated low ${\chi }_{\mathrm{RED}}^2$ and high correlation level values. Consistent values are found also in the daytime case.

Fig. 4

Radial profile (mean over all angles) of the observed PSF (black) and the reconstructed PSF (red). (a) Daytime and (b) nighttime. The associated standard deviations are indicated with shaded regions. For comparison, the radial profile of the diffraction limited PSF (D.L.) is also shown (blue dot-dashed line). The green profile on bottom shows the normalized deviation, as a function of distance, of the reconstructed PSF with respect to the observed one. The associated standard deviation is represented as a shaded area. We computed this quantity as the quadratic sum of the standard deviations associated to observed and reconstructed PSFs.

For what concerns the comparison with other PSFs, the diffraction limited PSF does not perform as well as the reconstructed one. All the values given in Table 4 indicate larger deviations than what has been obtained with the reconstructed PSFs. Furthermore, in the case of the two-Moffat fit, ${\mathrm{EE}}_{\mathrm{core}}$ and $R_{50}$ are well recovered, as expected, but the ${\chi }_{\mathrm{RED}}^2$ and correlation level values indicate that the reconstructed PSFs work equally well or even slightly better than the analytic profile.

4.1.

Impact of the PSF-R on an Ideal Science Case

The aim of this work is to obtain a reliable estimate of the differences between the reconstructed PSF and the observed one (as listed in Table 3) and to test the applicability of the reconstruction to an ideal science case. In this section, we aim to define a general method to assess how differences between observed and reconstructed PSF translate into errors on the scientific measurements. In fact, the same information should be used to define constraints on the reconstructed PSFs to meet the requirements for the precision on selected science measurements. For this reason, we limited as much as possible the effects of other factors related to specific observation. In particular, we envisaged a representative case where the reconstruction of the PSF is critical to reach the scientific goal and take advantage of the knowledge of the observed (true) PSF to assess quantitatively the potential differences that can be revealed as a comparison.

As an example, we consider the case of a distant (e.g., $z\sim 1.5$), bright and compact (e.g., effective radius comparable to the FWHM of the PSF) galaxy. We investigate what accuracy on the structural parameters is achieved when mode ling the 2D light profile of the galaxy with GALFIT using our reconstructed PSF model. For simplicity, we consider a galaxy whose light distribution can be described by a Sérsic profile. We did not simulate galactic substructures that, albeit more realistic, would complicate the interpretation of our results. Similarly, we chose to simulate a bright galaxy to assess the effects of the faint and outer region of the PSF, which would be hidden with a lower SNR but are fundamental because these are the most noticeable differences between the observed and the reconstructed PSF. The parameters of the simulated galaxy are: integrated magnitude $m_{\mathrm{Vega}}=17.5\mathrm{mag}$; Sérsic index $n=3.9$; effective radius $R_e=94\mathrm{mas}$ (6  px, corresponding to $\sim 2\mathrm{FWHM}$); axis ratio of 0.28. These values, with the exception of the integrated magnitude, are typical for a red nugget, an early type galaxy, at redshift $z\sim 1.5$.²⁶

To produce the mock observations, we make extensive use of SimCADO,^27,28 which we tuned accordingly to simulate SOUL + LUCI. We considered an exposure time of 54 000 s, subdivided into 180 frames (DIT) of 300  s each. To ensure a statistical significance of the results, we generated 100 mock observations of the same model galaxy, adding at each time random white noise.

The analysis is subdivided in two main parts: the simulation of the observations and the measurement. The simulated observations have been generated by convolving a pure Sérsic model with the observed PSF. A zoom-in of the central $128\times 128\mathrm{px}$ is presented in Fig. 5(a). The measurements have been performed with GALFIT giving different PSFs as input to the software: (i) the observed PSF, for the control measure; (ii) the reconstructed PSF, for the actual measure; (iii) the diffraction limited PSF, for comparison; and (iv) no PSF at all, again for comparison. The best-fit GALFIT model (a pure Sérsic profile) associated to the reconstructed PSF is shown in Fig. 5(b). The contours in the central region of the galaxy mark the same count levels as in the left panel: it can be noted the good agreement between the simulated galaxy and the model. This is emphasized by the normalized residuals, presented in Fig. 5(c). No evident residual structures can be detected.

Fig. 5

(a) Simulated galaxy. A $3.7{\mathrm{arcsec}}^2$ ($128\mathrm{px}\times 128\mathrm{px}$) window is shown. The nighttime observed PSF is used as input to generate the simulated observation. (b) The best-fit GALFIT Sérsic model. In this case, the reconstructed PSF is used as the input PSF for GALFIT. (c) Normalized residuals. The color scale is linear, the colorbar is shown on the bottom. The normalization has been done dividing the residuals by the simulated observation. It can be noted the excellent agreement between observation and model, despite the two being obtained using the observed and reconstructed PSF as input, respectively. The contours in the left and central panels refer to the same count levels.

Quantitatively, all the morphological parameters are in good agreement. We report in Table 5 the mean differences between recovered and input values along with the associated standard deviations. The table includes the results also for the control measure, the diffraction limited case, and the case where no PSF has been provided to GALFIT. The associated box plots are shown in Fig. 6, representing the distribution of the deviations from input values for the 100 simulations. The main parameters are shown, namely the integrated magnitude, the effective radius $R_{\mathrm{e}}$, the Sérsic index $n$, and the axis ratio $b/a$ (from top to bottom). The loss of accuracy is more noticeable for the case where no PSF at all is provided to GALFIT, confirming the need of the observed PSF in this particular case.

Table 5

Scientific evaluation results. The mean deviation from input values and associated standard deviation are reported for the main morphological parameters. Specifically, they are: the centre position (x and y), the integrated magnitude, effective radius, Sérsic index, axis ratio, and position angle of the galaxy. The input galaxy is obtained using the observed PSF, and each column refers to the different PSF used by GALFIT to recover the input parameters. In particular, the control measure has been obtained using the observed PSF also in the fit, whereas the one labeled “REC.” has been obtained using the reconstructed PSF. For completeness, we reported also the results obtained using the diffraction limited PSF and those without providing any PSF at all as input to GALFIT (“NONE”).

Parameter	OBS. CONTROL	REC.	D.L.	NONE
$X$ (mas)	$0.000\pm 0.164$	$-0.194\pm 0.179$	$-0.164\pm 0.164$	$-0.120\pm 0.284$
$Y$ (mas)	$-0.045\pm 0.105$	$-0.030\pm 0.105$	$-0.105\pm 0.105$	$0.000\pm 0.149$
Int. mag. (mag)	$-0.002\pm 0.019$	$-0.009\pm 0.016$	$0.020\pm 0.034$	$0.119\pm 0.017$
$R_{\mathrm{e}}$ (mas)	$-0.807\pm 3.005$	$2.497\pm 2.272$	$41.845\pm 11.168$	$69.5773\pm 6.055$
$n$	$-0.105\pm 0.324$	$-0.841\pm 0.188$	$0.421\pm 0.486$	$-1.640\pm 0.087$
b/a	$0.002\pm 0.005$	$0.042\pm 0.005$	$0.103\pm 0.006$	$0.217\pm 0.005$
PA (deg)	$0.095\pm 0.167$	$-0.621\pm 0.150$	$-0.896\pm 0.232$	$-0.737\pm 0.229$

Fig. 6

Box plots of the measured deviations from input values for integrated magnitude, effective radius $R_e$, Sérsic index $n$, and axis ratio $b/a$ from top to bottom. The control measure has been obtained using the observed PSF (the same used to produce the simulated observations) as input for GALFIT. On the contrary, the actual measure (REC.) has been obtained using the reconstructed PSF as input for GALFIT. For completeness, the box plots are shown also for the case where the diffraction limited PSF is given as input to GALFIT and when no PSF at all was provided as input to the software.

From Table 5 and Fig. 6, it can be noted that the recovered Sérsic index is somewhat biased toward smaller values. The best-fit model obtained providing the reconstructed PSF to GALFIT has an $n$ value $2\sigma$ lower than expected. As a general consideration, although an optimization of the analysis is possible (but beyond the aim of this study), our two-dimensional fit is performed successfully and the galaxy can be classified in all the cases as an early type object (having $2.5<n<4$).

To further assess whether the discrepancy between the measured and intrinsic Sérsic index is due to our PSF reconstruction, we also applied a mask on both observed and reconstructed PSFs selecting only the pixels inside a given cutting radius ($R_{\mathrm{CUT}}$). We then performed the same analysis described above, using the masked observed PSF to simulate the observations and giving the masked reconstructed one as input to GALFIT to recover the structural parameters of the galaxy. While gradually increasing $R_{\mathrm{CUT}}$, we record the variations in the measured morphological parameters.

We find that the only morphological parameter affected is indeed the Sérsic index. The behaviour of the measured $n$ as a function of $R_{\mathrm{CUT}}$ is shown in Fig. 7. If the central 29  px of the PSFs ($R_{\mathrm{CUT}}=3\mathrm{px}=44.85\mathrm{mas}$) are considered, the differences between the reconstructed and observed PSFs are negligible. Increasing $R_{\mathrm{CUT}}$, the recovered Sérsic index progressively shifts toward lower values, suggesting that the bias in the measure of this morphological parameter is due to the non-perfect reconstruction of the wings of the observed PSF. Specifically, for $R_{\mathrm{CUT}}>45\mathrm{mas}$, the difference between measured and input $n$ values steeply increases up to $-0.3$ at $R_{\mathrm{CUT}}\lesssim 90\mathrm{mas}$. The recovered $n$ value remains almost constant up to $R_{\mathrm{CUT}}\lesssim 0.4\mathrm{arcsec}$, from where a shallow increase in the deviation from input values is resumed.

Fig. 7

Mean deviation from the input value of the measured Sérsic index as a function of the cutting radius ($R_{\mathrm{CUT}}$) used to mask the external regions of both oberved and reconstructed PSFs. While the PSFs cores are very similar, the small differences in the wings of the PSFs translate into a deviation of the measured Sérsic index from the input value. After an initial settling to $-0.3$, the deviation resumes a shallow increase for $R_{\mathrm{CUT}}\gtrsim 375 \mathrm{mas}$.

5.

Summary and Conclusions

PSF reconstruction is a fundamental tool to analyze the data from AO instrumentation. It will support state-of-the-art scientific analysis of imaging and spectroscopic data. It is worth stressing here the fact that having a number of bright and isolated point sources in the observed field allows the PSF characterization directly from the science frames. This, however, can be tricky, in the case of AO observations where the PSF changes both with distance from the AO reference and in time.

The PSF-R technique allows to obtain the PSF of an astronomical observation simply from synchronous WFS telemetry data, without any access to the focal plane data. This is an unavoidable solution when there is no bright and isolated star in the instrumental FOV, suitable for deriving the observed PSF.

We have attempted here to carry out a PSF reconstruction of SOUL + LUCI LBT observations of bright, on-axis point sources to test the PSF-R code we are developing. To the best of our knowledge, this is the first published work to apply the PSF-R technique to real AO data acquired with a pyramid WFS. We have also carried out the PSF-R of an artificial source (an illuminated fiber), with optical turbulence disturbance generated and corrected using the Adaptive Secondary of LBT: these results have been already discussed in a preliminary stage by Ref. 16. The controlled atmospheric conditions allow to check that the tool is able to successfully recover the injected perturbations.

The PSF-R software we have used is able to reproduce the observed FWHM, EE, and SR with an accuracy of 2% to 4%, which is comparable to the level of accuracy reached by hybrid or calibrated techniques (see e.g., PRIME^29,30). Moreover, the reconstructed PSF provides a better fit to the observed data than the diffraction limited one. This level of agreement points to a successful approach of the adopted PSF-R algorithm, which matches a technology readiness level equal to 7, following the international ranking defined by the ISO 16290/2013 standard.

With the present work, we aim at setting the foundations for forthcoming analysis and refinements of the adopted PSF-R method. An effort is needed to improve the representation of the wings (at $>0.3\mathrm{arcsec}=20\mathrm{px}$) of the reconstructed PSF. It is worth mentioning that even when bright stars are present in the FOV, empirical PSF methods would not be able to provide an accurate description of the PSF wings. In the near future, we plan focussing on the treatment of off-axis (anisoplanatism) corrections and testing the PSF-R algorithm in the case of crowded stellar fields, e.g., in the case of a globular cluster.

This kind of studies paves the way for exploiting AO data reduction for ELT-size telescope observations of dense stellar fields or individual extra-galactic objects. The PSF-R software shown here will be an important starting point for developing a similar software for MICADO, the first-light camera of ELT, which will adopt a pyramid WFS similar to SOUL in its SCAO observing mode. In light of the future implementation of the MAORY module on MICADO, our group is also contributing in extending the algorithm to MCAO systems. It is worth noting that in the ELT the halo of a star will be a factor of $\sim 5$ larger than in LBT, and thus it will have a much lower surface brightness (a factor of $\sim 25$ fainter). On the other hand, the core of the PSF is a factor $\sim 5$ narrower in diameter. This means that, especially at longer wavelengths, even for bright stars the halo will be barely detected as a distinct feature, and that multiple halos will be blended into a faint distributed background. This indicates that the PSF-R approach will be even more important in the future for 30 to 40 m class telescopes. Only for stars that are bright enough to saturate, the outer halo becomes visible, thus the hybrid or calibrated PSF-R techniques will not be suitable in this case.

A first stage of scientific evaluation of the reconstructed PSFs has been carried out by simulating a long SOUL + LUCI LBT observation of a compact early type galaxy (red nugget) at $z\sim 1.5$. The reconstructed PSF allows us to derive the physical parameters (position, magnitude, half light radius) of the galaxy within 1 $\sigma$ level with GALFIT, whereas the Sérsic index is reproduced with an accuracy of $\sim 2 \sigma$. The axis ratio $b/a$ and position angle are reproduced within 5% and 0.6 deg, respectively. We propose that the possible systematic shift of the Sérsic index could be due to a non perfect representation of the wings of the reconstructed PSF with respect to the observed one. In summary, PSF reconstruction from telemetry data only is strongly required if one needs to obtain precise morphological measurements since typical high-$z$ galaxies are usually observed in blank extra-galactic fields. In this respect, the results of the presented analysis are extremely encouraging.

Lastly, we discuss, in the Appendix, an interesting relation linking EE, the FWHM, and SR of the observed PSF that can be useful, e.g., when one is interested in computing an instrument sensitivity.

We can conclude that the PSF-R approach presented here is an innovative method for pyramid WFS.

6.

Appendix A: A Practical Relation Linking SR, FWHM, and EE

In this appendix, we present a useful relation that links the observed EE (${\mathrm{EE}}_{\mathrm{OBS}}$) to the measured SR. Its application is of particular interest i.e. in the case that, for an AO systems, one is interested in the energy concentration in the core of the observed PSF knowing only the associated SR and FWHM with the final goal of computing an instrument sensitivity. The SR is defined as the ratio of the punctual peak values of the observed PSF and the diffraction limited one. A typical method for sensitivity computation is to consider the SR for the calculation of the EE within the PSF core:

To overcome this problem, a correction can be made in the above equation by means of a factor that accounts for the actual PSF shape. On a first approximation, this shape factor corresponds to the combined effect of the very-low order modes corrected by the AO. The shape of the core is dominated by the residuals of these few modes. In real cases, tip–tilt residuals represent the majority of the uncorrected low-order fraction for two main reasons: (i) because Kolmogorov power spectrum puts there 90% of the energy and (ii) because the largest non-turbulence wavefront error derives from telescope and instrument vibration (that may be as large as the atmospheric jitter and much faster). Moreover, the correlation time of the atmospheric tip–tilt—proportional to the ratio of wind velocity and telescope diameter—puts it in the 1 to 10 Hz,³¹ whereas the typical vibration power spectrum has peaks in the 1 to 100 Hz region. When the PSF peak is similar to the diffraction limited regime, corresponding to relatively high SR ($\mathrm{SR}>30\%$ to 50%), one may correct the equation above by adding the shape-factor:

If applied to the nighttime dataset, Eq. (6) gives ${\mathrm{EE}}_{\mathrm{core},\mathrm{OBS},\mathrm{uncorrected}}=0.36$, in contrast with the observed value of ${\mathrm{EE}}_{\mathrm{core},\mathrm{OBS}}=0.43$. The shape factor in this case results:

Equation (7) provides accurate results only under particular conditions. The major requirement is that the low-order aberration effect needs to be smaller than the EE aperture. Indeed, the higher is the broadening effect, the lower is EE inside a fixed aperture, therefore invalidating the normalization of the PSF. Specifically, the shape of the post-AO PSF inside the control radius matches as close as possible that of the diffraction limited PSF, but with intensity hampered by the amount of light outside the control radius associated to uncorrected modes.

We propose that a simple solution to avoid this issue is to take care in defining the aperture radius for ${\mathrm{EE}}_{\mathrm{core}}$ of the order of ${\mathrm{FWHM}}_{\mathrm{OBS}}$. In other words, we encourage the use of aperture referring to actual, post-AO light distribution that implicitly takes into account the broadening due to residual vibrations.

To quantify the limits of applicability of the relation in Eq. (7), we performed a simple simulation. We convolved the diffraction limited PSF, oversampled by a factor 5, with a 2D Gaussian functions of varying standard deviations (from 1 to $20 \mathrm{mas}$). In this case, as the observed PSF is obtained directly from the diffraction limited one, the lowering of SR values is only due to the residual vibration effects, which do not affect sensibly the wings of the PSF. Therefore, the effect of the aperture choice can also be estimated. We then measured the associated ${\mathrm{EE}}_{\mathrm{core},\mathrm{OBS}}$ and ${\mathrm{EE}}_{\mathrm{core},\mathrm{DL}}$ inside an aperture of radius $R_{\mathrm{EE},\mathrm{core}}{=\mathrm{FWHM}}_{\mathrm{OBS}}$. We also measured the SR directly from the resulting PSF. We then compared the measure ${\mathrm{EE}}_{\mathrm{core},\mathrm{OBS}}$ with those computed using Eqs. (6) and (7). The results are reported in Fig. 8. It can be observed that, after properly defining the aperture, Eq. (7) correctly estimates the observed ${\mathrm{EE}}_{\mathrm{core},\mathrm{OBS}}$ within 10% accuracy for a wide range of residual vibrations.

Fig. 8

Application of Eq. (6) (UNCORRECTED) and Eq. (7) (CORRECTED) to mock diffraction limited PSFs, each convolved with a 2D Gaussian function of varying standard deviations. The relative error in the resulting EE are shown as a function of the standard deviation of the convolving 2D Gaussian (${\sigma }_{\mathrm{VIBRATION}}$).