SPHERE adaptive optics performance for faint targets

Context: High contrast imaging is a powerful technique to detect and characterize planetary companions at large orbital separations from their parent stars. Aims: We aim at studying the limiting magnitude of the VLT/SPHERE Adaptive Optics system and the corresponding instrument performance for faint targets (G $\ge$ 11.0 mag). Methods: We computed coronagraphic H-band raw contrast at 300 [mas] and FWHM of the non-coronagraphic PSF, for a total of 111 different stars observed between 2016 and 2022 with IRDIS. For this, we processed a large number of individual frames that were obtained under different atmospheric conditions. We then compared the resulting raw contrast and the PSF shape as a function of the visible wave front sensor instant flux which scales with the G-band stellar magnitude. We repeated this analysis for the top 10\% and 30\% best turbulence conditions in Cerro Paranal. Results: We found a strong decrease in the coronagraphic achievable contrast for star fainter than G $\sim$ 12.5 mag, even under the best atmospheric conditions. In this regime, the AO correction is dominated by the read-out noise of the WFS detector. In particular we found roughly a factor ten decrease in the raw contrast ratio between stars with G $\sim$ 12.5 and G $\sim$ 14.0 mag. Similarly, we observe a sharp increase in the FWHM of the non-coronagraphic PSF beyond G $\sim$ 12.5 mag, and a corresponding decrease in the strehl ratio from $\sim$ 50\% to $\sim$ 20\% for the faintest stars. Although these trend are observed for the two turbulence categories, the decrease in the contrast ratio and PSF sharpness is more pronounced for poorer conditions.


Introduction
High-contrast direct imaging (DI) is a powerful technique for detecting and characterizing young and self-luminous giant planets, reflected light planets, highly polarized disks, Solar System objects, among other interesting objects. Thanks to its extreme Adaptive Optics (AO) system, called SAXO (Fusco et al. 2006), the Spectro-Polarimetric High-contrast Exoplanet REsearch instrument (SPHERE; Beuzit et al. 2019) is one of the most powerful instrument in the world for direct imaging science, and after more than 6 years of operations, it has led to the detection and characterization of a variety of (forming) planetary systems and protoplanet disks. Notable examples of these include β Pictoris b (Lagrange et al. 2018), the PDS 70 planet forming system , the multiplanet systems around HD 8799 (Zurlo et al. 2016) and YSES 1 (Bohn et al. 2020), the TW Hydrae disk (van Boekel et al. 2017), among others.
While different long-term SPHERE DI programs have mainly studied relatively bright stars (e.g. SHARDDS: Milli et al. 2017a; SHINE: Desidera et al. 2021;YSES: Bohn et al. 2021; BEAST: Janson et al. 2021), there is still a large population of stars located in young (age 10-20 Myr) stellar associations, which are ideal laboratories to study in-situ planet formation inside protoplanetary disks . However, these stars are embedded in gas and dust rich disks, which absorb a large fraction of the stellar visible light, making these object very red (see a discussion in Boccaletti et al. 2020). This effect is even stronger for M-type stars, which are intrinsically fainter in the optical when compared to the near infrared (nIR). The main dis-advantage of targeting these kinds of stars with SPHERE is that the instant flux recorded in the visible wave front sensor (WFS) is strongly reduced, and thus the AO correction is highly degraded. As more and more such faint stars (in visible light) are routinely being observed with SPHERE, it is mandatory to characterize the SAXO limiting magnitude. Motivated by this idea, we studied the SPHERE performance for targets fainter than G = 11.0 mag, where G corresponds to the Gaia Early Data Release 3 (Gaia EDR3) magnitude (Riello et al. 2021).
This paper is organized as follows: in section 2 we present the data used for the analysis we describe in details our methods. In section 3 we present the results of the raw contrast (from coronagraphic images) and the FWHM of the non-coronagraphic PSF, as a function of the WFS flux and the stellar G-band magnitude. Similarly, we computed the strehl ratio (SR) for a handful of selected stars and we compare it with the FWHM previously measured. In section 4 we compare the contrast from the raw images to that obtained in post-processing data. Finally, the discussion and conclusions are presented in section 5.

Data and method
We analyzed a total of 111 different stars observed with SPHERE between 2016 and 2022, all of them fainter than G = 11.0 mag, which are listed in Table A.1. Figure 1 shows the Gaia EDR3 color-magnitude diagram (CMD) with the position of all of the targets. Similarly, Figure 2 shows a 2MASS All Sky Catalogue (Cutri et al. 2003) nIR CMD for our stars. The median G-band and H-band magnitude of the sample is 12.4  We also included both imaging and polarimetric observations. In the case of the coronagraphic observations we included only individual frames with detector integration time (DIT) longer than 2 seconds. The list of header keywords used to filter the data are summarized in Table 1. We note that all of the raw data (noncoronagraphic FLUX and coronagraphic OBJECT data cubes) and their corresponding AO files (FITS binary tables containing atmospheric and SAXO telemetry data; see Milli et al. 2017b) were downloaded from the ESO archive. For the analysis, we developed an IDL-based code which computes the contrast for each individual raw frame, following the steps listed below. First, we fitted a 2-D Gaussian function to each individual non-coronagraphic FLUX frames. From the fit we estimated the peak flux ( f nc ) and the corresponding FWHM in the x-and y-pixels direction. The final f nc value was computed from the mean of all individual flux frames. We then used the mean f nc to predict the peak flux ( f c ) in each individual coronagraphic OBJECT frames, by scaling the DIT, neutral density filters and other optical elements with different transmissions (e.g. when a sequence is obtained with and without the apodizer). We then estimated the noise level (σ 300 ) at 24 pixel from the coronagraph center (corresponding to ∼ 300 mas), by computing the standard deviation in a 5 pixels wide annulus. The resulting 5-σ raw contrast at 300 mas (RC 300 ) was obtained from: We note that for broad band observations (B_H filter) and differential polarimetry (DP_0_BB_H filter) we used only the left side of the Infra-Red Dual-beam Imaging and Spectroscopy (IRDIS; Dohlen et al. 2008) detector, while for the dual-imaging mode (DB_H23 filter) we analyzed independently the left (λ c ∼ 1.59 µm) and right (λ c ∼ 1.67 µm) side of the detector, and we computed then mean value from the two channels. Finally, to study the effect of the atmospheric conditions on the contrast and the PSF, we retrieved the Differential Image Motion Monitor (DIMM) seeing and coherence time (τ 0 ) at the beginning of each exposure, directly from the image header.

Instant flux on the visible WFS
As part of the SAXO system, SPHERE is equipped with a 40 x 40 Shack-Hartmann WFS, (with a total of 1240 effective subapertures imaging the SPHERE entrance pupil). The light reaching the WFS is recorded by a 240 x 240 pixels electron multiply-ing charge-coupled device (EMCCD) detector, sensitive to visible light (∼ 0.45-0.95 µm). The main advantage of this detector is that it can operate at very high frequency (up to 1380 Hz), low read-out noise (< 0.2 e − /pixel/frame) and high multiplication gain (up to M = 1000; Sauvage et al. 2016b). However, as mentioned before, for our analysis we included only data taken with the EMCCD operating at 300Hz and M = 1000, which is the mode automatically selected by the real time computing system (SPARTA; Fedrigo et al. 2006) for faint targets (R > 10.5 mag). The instant flux in the WFS EMCCD detector ( f WFS ) is recorded by SPARTA, and it is saved in the AO files. These files (among other telemetry data) contain the total flux per pupil, i.e., integrated in all 1240 subapertures. This flux (tagged FLUX_AVG in the binary tables) is thus expressed in units of [ADU/frame/pupil]. To convert FLUX_AVG into [e −1 /subaperture/frame], we used: where g = 16.5 [e − /ADU], n s = 1240 [subapertures/pupil] and M = 1000 (unit-less).

Raw contrast and FWHM as a function of the WFS flux
After computing the raw contrast, we cleaned-up the sample by removing those frames leading to very deviant values. Typical examples of these include cases with the stars out of the coronagraph, AO loop open, contamination from close binary companions, extreme cases of low-wind effect (Sauvage et al. 2016a;Milli et al. 2018), among others. We then compared the raw contrast as a function of f WFS . Figure 3 shows the RC 300 versus the f WFS binned in 1 [e − /subaperture/frame] flux bins. We note that only frames with header values listed in Table 1 were included (as explained in section 2). The black dots correspond to atmospheric turbulence category TCAT10 (seeing < 0.6 arcsec; τ 0 > 5.2 ms) and the red dots to TCAT30 (seeing < 0.8 arcsec; τ 0 > 4.1 ms). The error bars correspond to 1σ equaltailed confidence interval for each bin. As can be seen, there is a sharp decrease in the achievable contrast for f WFS below ∼ 3 [e − /subaperture/frame]. In this regime the AO correction is degraded due to the lack of photons collected by the WFS. We also note that the RC 300 achieved under TCAT10 conditions is slightly better than that under TCAT30 . Similarly, Figure 4 shows the mean FWHM (from the x-and y-direction) measured on noncoronagraphic frames as a function of f WFS . Again this quantity is computed on the left side of of the IRDIS detectors in the dualband mode (see Table 1). It can be seen that for f WFS level above ∼ 3 [e − /subaperture/frame] the point-spread-function (PSF) is nearly flat with a FWHM of ∼ 53 [mas], for both TCAT10 and TCAT30 . On the other hand, for lower f WFS values, the AO correction degrades and thus the FWHM increases with decreasing flux level, reaching a mean value larger than ∼ 100 [mas], when reaching the sub-electron per frame regime, for both TCAT10 and TCAT30 .

Correlation between stellar magnitude and WFS flux
To understand the SPHERE performance not only as a function of f WFS but also as a function of the target magnitude, we investigated the correlation between f WFS and the stellar G-band magnitude. We choose this particular band because the sensitivity curve of the WFS detector is very close to the G-band trans-  mission curve and also since nearly all of the targets observed by SPHERE show a G-band entry in the EDR3 catalogue. Figure 5 shows the normalized transmission of the Gaia EDR3 bands 1 , the Johnson R broad band filter Bessell (1990) and the WFS response, including all of the upstream optical elements and with the WFS spectral filter in the OPEN position (which is the case for faint stars). It can be seen that the G-band transmission function is the closest one to the WFS response function, although the latter one is slightly more sensitive to longer wavelength. For this reason, to predict the expected f WFS value during an observing sequence of a faint target, the G-band should be used, instead of the R-band (as currently implemented in SPARTA). In fact, even though in general the value of these two magnitudes are very close, there are cases where the difference is as large as one magnitude (corresponding to a factor ∼ 2.5 in flux). For this analysis, we first corrected the flux in the WFS for the target airmass, using an extinction coefficient of 0.12 (Patat et al. 2011). We then compared this value with the stellar G-band magnitude. We note that to minimize the stellar color effect, we only included stars with 1 < G B -G R < 3 (see Figure 1). The results are shown in Figure 6. The error bars in the WFS flux correspond to the standard deviation from all individual measurements. The solid line corresponds to the WFS magnitude scale (m WFS ), given by: We note that we removed outliers beyond 2.5 σ from the fit, leading to a scatter of 0.18 mag. This value is mainly explained by observations performed under different observing conditions (sky transparency, seeing, etc) and more importantly due to the fact that we included observations using different spatial filters (Fusco et al. 2016), as listed in Table 1, meaning that a different amount of the light is blocked for different observations. Ideally, this analysis should be repeated with data taken under photometric conditions and using the WFS spatial filter in the OPEN position, as opposed to the majority of the observations analyzed here that were done using either the SMALL or MEDIUM spatial filter.
Using this correlation we can directly compare the stellar G magnitude and the expected raw contrast, for different atmospheric conditions. We therefore replicated Figure 3 and Figure  4, but this time we replaced f WFS by the stellar G magnitude. The results are shown in Figures 7 and 8. As can be seen, there is roughly one order of magnitude decrease in the contrast ratio beyond G ∼ 12.5

Strehl ratio
In addition to the FWHM values, we investigated the dependence of the strehl ratio (SR) with the f WFS and the stellar magnitude. For this, we first computed a theoretical diffraction limited PSF for the VLT as seen by SPHERE, (i.e, including the effect of the Lyot stop and the apodizer in the optical path) at 1.63 µm . The PSF was computed over an area of 29x29 pixels (including light up to the third airy ring) with the IRDIS plate scale of 12.27  The solid line corresponds to the best linear fit between these two quantities. The position of AX cha is also shown (red star). Fig. 7. 5-σ raw contrast at 300 mas as a function of the stellar G-band magnitude. The black and red lines correspond to atmospheric turbulence category TCAT10 and TCAT30, respectively.
[mas/pix], and was then normalized by the total integrated flux in all 29x29 pixels area. A 3D view of the resulting normalized PSF is shown in Figure 9. We then fitted a two-dimensional Moffat profile (Moffat 1969) of the form: where U corresponds to: The resulting coefficients are listed in Table 2. We then repeated the Moffat fit, but this time to a total of 332 observed FLUX frames. We note that in this case we included also brighter stars (leading to higher SR) that were not part of the sample presented in section 2. The resulting SR was simply estimated by taking the ratio between the peak of the Moffat fit (corresponding to I 0 in eq. 4) of the observed FLUX frames and that from the theoretical PSF. An example of an observed PSF leading to SR = 0.78 and its corresponding Moffat fit is shown in Figure 10.    In addition, we measured the FWHM to the FLUX frames (as explained in section 3.1), and we compared them with the measured SR. The results are presented in Figure 11. We finally Table 3. Resulting coefficients for the SR-FWHM correlation presented in equation 6.

Coefficient
Value a -0.099 b 5.296 c -33.236 α 0.716 obtained a correlation between the SR and FWHM by fitting a function of the form: The resulting coefficients are listed in Table 3. The RMS of the post-fit residuals is 0.023. The best-fit is also over plotted in Figure 11 Using this correlation, we then converted the FWHM values into SR. We therefore replicated Figures 4 and 8

Final contrast
To investigate what is the achievable final contrast as a function of the stellar magnitude, we compared the raw contrast with that measured after post-processing. We included in this analysis 19 different stars, leading to RC 300 values covering roughly two orders of magnitude (between ∼ 2×10 −2 and 2×10 −4 ). For this we retrieved the post-processing contrast curves (corrected for the throughput of the algorithm) produced by the SPHERE Data Center (Delorme et al. 2017), for the classical Angular Differential Imaging (cADI; Marois et al. 2006) and the Template Local Optimized Combination of Images (TLOCI; Marois et al. 2010;Galicher et al. 2018) techniques. We note that we only included relatively long observing sequences, leading to a total field of view (FoV) rotation > 15 degrees, to avoid strong self-subtraction and thus low throughput. These results are presented in Figure 14. The black dots and blue stars correspond to cADI and TLOCI reduction techniques. The dashed black and blue lines correspond to the median conversion factor, corresponding to 2.7 and 12.4 for the cADI and TLOCI techniques respectively. We can then use this result to estimate the final contrast on post-processing data by scaling RC 300 using these two values (e.g. to convert Figures 3 and 7 into final contrast). Finally, we repeated this conversion, but this time from RC 300 to final contrast at 500 [mas]. We obtained median conversion factors of 11.0 and 43.5 for the cADI and TLOCI techniques, respectively.

Summary and discussion
In this paper we have analyzed a sample of 111 relatively faint stars (G ≥ 11.0 mag) observed with IRDIS/SPHERE between 2016 and 2021, with the aim of investigating how the AO correction (and thus the final performance in terms of achievable contrast) degrades when approaching the limiting magnitude regime. For this we measured the H-band RC 300 from coronagraphic frames and we compared this quantity to the instant flux received by the visible WFS. We found that the RC 300 is at the ∼ 10 −3 level for f WFS 3 [e −1 /subaperture/frame].
A&A proofs: manuscript no. sphere_performance  For f WFS below this value, the AO correction strongly degrades, decreasing the RC 300 by a factor ∼ 10 when reaching the 1 [e −1 /subaperture/frame] level. In terms of stellar magnitude, this turning point is observed at G ∼ 12.5 [mag]. At G ∼ 14 [mag] the measured RC 300 is at the 10 −2 level. We observed this behaviour in the data taken under both TCAT10 and TCAT30 atmospheric conditions, with the former one performing slightly better.
Similarly, we computed the FWHM and SR from the FLUX frames. We observed a very similar trend above and below the ∼ 3 [e −1 /subaperture/frame] limit. In particular, we measured a factor ∼ 2 and 3 increase in the FWHM from this limit and the ∼ 1 [e −1 /subaperture/frame] f WFS level, for TCAT10 and TCAT30 , respectively. Similarly, the SR drops from above 0.5 to less than 0.2, for both TCAT10 and TCAT30 .
Finally, we computed the final contrast at 300 [mas] and 500 Fig. 12. Strehl ratio as a function of the visible WFS flux. The black and red lines correspond to atmospheric turbulence category TCAT10 and TCAT30 , respectively.
[mas], in a selected sample of stars, using the post-reduction contrast curves. We obtained a median conversion factors at 300 [mas] of ∼ 3 and 12, for the cADI and TLOCI techniques, respectively. Similarly, we found scaling factors to 500 [mas] of ∼ 11 and 44 for cADI and TLOCI, respectively. Using these conversion factors, we can then estimate the final contrast achievable at 300 [mas] and 500 [mas], by scaling the results presented in this work.
As an example, let us consider a star with the exact same luminosity and colors of the well known planet host star PDS70 ). Let us also consider that such star is located at ∼ 213 [pc] from the Earth, so that it has: G = 13.0 [mag] and H = 11.2 [mag]. Following Figure 7, we can expect a RC 300 of ∼ 4×10 −3 under good atmospheric conditions Fig. 13. Strehl ratio as a function of the stellar G-band magnitude. The solid black and red dots (lines) correspond to atmospheric turbulence category TCAT10 and TCAT30 , respectively.

Fig. 14.
Comparison between the contrast measured at 300 [mas] from individual raw frames and post-reduction full observing sequence. The black dots and blue stars correspond to cADI and TLOCI reduction techniques, respectively.
(TCAT30 or better). If we perform a relatively long (FoV rotation 15 deg) ADI sequence, we can then estimate a 5σ final contrast at 500 [mas] using the TLOCI technique of ∼ 10 −4 . If we compare these number with the AMES-COND models (Allard et al. 2001), this limit corresponds to a 10 [Myr] planet with a ∼ 4.8 M J mass and located at a physical separation of ∼ 107 [AU] from the host star.
Improving the contrast on faint targets and reaching fainter or redder stars represent one of the primary scientific drivers for the instrument upgrade called SPHERE+ (Boccaletti et al. 2020). Indeed ALMA surveys such as DSHARP (Andrews et al. 2018) revealed dozens of star-forming disks with structures such as gaps, rings or spirals, suggesting planet formation is on-going. PDS70 represents one of the brightest examples of this population, the bulk of it being currently out of reach for SPHERE (see Fig. 12 of Boccaletti et al. 2020). Coupling the current AO correction stage to a more sensitive second-stage AO in cascade represents the main technical improvement of SPHERE+ (Cerpa-Urra et al. 2022) and the characterisation of the current performance on faint targets presented here is crucial to make informed technical choices in the design of the second stage AO for SPHERE+.