© The Authors 2023

1 Introduction

Since the discovery of the first giant planets around solar-type stars almost 30 yr ago, thousands of planets have since been discovered, with masses down to a few Earth masses. However, Solar System analogues¹ have not been detected thus far, and so, we still do not know if our own planetary system is indeed unique. Detections of Earth twins have still not been made possible, as this requires major on-going research to correct for stellar activity with radial velocity (RV) techniques at the appropriate level (Meunier & Lagrange 2019). Also, detecting Jupiter twins orbiting at 5 au is very challenging (and not possible at larger separations) with radial velocity because decade(s) of careful monitoring are needed; additionally, long-term stellar activity is also responsible for a long-term noise. As a consequence, the orbital parameters of the (still rare) RV planets announced beyond 5 au are poorly characterized (Wittenmyer et al. 2016; Fernandes et al. 2019; Fulton et al. 2021). These limitations stand in the way of making precise comparisons between the radial distribution of giant planets beyond their forming regions and predictions from population synthesis models meant to constrain formation scenarios (Lagrange et al. 2023).

Giant planets (GPs) have played a significant role in the building of the Solar System (see e.g., Levison & Agnor 2003; Raymond et al. 2014; Morbidelli et al. 2012) and of exoplanetary systems (see e.g., Quintana & Barclay 2016; Childs et al. 2019). Various mechanisms may be involved, among which dynamical interactions with lighter bodies (e.g., telluric planets, planetesimals) once the proto-planetary disk has dissipated and the dynamics is no longer controlled by gas. In particular, GPs may play a role in the development of life on Earth analogues (Horner & Jones 2010) and could even have driven the delivery of water on Earth (Morbidelli et al. 2012). From an observational point of view, remote GPs may also impact the detectability of lighter and closer-in planets with the RV technique and with astrometry because of the more complex RV or astrometric signals in case of multiple systems. Hence, knowing their giant planet population is key to modeling individual systems.

Detecting giant planets is therefore crucial to understanding planetary system formation and evolution. While RV or transit techniques are best suited to detect GPs orbiting typically within 5 au, they are not well adapted for detecting and characterizing more remote ones. Absolute astrometry is well adapted for giants in the 5–10 au (Perryman et al. 2014), even though such long period planets may be difficult to fully characterize (Ranalli et al. 2018), especially in the case of multiple systems. Micro-lensing will also be very useful to constrain the giant planet demographics in the 5–10 au range (Beaulieu & Bachelet 2021). High-contrast direct imaging (DI) is probably the most promising technique for detecting and characterizing analogues of our Solar System's giants planets in the future. Yet current high-contrast instruments like SPHERE (Beuzit et al. 2019) or GPI (Macintosh et al. 2014) are sensitive to massive young giants orbiting typically beyond 10 au. As an example, the SPHERE SHINE GTO survey had a 20% chance of detecting a 2 M_Jup at 20 au, and a 10% chance of detecting a 4 M_Jup planet at 5 au (Vigan et al. 2021). These poorer performances are due to 1) limitations at the instrumental and data processing levels and 2) a geometrical effect, namely: unless they are on pole-on orbits, short-period planets may be missed in a single observation because of a small projected separation at the time of the observation. For instance, due to geometrical effects, a single observation only allows us to explore half of the 5 au region around a star located 30 pc away when the planet orbit is seen edge-on. Fortunately, this geometrical effect can be easily overcome: the region explored is increased by more than 70% by observing the star twice, each time a few years apart (Lannier et al. 2017).

In this paper, we apply PACO (patch covariance), a promising detection algorithm (Flasseur et al. 2018, 2020a,b) on a small sample of stars representative of SPHERE targets, which are members of young close associations observed as part of the SPHERE/SHINE survey (Desidera et al. 2021) and observed under a wide range of atmospheric conditions. This analysis constitutes a test-bed for a forthcoming massive reduction of the SPHERE archive. Our aim is to find the best analysis strategy and to estimate the detection limits achievable on these stars. Moreover, we define our sample so as to address the astrophysical question of how far we are from detecting young giant planet siblings of our Solar System.

Our paper is organized as follows: the sample and the data are described in Sect. 2. The data reduction and analysis are described in Sects. 3 and 4, respectively. Section 5 describes the achievable performance and Sect. 6 presents the astrophysical results. A brief summary and a presentation of future work are provided in Sect. 7.

Fig. 1 Properties of the sample stars in terms of distance, age, spectral type, and magnitude.

2 Star sample

Our sample is a collection of all (24) young (≤150 Myr), close by (<60 pc), solar-type stars (FGK) observed during the SPHERE/SHINE survey early science release (Desidera et al. 2021) and previously analyzed with conventional postprocessing algorithms by (Langlois et al. 2021, hereafter, F150). The thresholds in age and distance were chosen to ensure the best detection limits (5 M_Jup, down to possibly 1 M_Jup) possibly down to 5–10 au from the stars, namely, at the locations of the Solar System giants². Figure 1 and Table A.1 show the properties of the stars in our sample directly extracted from Desidera et al. (2021). It can be noted that the histogram of the ages of the considered stars is bimodal, with one subset aged between 20 and 60 Myr and the other aged about 150 Myr. This bimodal distribution is caused by some stars belonging to young co-moving groups such as AB Doradus (ABDO, 150 Myr), Tucana-Horologium (TUC, 45 Myr), Carina (CAR, 45 Myr), or β Pictoris (BPIC, 24 Myr).

All our targets were observed in angular (and spectral) differential imaging (A(S)DI, Marois et al. 2006) using the telescope in pupil tracking mode. The standard observing mode of the SHINE survey, namely, the IRDIFS mode was used, with IRDIS dual band images in H2 and H3 (Dohlen et al. 2008) and IFS (Claudi et al. 2008) data covering the YJ bands³. Tables A.2 and A.3 provide descriptions of the observations and of the associated atmospheric conditions during these observations. The stars were observed around the time of the meridian crossing, to ensure the largest amplitude of parallactic angle variations (of at least 30 degrees). The atmospheric conditions were heterogeneous, with a seeing ranging from 0.5" to 1.0" most of the time and the coherence time ranging between 1 and 10 ms.

3 Data reduction and frame centering

The reduction pipeline from raw to centered datasets is similar to the one built for the SHINE survey, described in Delorme et al. (2017) and Langlois et al. (2021). To improve the centering of the IRDIS and IFS frames, a custom built routine that uses the waffle center calibration has been developed (Dallant et al. 2023). As part of the observing sequence, before and after the corona-graphic sequence, two coronagraphic images are recorded with a waffle pattern applied to the deformable mirror to create four replicas of the point spread function (PSF) at a separation of about 14 λ/D from the central star. These replicas, called "satellite spots", are used to determine precisely the position of the central star behind the coronagraphic mask before and after the long coronagraphic sequence.

To determine the accurate positions of the satellite spots, small circular regions are extracted around their theoretical positions and a bi-dimensional, non-istropic Gaussian fit is performed using a trust region reflective algorithm (Branch et al. 1999), particularly suited for large sparse problems with bound constraints. Estimates of the central star positions are then computed via the centroid of the resulting fitted satellite spots and their associated uncertainties. The frames are re-centered using the mean value of the two estimated centroids. This new routine is marginally more precise than the one currently implemented in the SPHERE data center, but its main advantage is a much faster computational time: assembling the 4D datacube takes only a few minutes, that is: more than one order of magnitude faster than the current pipeline, without any loss in precision.

We note that when precise astrometric measurements of known companions are needed, the satellite spots are generated during the whole coronagraphic sequence. In such cases, all frames are recentered individually using their own satellite spots.

4 Analysis pipeline

The analysis pipeline is based on the PACO A(S)DI pipeline described in Flasseur et al. (2020a,b). A few important upgrades were made, as follows: (1) an improvement in the PACO robustness, namely, the capability to run PACO on diverse and heterogeneous datasets, while consistently providing reliable results (see Sect. 4.1); (2) an optimization of spectral priors for PACO ASDI (see Sect. 4.2); (3) an automated and improved computation of astrometric and photometric error bars for each characterized source (see Sect. 4.3); and (4) an automated classification of the status of any identified candidate companion in case of multi epoch observations (see Sect. 4.4).

Finally, in view of the forthcoming massive re-reduction of all SPHERE data, a tool was developed to automatically identify any potential companion, and gather associated astrophysical information (astrometry, photometry, spectra, etc.) required for further analysis. These upgrades are described in the following subsections. Both the centering routines and the analysis pipeline are hosted on the COBREX data center, a modified and improved server based on the SPHERE data center.

4.1 Improvements of PACO robustness

The principle of the PACO algorithm is described in Flasseur et al. (2018). No fundamental modifications were made concerning the core and the technical elements of the methods. The main updates are: (1) a refinement of the PSF fitting routine, now implementing a robust strategy based on iteratively re-weighted least-squares (Huber 2011). This routine improves the robustness of the fit in the case of very low signal-to-noise ratio (S/N) in the measured off-axis PSF (e.g., in absorption bands); and (2) the management of time-variable missing data, with a time-variable mask, to account for possible evolution (during the sequence of acquisition) of the field-of-view with exploitable data.

In addition, an engineering effort has been made to validate, via numerical experiments, the faithfulness of the astrophys-ical quantities produced by the algorithm, particularly those concerning the astrometry and photometry, as well as their associated error bars. Due to residual noise, the flux estimate in the absence of a source is not completely zero on average. The average level of this effect was estimated for both ADI and ASDI. Thus, flux estimates for which the detection confidence is less than 1σ in ASDI and 2.5σ in ADI were not considered or used in this analysis or in the massive reduction. Along the same line, code upgrades (accelerations, automations, and case-specific handlings) have been implemented to allow for massive reductions performed on a computer server.

4.2 PACO ASDI spectral priors to increase the sensitivity

The ASDI mode of PACO offers the possibility to combine multi-wavelength datasets into a detection map, using specific weights to maximize the detection efficiency. These weights, {w_ℓ}_{ℓ=1:L ∈} [0; 1], are referred to as the spectral priors They are represented by vectors with as many components as wavelengths: L = 2 for IRDIS data and L = 39 for IFS. Since all the photometric measurements within PACO are expressed with respect to the target star, the priors should also be expressed as the expected companion contrast relative to the host star (shape-wise). They are then normalized between 0 and 1 (dynamic-wise). When simultaneously using multiple priors, the PACO algorithm computes the S/N of the detected source for each prior. As an illustration, we show in Fig. 2 the S/N measured on point sources considering various priors for IRDIS data in the case of an injected fake planet. In Fig. 3, we show the same kind of plot for an IFS dataset and a real point-like source: HD 206893b (Milli et al. 2017). A more classical ASDI spectral combination approach, as the one implemented in the PCA and TLOCI versions of the SPHERE data center, is somewhat similar⁴ to the "flat prior" combination (i.e., assuming that the sought-for exo-planets have the same spectral energy distribution (SED) as their host stars), which is highlighted in Figs. 2 and 3.

As expected, the S/N is higher when the spectral prior is similar to the planet spectrum. Thus, we must define sets of spectral priors representative of the variety of the spectra of the potential exoplanets to optimize the detection capabilities. Besides, while increasing the number of priors improves the sensitivity to different types of objects⁵, it also significantly increases the computational time and (moderately) increases the number of non-redundant false positives at a given detection threshold (see Sect. 4.2.2). A trade-off must then be found.

Fig. 2 Evolution of the S/N (blue) for various priors for a point source observed with IRDIS, and whose spectrum corresponds to the red mark. The x-axis gives various values of the priors. The S/N corresponding to the [0.5, 1] prior (significantly different from the normalized contrast of the injected source) is below the 4σ threshold used for this analysis, while the S/N reached for spectral prior [1, 0.1] (similar to the normalized contrast of the injected source) leads to a clear detection with a significance above the 5σ detection threshold (dashed black line).

Fig. 3 Evolution of the S/N (blue) for various priors for HD 206893b (2017-07-13, IFS, YJ band). The x-axis gives the index of the priors. The S/N corresponding to the 10 last priors are below the 5σ detection threshold. Priors are ordered by decreasing S/N for clarity purposes. The S/N reached for spectral prior close to the spectrum (in contrast) leads to a detection with a significance above the 5σ detection threshold.

Fig. 4 Integrated H2/H3 flux ratio computed for a planet with a solar metallicity, C/O=0.50, and with different values of log(g) and Teff. ExoREM spectra were used.

4.2.1 Selection of spectral priors using fake planets injections

To select the spectral priors, we used a set of four targeted stars that were used in a recent internal blindtest conducted by the SHINE consortium that was aimed at comparing the performance of various detection and characterization algorithms. The targets properties are provided in Table A.4, and the observing and atmospheric conditions are given in Table A.5. Several hundreds fake planets with various properties were randomly injected between 0.12" and 5.5" in the case of IRDIS data -while still avoiding any blending among injected sources. For IFS data, about 100 fake point sources (FPSs) were injected. The injected FPS spectra were taken from the BT-Settl grid (Allard 2014) and the priors were built using the ExoREM (Charnay et al. 2019) spectra. We purposely used different models to inject the fake planets and to build the priors, to avoid biases that could occur when using the same library for both the injection process and the prior definition.

4.2.2 Priors selection for IRDIS

To find the trade-off between sensitivity and computation time, we considered all possible available spectra, and for each raw spectrum, we computed the ratio (or ), where H2 and H3 (respectively K1 and K2) represent the integrated fluxes in these spectral bands. Figures 4 and 5 show examples of such ratios for a source with various effective temperatures, a solar metallicity and a C/O ratio of 0.5 for both filter combination. Two regimes can be identified: (i) a "cold" one, corresponding to T_eff roughly below 1000 K, where the (resp ) ratio varies from several tens even hundreds down to 1 mostly because of strong CH₄ absorption in H3 and K2 filters, and (ii) a "hot" one, corresponding to T_eff above 1000 K, where this ratio is almost constant between 1 and 0.6. A similar behavior is observed for all metallicities and C/O ratio.

The goal is to find the optimal number of priors to explore these features. To do so, we proceeded as follows:

We first used all the injected FPSs and measured in each case the S/N considering different priors. These priors correspond to the 15 priors showed in the x-axis of Fig. 2. In this example, the maximum S/N is 6.7. We then repeated the process using instead eight spectral priors uniformly spread over the considered parameter space (by steps of 0.2), five (by steps of 0.3), and four (by steps of 0.4). We then considered the evolution of the maximum S/N as a function of the number of priors (for each injected FPS). In the case of relatively faint FPSs (sources that cannot be identified on a single frame), using five priors (or more) does not significantly degrade the S/N (less than 1%), while using four (or less) does.

Second, we then studied the impact of the number of priors (between 1 and 8) on the rate of false positives. To do so, we computed the average number of false positives identified by PACO in the IRDIS FoV (limited to 5.5" to avoid edge effects) on a total of 20 datasets and compared it to the theoretical value. We can compute the theoretical expected number of false positives as follows: because the pixel distribution on the S/N maps is Gaussian, the number N_fp of false positive per map according to the number n_pixel of pixels processed and the probability of false alarm P_FA(τ) at a given detection threshold τ: (1)

As an illustration, the numberofpixels to process in each IRDIS dataset is approximately equal to one million and, with a 5σ detection confidence (i.e., τ = 5, P_FA(5) ≃ 2.87 × 10⁻⁷), we have:

with false alarms expected in each detection map.

Figure 6 shows the results for τ = 5 with an increasing number of spectral priors used. The false positive rate when considering a single S/N map corresponding to a given prior (no matter which prior is used) is in good agreement with what is expected from a Gaussian noise distribution (see Eq. (1)). When working with several priors, the number of detections will only increase if a new independent source (i.e., detected for the first time in the current prior) is detected. Redundant detections will only be accounted for one source. Because false positives are often redundant, the empirical cumulative false positive rate is lower than the theoretical cumulative one (i.e., if all false positives were independent).

This confirms the Gaussian nature of the S/N map produced by PACO as well as the associated statistical guarantees (i.e., control of the probability of false alarms and of detections). Based on this study, we chose to include five spectral priors in our library Ω_IRDIS ∈ ℝ^2×5, achieving the desired trade-off between maximizing our detection performance and lowering the number of false positives:

Fig. 5 Integrated K1/K2 flux ratio computed for a planet with a solar metallicity, C/O = 0.50, and with different values of log(g) and Teff. ExoREM spectra were used.

Fig. 6 Average number of false positives for a detection threshold τ = 5 with one to eight spectral priors with IRDIS. Twenty observations were considered for these experiments. The red dashed line shows the theoretical false positive rate on a single S/N map. The mean experienced false positive rate for each S/N considered independently map is represented by the yellow dashed rectangles and the mean experienced cumulative false positive rate by the blue rectangles.

4.2.3 Priors selection for IFS

In this section, we describe how we build a library Ω_IFS ∈ ℝ^N×L of N spectral priors for processing the IFS data with PACO ASDI.

Following the notation introduced in Sect. 4.2, each element w ∈ ℝ^L of Ω_IFS is a vector with as many components as wavelengths (i.e., L = 39 for the IFS). In practice, we build a different set Ω_IFS for each stellar spectral type since each element w should be expressed in contrast units, so that it depends on the spectral type of the star explicitly. For the purposes of illustration, the method is described for any given stellar spectral type without loss of generality. The set Ω_IFS is built from a set of sub-stellar spectra provided by the ExoREM models. The spectral resolution L_ER of each element is much larger than L, with a typical value of L_ER = 500 for this study.

For a given stellar spectrum, , and a sub-stellar spectrum, , at the same spectral resolution, L_ER, we obtained an (intermediate) spectral prior by dividing the two elements component-wise, namely: , . Each intermediate spectral prior w' is then normalized by the maximum value over its components. The last operations aim to reshape w^' from L_ER = 500 to the spectral resolution, L, of the measurements by first applying a convolution with a Gaussian kernel of standard-deviation in the order of magnitude of L (typically between 30 and 50) and then re-sampling the results at the targeted spectral resolution L to get a final spectral prior w.

We aimed to include in Ω_IFS the minimum number N of spectral priors needed to represent the diversity of the observations. For that purpose, we build Ω_IFS from a large set Ω_ΕR (i.e., N_ER ≫ N), and we progressively add non-redundant atoms w in Ω_IFS. We first started with a single spectral prior w^ref of reference in the set Ω_IFS. This spectral prior is built from a sub-stellar model of reference obtained from the ExoREM simulator with the following parameters: T_eff = 1000 K, Fe/H = 1.0, C/O = 0.50, and log(g) = 4.0. By looping over the (fixed) elements of Ω_ER, we compute the Euclidean distance between each resulting candidate spectral prior and the spectral priors already present in Ω_IFS. To the set Ω_IFS, we then added the candidate spectral prior that maximized the distance, namely, the spectral prior that differs the most from the already selected ones. For practical reasons, we preset the number, N, of targeted elements in Ω_IFS to perform the above selection procedure. We repeated this procedure for various values of N and selected the number that gives a satisfying trade-off between precision and recall, while simultaneously leading to a manageable computation time at data reduction time. As an illustration, we find that reducing the number N of spectral priors in Ω_IFS from 31 to 20 decreases the detection capabilities only marginally (by less than 1% in terms of S/N loss). In addition, we find that decreasing the number N of spectral priors in the same proportion does not significantly impact the false positive rate (see Fig. 7). To more significantly decrease the false positive rate, it seems better to select N ≤ 13. However, doing so would result in a significant decrease in the detection capabilities (by more than 15% in terms of S/N loss). As a conclusion of this study, we choose N = 20 spectral priors for the IFS instrument. This number is driven by a trade-off between precision and recall, namely, in order to keep the S/N loss sufficiently small, while limiting the number of false positives to a value similar to the IRDIS one.

The whole optimization process for a given stellar spectral type and for a targeted number N of spectral priors is described in the form of a pseudo-code by Algorithm 1. This selection procedure was repeated for all stellar spectral types considered in this work. Figure 8 shows an example of such built library of spectral priors for a G2 star observed in YJ bands.

Based on the built library of spectral priors, we can now, as for IRDIS, compare the empirical false positive rate with the theoretical value. We can again use Eq. (1) with the number of processed pixels for IFS n_pixel = 140 000. We find N_fp = 0.04 for IFS. We can see with Fig. 7 that (as for IRDIS) the false positive rate when considering individual S/N map is in good agreement with what is expected from a Gaussian noise distribution.

Fig. 7 Number of false positives experienced at detection time when using between N = 1 and N = 31 spectral priors with IFS. The detection threshold is set at τ = 5, and the results are averaged over 32 observations. The red dashed line shows the theoretical false positive rate on a single S/N map. The mean experienced false positive rate for each S/N considered independently map is represented by the yellow dashed rectangles and the mean experienced cumulative false positive rate by the blue rectangles.

Fig. 8 Example of the library ΩIFS ∊ ℝ20×39 of spectral priors built for a G2 star in YJ bands. Priors are expressed in contrast unit and are normalized between 0 and 1. We systematically included in the library a "flat" spectral prior giving the same weights to all spectral channels.

4.3 Refined astrometric and photometric error budgets

The PACO algorithm provides only fitting errors for both photometric and astrometric measurements. To get a complete error budget, we need to take into account other sources of errors.

For the astrometry error budget, we use both the results of the F150 analysis with the SpeCal pipeline (Langlois et al. 2021) and the calibration obtained by Maire et al. (2021). However, the F150 analysis used an average value of the typical centering error. Thanks to the improved frame centering routine described in Sect. 3, we are able to derive a precise estimate of this error on each data set and we propagate it through the whole pipeline.

For the photometry error budget, we used, for the first time, the differential tip-tilt sensor (DTTS, Baudoz et al. 2010, see Sect. 4.3.2 for the details) measurements to derive proper photometric error bars. Although such a precise analysis has already been done on specific targets for some particular studies, nothing was implemented routinely to perform massive analysis. This new analysis allows us to perform a complete, accurate, automated, and homogeneous error estimation for both astrometry and photometry.

4.3.1 Astrometric error budget

The PACO algorithm provides an astrometric fitting error term, hereafter denoted σ_{sep, PACO} and σ_{PA, PACO}, associated, respectively, to the angular separation (sep) and to the parallactic angle (PA) of a given signal. Several additional sources of errors induced by pre-processing steps, such as the recenter-ing of the individual frames, or systematics related to SPHERE itself must be considered. We therefore combine several additional terms to refine the global error budget. The uncertainties associated with the separation and PA are found as follows. For the uncertainties on the angular separation, we combined four terms:

A distortion error of 0.4 mas at 1 as (Maire et al. 2021), scaling linearly with the separation of a source: (2)

A plate scale error, scaling linearly with the separation of a source. This plate scale and the associated error bars are measured during each observing run (astrometric calibration, see Langlois et al. (2021): (3)

An error on the re-centering of the individual frames, as estimated by the re-centering procedure described in Sect. 3: (4)

PACO internal error σ_sep,PACO.

Those 4 terms are quadratically combined to obtain the full error budget σ_sep,tot(as).

For the uncertainties on the PA, we also combine four terms: An error on the pupil angle equal to 0.52 mas at 1 as (Maire et al. 2021), scaling linearly with the separation: (5)

An error associated with the true North, as measured using the astrometric calibrations, σ_TN(°). An error on the re-centering of the individual frames, as estimated by the re-centering procedure (see Sect. 3): (6)

PACO internal error σ_PA,PACO(°).

Those four terms are also quadratically combined to obtain the full error budget σ_PA,tot(°).

4.3.2 Photometric error budget

Our aim here is to estimate the relative photometric uncertainties using SPARTA data (Suárez Valles et al. 2012) as well as using information from the DTTS. The DTTS is a control organ of SPHERE that ensures the star is always well centered on the coronagraph. It diverts a small fraction of the stellar light to produce an image of the star, allowing us to have a direct access to a PSF during the observation. While this PSF is not exactly the same as it would be on the science cameras due to non-common aberrations, it can still be used to monitor the photometric variability during the observing sequence. In the following, these series are denoted by DTTS(t) as a function of the time t of observation. SPARTA is the real time control computer of the adaptive optics system. Over the course of an observation, SPARTA collects information on the observing conditions that are then stored. As for the astrometry, the PACO algorithm provides a fitting photometric error σ_PACO, but additional terms are needed to estimate the global error budget, as follows.

First, there is the error associated with the flux calibration of the coronagraphic frames using the PSF. Because the observing conditions vary during the observing sequence, this error is time-dependent. Using datasets with bright background companions, detectable with a high S/N on each individual frame, we found that the flux variations are well correlated with the Strehl ratio (SR) variations as provided by SPARTA⁶. Besides, the photometry of our faint sources cannot be estimated on each frame. Hence, we used the time series SR(t) data taken during the observation to estimate an average photometric error over the whole coronagraphic sequence.

Because SPARTA provides the SR at 1.6 μm, we compute the Strehl ratio at the working wavelength of λ_i by using the following approximation – based solely on the adaptive optics fitting error and the Maréchal approximation (Maréchal 1948) – to capture the wavelength dependency for good conditions: (7)

Then, we computed the standard deviation of the SR, after removing the values for rejected frames, as follows: (8)

where SR₀ is the SR of the PSF used. In PACO, it is the average between the off-axis PSF taken before the observing corona-graphic sequence and the one after. The error associated to the extracted spectra at each wavelength λ_i is therefore: (9)

It sometimes happens that no SPARTA data are available. In such a case, the error is computed using the difference between the two available PSFs.

Second, there is the sky transparency: this term is measured thanks to the DTTS. While the peak of the DTTS PSF is also linked to the SR variation recorded by SPARTA, the total integrated flux from the DTTS PSF directly relates to the evolution of the sky transparency, but does not provide an absolute photometric measurement. By estimating the DTTS flux, we estimate the median and standard-deviation of the sky transparency variations over the coronagraphic sequence: (10) (11)

Lastly, we have the PACO internal error σ_PACO. Those three terms are then quadratically summed to obtain the full photometric error budget.

4.4 Multi-epoch discrimination tool

Second epoch observations are crucial to test whether a signal is due to a source gravitationally bound to the target star. We developed a tool that automatically performs this analysis. In cases where the candidate is not recovered, it computes the false positive probability, given the detection limits achieved for the second epoch.

In practice, we first checked whether a source is detected in the second epoch data at the position expected from a background object, knowing the proper motion of the star⁷. This step allows us to identify the background sources (providing that the detection limits of the second epoch data set are good enough). If the signal cannot be associated with a background source, we search for a fainter signal (detected at τ ≥ 4) within a disk D centered on the position of the source at the first epoch, and with a radius corresponding to the motion of a gravitionally bound object orbiting on a circular pole-on orbit (corresponding to the maximum possible motion in projected separation). If a signal is found, we attribute it to a possible companion. In case of an ambiguous choice (i.e., the motion of the object can be associated either with a background source or a gravitionally bound source) a flag is raised to report the ambiguity. If no signal is found, we check whether the non-detection in the second epoch is due to poorer conditions or to the fact that the first detection was a false positive. To do so, we carried out the following steps.

First, we compute the S/N (labeled as S/N₂) that the signal should have in the second epoch using the second epoch contrast⁸ map. Next, we measured the maximum S/N in the disk area in the second epoch data (S/N_max). Finally, as the S/N maps follow a centered Gaussian distribution with a unit variance, we can estimate the probability p of the source to be a real signal: (12)

5 Contrast comparison with other algorithms

Here, we present the contrast performance achieved by PACO on both IRDIS and IFS and compare them with the performance achieved by TLOCI (Marois et al. 2014) and PCA ADI (Soummer et al. 2012; Amara & Quanz 2012) for IRDIS, and TLOCI and PCA ASDI for IFS (Mesa et al. 2015). Those algorithms were used for the analysis of the SHINE F150 survey (Langlois et al. 2021). They are implemented in the SpeCal package (Galicher et al. 2018) dedicated to the analysis of the SHINE data with the following characteristics:

For TLOCI, the stellar profile is estimated frame by frame for each pixel in the field of view. The estimation uses a linear combination of all data to minimize the residuals after subtraction. The area on which the optimization is computed is much bigger than the subtraction area, thus mitigating as best as possible the self-subtraction of point-like sources. The Specal implementation of TLOCI also assumes a flat planet spectrum in contrast. The parameters for TLOCI were set as they were in the F150 reduction, as follows.

First, the optimization zone was separated by 0.5 full width at half maximum (FWHM) from the region of interest to avoid bias in the linear combination in case of the presence of a source in the region of interest. Next, the radial width (in radius) of the subtraction zone was set to one FWHM; Then, the radial to azimuthal ratio of the subtraction zone was set to 1.5; Finally, the optimization zone area was set to 20 PSF FWHM.

Both PCA algorithms are based on the equations described in Soummer et al. (2012). For IRDIS (ADI processing), the principal components (PCs) are computed independently for each spectral channel. For IFS (ASDI processing), the PCs are computed using the spatial and spectral dimensions simultaneously. For IRDIS, 5 PCs were used, and for IFS, 50, 100 and, 150 PCs were used.

The throughput of both TLOCI and PCA was estimated internally by SpeCal at each location in the field of view by generating a datacube with fake planets. The 1D throughput curve is then applied to the residual maps of both algorithms.

The contrast curves and maps were estimated for each spectral channel by computing the pixel by pixel azimuthal standard deviation in an annulus of 0.5 FWHM on the residual maps, once it had been corrected from the throughput. The 5σ detection limits are derived from this estimation by taking into account several corrections: the flux loss from ADI subtraction, the coronograph transmission, and the neutral density of the off-axis PSF. Finally, these detection limits were normalized by the off-axis PSF flux. The S/N maps were directly derived from the estimated flux and its associated standard-deviation (i.e., contrast at 1σ).

In Sect. 5.1, we present the contrast comparison between PACO and the algorithms described above. Since this direct comparison of contrast is biased by the diverse hypotheses made by each algorithm, we then present in Sect. 5.2 a set of numerical experiments resorting to massive injections of FPSs. This demonstrates the reliability of the contrast curves obtained with PACO, as well as the gain in sensitivity and control of the probability of false alarms compared with the two other algorithms.

5.1 Contrast performance

Figure 9 shows the predicted 5σ contrast as a function of the angular separation on IRDIS (with Fig. 10 showing the same but for IFS) obtained with the three considered algorithms for all epochs of the 24 stars considered in this study. For this comparison, we used the combined contrast of TLOCI and PCA, as well as the contrast obtained with the flat prior with PACO (hereafter noted as PACO-flat) to make the most direct comparison possible. Although it was not used in this study, we also included the median PACO ADI contrast curve (i.e., obtained without joint processing of the spectral channels) in Fig. 9 for reference. The gain offered by the ASDI mode is very close , which corresponds to the expected theoretical value when combining the information of two (independent) channels. In order to compare results from the three considered algorithms, it has to be noted that given the non-statistical nature of TLOCI and PCA, the 5σ detection limits are not statistically grounded; that is to say, we experienced in practice many more false alarms as theoretically expected for the targeted confidence level, especially at short angular separations from the target star. Moreover, the flat prior represents the most difficult case for PACO (especially for IFS) because we try to detect a signal with the same spectra as the host star, which means that the spectral prior does not explicitly help in disentangling the two components.

For IRDIS, at close angular separations, PACO shows a better performance than TLOCI (resp. PCA) by a factor of about 7 (resp. 5) at 0.5'', and by 5 at 1'' and beyond compare to both TLOCI and PCA. Moreover, due to the statistical nature of PACO, the number of false positives follows what is expected at a 5σ confidence under a multivariate Gaussian hypothesis (see Fig. 6), unlike TLOCI and PCA.

For IFS, the gain ranges between a factor of 3 and 5, depending on the separation, compared to PCA. It is much larger compared to TLOCI by about a factor of 10 for all separations, but this result is expected as TLOCI performs worse on IFS as compared to PCA. However, the achieved performance in terms of contrast is much more consistent with PACO than it is with PCA. We remind that using the flat prior as a benchmark allows us to compare PACO with PCA/TLOCI as fairly as possible. It does not however represent the full capability of PACO to detect faint sources because this is the "worst" possible case, as we are trying to detect a highly correlated planetary spectrum with respect to the star spectrum.

Fig. 9 Contrast comparison at 5σ between PACO, TLOCI, and PCA for IRDIS for the 24 science stars considered in this paper. A flat prior was used for PACO. Dashed lines show the 95% completeness interval. The grey area represent the coronographic mask. Contrast curves provided by PCA and TLOCI do not strictly correspond to a 5σ false alarm rate contrarily to the contrast curves of PACO. The achievable contrasts are thus significantly over-optimistic for PCA and TLOCI (see discussion in the main text).

Fig. 10 Contrast comparison at 5σ between PACO, TLOCI, and PCA for IFS for the 24 science stars considered in this paper. Dashed lines show the 95% completeness interval. The grey area represents the coronographic mask. Contrast curves provided by PCA and TLOCI do not correspond to a 5σ false alarm rate contrarily to the contrast curves of PACO. The achievable contrasts are thus significantly over-optimistic for PCA and TLOCI (see discussion in the text).

Fig. 11 Comparison of the contrast curves at 5σ obtained with PCA (5 modes), TLOCI, and PACO-flat for IRDIS on HD 377. The contrasts of the injected fake planets were computed using 2D contrast maps, thus, there are differences seen with the 5σ curve: local variations of the achieved contrast are averaged azimuthally.

5.2 Validation of the reliability of the contrast curves

Each algorithm uses different hypothesis to compute the contrast limits. To further assess the gain in contrast obtained with PACO, and the comparison with TLOCI or PCA, we present a set of numerical experiments. With PACO, the contrast is estimated assuming that the statistical parameters (mean and covariance matrices) characterizing the stellar leakages are computed from pure noise realizations. In practice, the underlying presence of the (unknown) sought objects corrupts these estimates, which leads to a (slight) bias in the estimated performance. In this section, we aim to quantify this bias via numerical experiments. For that purpose, we resorted to massive injections of FPSs at contrast levels predicted by PACO and we re-ran the algorithm to quantify the real detection confidence experienced for such levels of contrast. Ensuring such confidence is key to allow an unsupervised selection of candidate companions through simple thresholding of the derived S/N map. Given computational constraints, injection tests (as presented in Sects. 5.2.1 and 5.2.2) were performed on a dataset of HD 377 (star included in this study) for both IRDIS and IFS.

Fig. 12 S/N map provided by PACO (top right), residual map provided by SpeCal PCA using 10 modes (bottom middle), associated S/N map (top middle), residual map from TLOCI (bottom left) and associated S/N map (top left). The injected fake planets are clearly visible on the S/N maps from PACO. None of the high S/N detections on the PCA/TLOCI map corresponds to injected sources. The locations of injected sources are highlighted by the red boxes.

Fig. 13 Retrieved S/N for the injected sources for IRDIS. The detection threshold τ was set to 4. The red line shows the median S/N of the injected sources and the grey area the 1σ interval containing 68% of the detected injected sources.

5.2.1 IRDIS data

For IRDIS, two sets were created: one with injections at close separations (<1.25 as) and one at large separations (>1.5 as), as described in Appendix B. For both, from the 2D contrast maps provided by PACO, we set the contrast of the FPSs in order to achieve (theoretically) a detection slightly above the 5σ threshold on the S/N map. More detailed information on the close-in injected sources parameters can be found in Table B.1. Figure 11 shows the injected sources contrast (red dots) compared to the contrast curve provided by PACO, as well as the contrast curves provided by TLOCI and PCA for this particular target in the first case (close separation). Figure 12 shows the S/N maps obtained with the three algorithms as well as residual maps⁹ obtained after the subtraction of the estimated stellar component for PCA and TLOCI for the injections at close separations.

We also created a set of injections of companions at larger separations. Detailed information on the injected sources parameters can be found in Table B.2. Figure B.1 shows the contrast of the injected FPSs compared to the contrast curve of the three algorithms. The corresponding S/N and residual maps obtained after subtraction of the estimated stellar component for PCA and TLOCI are shown in Figs. B.2 to B.6.

Figure 13 shows the retrieved S/N of FPSs in both cases (for a total of 90 injected sources)¹⁰. The median S/N is 5.0 which is in very good agreement with the contrast of injection. We also calculated the empirical standard deviation value σ = 1.15. This result is also close to the theoretical value of 1 predicted by a Gaussian distribution.

These results confirm that the contrast estimates produced by PACO (e.g., shown in Fig. 9) are reliable and statistically grounded, and that the detection sensitivity is improved compared to TLOCI and PCA for the whole range of angular separations.

Fig. 14 S/N and residual maps for the injected fake sources on IFS. Top row: PACO ASDI S/N maps corresponding to the various injected spectra (see Fig. B.7). Middle row: Residual and S/N maps using PCA. Bottom row: Residual and S/N maps using TLOCI. The locations of injected sources are highlighted by red boxes.

5.2.2 IFS data

For IFS, we injected 72 sources with three different shapes of spectra (24 sources per shape) corresponding to three of the spectral priors (flat, L-type, T-type) used during the reduction (hereafter denoted as cases), for which the contrast injected was also computed using the 2D contrast maps. We also considered a "flat" injection that corresponds (as the priors are expressed in contrast unit) to a SED with a same shape than that of the star. The normalized spectra injected for the three cases can be found in Fig. B.7 as well as complete information on the injected sources in Table B.3. Corresponding S/N and residual maps are given in Fig. 14.

The S/N of the detected injections are shown in Fig. 15. As done for IRDIS, we can compute the median S/N of the injected sources, which is 4.3, with a standard deviation of 1.16. This suggests that the estimated contrast is slightly optimistic by about 15%.

As expected, we conclude that the detection limits in contrast derived by PACO for IFS are slightly optimistic without questioning the previously mentioned results, because the equivalent number of independent spectral channels that are recombined is about L/2. The detection sensitivity is still improved at all angular separation, coupled with the false positive rate consistent with the chosen S/N threshold.

Fig. 15 Retrieved S/N for the injected sources for IFS. The detection threshold τ was set to 4. The red line shows the median S/N of the injected sources and the grey area the 1σ interval containing 68% of the detected injected sources.

6 Results of the mini-survey

6.1 Identified point sources and status

Running the PACO algorithm along with the analysis tools described in Sects. 4.2 and 4.4 over the 37 datasets of this survey allows us to identify 61 (58 IRDIS, 3 IFS) point-like sources with an S/N above the 5σ detection confidence. We note that one source, PZ TEL B, is detected with both instruments. For comparison, only 40 sources were detected by the F150 analysis. This again illustrates the enhanced detection capabilities of PACO with respect to more classical algorithms. We classify the sources as follows:

BCKG (background source): a source classified as background in the F150 or in this analysis using a proper motion analysis.

KC (known companion): either a planet, brown dwarf or stellar companion.

SUSP BCKG CMD (suspected background using a color magnitude diagram): source detected in the present analysis, which was not detected in the F150, with only one epoch available and a color consistent with a background nature.

CC (candidate companion): a source detected in the present analysis that was not detected in the F150, with only one epoch and with a color compatible with a planetary or brown dwarf companion.

FP: a false positive identified using the multi-epoch tool described in Sect. 4.4. This classification is only possible for stars with multiple epochs with similar quality data.

Figure 16 shows a pie-chart diagram of the sources classification. Individual data are to be provided in a VizieR table (see Table 1 for an example of the parameters provided). Finally, the targets are plotted on a CMD in Fig. 17.

Among the 60 identified signal of interest, eight are not characterized because of a lack of a second epoch of observation. Among those eight sources, six fall into the domain of probable background sources (grey area in Fig. 17). We classified them as SUSP BCKG CMD. Two are more promising giving their colors. They require additional observations to definitely distinguish between a background source, a false positive, or a bound companion. We therefore classified them as candidate companions (CC) at this stage.

Finally, we found twelve false positives over the course of this survey. Given the number of datasets for IRDIS (37) and the number of priors used (5), we should expect around 16 false positives (see Fig. 6). Taking also into account that FPs could be present amongst the eight sources classified as candidates companions (either "CC" or "SUSP BCKG CMD"), we conclude that this number is in good agreement with the theory predictions, confirming the Gaussian nature of noise in the detection maps of PACO. One IFS false positive was also found.

Fig. 16 Classification of the 60 sources detected during the survey. Note: PZ Tel B is visible in both IRDIS and IFS, thus accounted for only one source in this plot.

Fig. 17 CMD plot with the classification of the 60 sources detected during the survey. Note: individual error bars are not shown for clarity purposes except for the two most interesting candidates (yellow circles). Typical error bars are indicated in the top right. The grey dots represent the background sources identified during the SHINE F150 survey.

Fig. 18 Detection maps using the MESS2 (Lannier 2015) tool for HD 202917, DI only (left), DI+RV (right). The eccentricity is assumed between 0 and 0.5 and the inclination between 0 and 90 degrees. Contrast maps are converted to mass maps using the COND atmospheric model (Allard et al. 2001) and the stellar parameters of the system are extracted from Desidera et al. (2021). Two epochs were used here.

Fig. 19 Comparison between the detection limits obtained with PACO ASDI (left) and TLOCI (right) for HD 202917 using the MESS2 tool (Lannier 2015). The eccentricity is assumed between 0 and 0.5 and the inclination between 0 and 90 degrees. Contrast maps are converted to mass maps using the COND atmospheric model (Allard et al. 2001) and the stellar parameters of the system are extracted from Desidera et al. (2021). Two epochs were used.

6.2 Searching for additional companions

Here, we want to put constraints on the properties of possible, yet unseen companions. We used the MESS2 tool (Lannier 2015) that uses, for each target, all the detection limit maps derived from the PACO analysis (expressed in terms of mass) as well as the radial velocity data, whenever they are available in the ESO HARPS (Mayor et al. 2003) archive. Indeed, the combination of direct imaging data with RV data provides a large exploration of the star's environments, going from a fraction of au out to 100 au. To convert the PACO contrasts into detection limits expressed in masses, we use luminosity-mass relationships given by the COND atmospheric model (Allard et al. 2001), and the stars' ages and masses provided by Desidera et al. (2021). Finally, to run MESS2, we assumed a uniform distribution of eccentricities between 0 and 0.5 and a uniform orbital plane inclination between 0 and 90 degrees.

An example of detection limits obtained using the imaging data alone, on the one hand, and using both the imaging and RV data, on the other hand, is provided in Fig. 18. The results obtained for all stars of the sample are presented in Appendix C.

The detection limits in terms of masses are significantly improved with respect to previous analyses, thanks to improved detection limits. For instance, we get a detection limit of about 5 M_Jup (68% probability) at 5 au for HD 202917, while the detection limit with TLOCI is 10 M_Jup at the same separation (see Fig. 19). At 10 au, 3 MJup planets could be detected, compared to 6 MJup with TLOCI. This represents a substantial improvement in the detection limits. Nonetheless, Jupiter siblings are still out of reach in the present data. Improved adaptive optics systems, such as that of the SPHERE+ project (Boccaletti et al. 2020) on the VLT, are needed to reach such objectives. Finally, we see that, in most cases, combining the radial velocity (RV) and direct imaging (DI) allows us to bridge the gap between the two techniques. The detection limits will be further improved using HIPPARCOS-Gaia data in the typical 3–10 au region.

7 Conclusion

Here, we present a number of upgrades made to our reduction and analysis pipeline using PACO that are aimed at improving the sensitivity and the characterization of the detected sources. In particular, we have improved the precision and robustness of the astrometric and photometric error bars for the detected sources. We have developed custom built spectral prior libraries to optimize the detection capability of the ASDI mode for both IRDIS and IFS. The contrast performances are significantly improved, compared to what can be obtained with classical algorithms such as TLOCI and PCA. Also, we have shown that PACO provides statistically meaningful S/N maps. This work paves the way toward an end-to-end, homogeneous, and unsupervised massive re-reduction of archival SPHERE direct imaging data in the quest for detecting exoJupiters.

We used PACO ASDI to search for exoJupiters in a sample of 24 selected young, solar-type targets observed with SPHERE IRDIS and IFS that are part of the SHINE survey. This new analysis allowed us to identify two candidate companions in this small sample. Second epochs are necessary to unveil their nature.

Acknowledgements

This project has also received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (COBREX; grant agreement no. 885593). SPHERE is an instrument designed and built by a consortium consisting of IPAG (Grenoble, France), MPIA (Heidelberg, Germany), LAM (Marseille, France), LESIA (Paris, France), Laboratoire Lagrange (Nice, France), INAF - Osservatorio di Padova (Italy), Observatoire de Genève (Switzerland), ETH Zürich (Switzerland), NOVA (Netherlands), ONERA (France) and ASTRON (Netherlands) in collaboration with ESO. SPHERE was funded by ESO, with additional contributions from CNRS (France), MPIA (Germany), INAF (Italy), FINES (Switzerland) and NOVA (Netherlands). SPHERE also received funding from the European Commission Sixth and Seventh Framework Programmes as part of the Optical Infrared Coordination Network for Astronomy (OPTICON) under grant number RII3-Ct-2004-001566 for FP6 (2004-2008), grant number 226604 for FP7 (2009-2012) and grant number 312430 for FP7 (2013-2016). This work has made use of the SPHERE Data Centre, jointly operated by OSUG/IPAG (Grenoble), PYTHEAS/LAM/CeSAM (Marseille), OCA/Lagrange (Nice), Observatoire de Paris/LESIA (Paris), and Observatoire de Lyon (OSUL/CRAL). This research has made use of the SIMBAD database and VizieR catalogue access tool, operated at CDS, Strasbourg, France. This work is supported by the French National Research Agency in the framework of the Investissements d'Avenir program (ANR-15-IDEX-02), through the funding of the "Origin of Life" project of the Univ. Grenoble-Alpes. This work was supported by the Action Spécifique Haute Résolution Angulaire (ASHRA) of CNRS/INSU co-funded by CNES. This work is based on observations collected at the European Southern Observatory under ESO programme(s) 096.C-0241(A/B/C/F/G), 095.C-0298(A/B/D/H), 1100.C-0481(F/G), 097.C-0865(A/B/C/D), 198.C-0209(A/E/G).

Appendix A Target parameters and observing logs

Appendix B Compared detection and contrast maps between PACO ASDI, TLOCI and PCA for 5σ injected sources from PACO detection limits

Fig. B.1 Contrast comparison between PCA (10 modes), TLOCI and PACO. The contrast of the injected fake planets were computed using 2D contrast maps, hence, we can see the differences with the 5σ curve: local variations of the achieved contrast are averaged azimuthally.

Fig. B.2 PACO ASDI S/N map assuming the SED of sought sources is flat. The locations of injected sources are highlighted by the red boxes.

Fig. B.3 PCA residual map. The locations of injected sources are highlighted by red boxes.

Fig. B.4 PCA S/N map. The locations of injected sources are highlighted by red boxes.

Fig. B.5 TLOCI residual map. The locations of injected sources are highlighted by red boxes.

Fig. B.6 TLOCI S/N map. The locations of injected sources are highlighted by red boxes.

Fig. B.7 Normalized injected contrasts for IFS for the three cases considered.

Appendix C MESS2 results

Fig. C.1 Results for HD 105, DI only (left), DI+RV (right). One epoch was available for DI.

Fig. C.2 Results for HD 377, DI only (left), DI+RV (right). One epoch was available for DI.

Fig. C.3 Results for HD 987, DI only (left), DI+RV (right). One epoch was available for DI.

Fig. C.4 Results for HD 1466, DI only (left), DI+RV (right). Two epochs were available for DI.

Fig. C.5 Results for BD-12 243, DI only (left), DI+RV (right). One epoch was available for DI.

Fig. C.6 Results for HD 8558, DI only (left), DI+RV (right). Three epochs were available for DI.

Fig. C.7 Results for HD 17925, DI only (left), DI+RV (right). One epoch was available for DI.

Fig. C.8 Results for HD 43989, DI only (left), DI+RV (right). Two epochs were available for DI.

Fig. C.9 Results for HD 44627, DI only (left), DI+RV (right). One epoch was available for DI.

Fig. C.10 Results for HD 45270, DI only (left), DI+RV (right). One epoch was available for DI.

Fig. C.11 Results for CD-61 1439, DI only (left), DI+RV (right). One epoch was available for DI.

Fig. C.12 Results for HD 49855, DI only (left), DI+RV (right). One epoch was available for DI.

Fig. C.13 Results for HD 61005, DI only (left), DI+RV (right). One epoch was available for DI.

Fig. C.14 Results for HD 118100, DI only (left), DI+RV (right). One epoch was available for DI.

Fig. C.15 Results for HD 1642490, DI only (left), DI+RV (right). Four epochs were available for DI.

Fig. C.16 Results for HD 174429, DI only (left), DI+RV (right). Five epochs were available for DI.

Fig. C.17 Results for HD 181327, DI only (left), DI+RV (right). One epoch was available for DI.

Fig. C.18 Results for HD 189245, DI only (left), DI+RV (right). One epoch was available for DI.

Fig. C.19 Results for HD 218860, DI only (left), DI+RV (right). One epoch was available for DI.

Fig. C.20 Results for HD 224228, DI only (left), DI+RV (right). One epoch was available for DI.

Fig. C.21 Results for HD 90712 (left) and HD 197890 (right). No radial velocity data were available. One epoch was available for DI.

Fig. C.22 Results for CD-31 16041. No radial velocity data were available. Two epochs were available for DI.

References

All Tables

All Figures

	Fig. 1 Properties of the sample stars in terms of distance, age, spectral type, and magnitude.
In the text

Fig. 2 Evolution of the S/N (blue) for various priors for a point source observed with IRDIS, and whose spectrum corresponds to the red mark. The x-axis gives various values of the priors. The S/N corresponding to the [0.5, 1] prior (significantly different from the normalized contrast of the injected source) is below the 4σ threshold used for this analysis, while the S/N reached for spectral prior [1, 0.1] (similar to the normalized contrast of the injected source) leads to a clear detection with a significance above the 5σ detection threshold (dashed black line).

In the text

Fig. 3 Evolution of the S/N (blue) for various priors for HD 206893b (2017-07-13, IFS, YJ band). The x-axis gives the index of the priors. The S/N corresponding to the 10 last priors are below the 5σ detection threshold. Priors are ordered by decreasing S/N for clarity purposes. The S/N reached for spectral prior close to the spectrum (in contrast) leads to a detection with a significance above the 5σ detection threshold.

In the text

	Fig. 4 Integrated H2/H3 flux ratio computed for a planet with a solar metallicity, C/O=0.50, and with different values of log(g) and Teff. ExoREM spectra were used.
In the text

	Fig. 5 Integrated K1/K2 flux ratio computed for a planet with a solar metallicity, C/O = 0.50, and with different values of log(g) and Teff. ExoREM spectra were used.
In the text

Fig. 6 Average number of false positives for a detection threshold τ = 5 with one to eight spectral priors with IRDIS. Twenty observations were considered for these experiments. The red dashed line shows the theoretical false positive rate on a single S/N map. The mean experienced false positive rate for each S/N considered independently map is represented by the yellow dashed rectangles and the mean experienced cumulative false positive rate by the blue rectangles.

In the text

Fig. 7 Number of false positives experienced at detection time when using between N = 1 and N = 31 spectral priors with IFS. The detection threshold is set at τ = 5, and the results are averaged over 32 observations. The red dashed line shows the theoretical false positive rate on a single S/N map. The mean experienced false positive rate for each S/N considered independently map is represented by the yellow dashed rectangles and the mean experienced cumulative false positive rate by the blue rectangles.

In the text

	Fig. 8 Example of the library ΩIFS ∊ ℝ20×39 of spectral priors built for a G2 star in YJ bands. Priors are expressed in contrast unit and are normalized between 0 and 1. We systematically included in the library a "flat" spectral prior giving the same weights to all spectral channels.
In the text

Fig. 9 Contrast comparison at 5σ between PACO, TLOCI, and PCA for IRDIS for the 24 science stars considered in this paper. A flat prior was used for PACO. Dashed lines show the 95% completeness interval. The grey area represent the coronographic mask. Contrast curves provided by PCA and TLOCI do not strictly correspond to a 5σ false alarm rate contrarily to the contrast curves of PACO. The achievable contrasts are thus significantly over-optimistic for PCA and TLOCI (see discussion in the main text).

In the text

Fig. 10 Contrast comparison at 5σ between PACO, TLOCI, and PCA for IFS for the 24 science stars considered in this paper. Dashed lines show the 95% completeness interval. The grey area represents the coronographic mask. Contrast curves provided by PCA and TLOCI do not correspond to a 5σ false alarm rate contrarily to the contrast curves of PACO. The achievable contrasts are thus significantly over-optimistic for PCA and TLOCI (see discussion in the text).

In the text

	Fig. 11 Comparison of the contrast curves at 5σ obtained with PCA (5 modes), TLOCI, and PACO-flat for IRDIS on HD 377. The contrasts of the injected fake planets were computed using 2D contrast maps, thus, there are differences seen with the 5σ curve: local variations of the achieved contrast are averaged azimuthally.
In the text

Fig. 12 S/N map provided by PACO (top right), residual map provided by SpeCal PCA using 10 modes (bottom middle), associated S/N map (top middle), residual map from TLOCI (bottom left) and associated S/N map (top left). The injected fake planets are clearly visible on the S/N maps from PACO. None of the high S/N detections on the PCA/TLOCI map corresponds to injected sources. The locations of injected sources are highlighted by the red boxes.

In the text

	Fig. 13 Retrieved S/N for the injected sources for IRDIS. The detection threshold τ was set to 4. The red line shows the median S/N of the injected sources and the grey area the 1σ interval containing 68% of the detected injected sources.
In the text

	Fig. 14 S/N and residual maps for the injected fake sources on IFS. Top row: PACO ASDI S/N maps corresponding to the various injected spectra (see Fig. B.7). Middle row: Residual and S/N maps using PCA. Bottom row: Residual and S/N maps using TLOCI. The locations of injected sources are highlighted by red boxes.
In the text

	Fig. 15 Retrieved S/N for the injected sources for IFS. The detection threshold τ was set to 4. The red line shows the median S/N of the injected sources and the grey area the 1σ interval containing 68% of the detected injected sources.
In the text

	Fig. 16 Classification of the 60 sources detected during the survey. Note: PZ Tel B is visible in both IRDIS and IFS, thus accounted for only one source in this plot.
In the text

	Fig. 17 CMD plot with the classification of the 60 sources detected during the survey. Note: individual error bars are not shown for clarity purposes except for the two most interesting candidates (yellow circles). Typical error bars are indicated in the top right. The grey dots represent the background sources identified during the SHINE F150 survey.
In the text

Fig. 18 Detection maps using the MESS2 (Lannier 2015) tool for HD 202917, DI only (left), DI+RV (right). The eccentricity is assumed between 0 and 0.5 and the inclination between 0 and 90 degrees. Contrast maps are converted to mass maps using the COND atmospheric model (Allard et al. 2001) and the stellar parameters of the system are extracted from Desidera et al. (2021). Two epochs were used here.

In the text

Fig. 19 Comparison between the detection limits obtained with PACO ASDI (left) and TLOCI (right) for HD 202917 using the MESS2 tool (Lannier 2015). The eccentricity is assumed between 0 and 0.5 and the inclination between 0 and 90 degrees. Contrast maps are converted to mass maps using the COND atmospheric model (Allard et al. 2001) and the stellar parameters of the system are extracted from Desidera et al. (2021). Two epochs were used.

In the text

	Fig. B.1 Contrast comparison between PCA (10 modes), TLOCI and PACO. The contrast of the injected fake planets were computed using 2D contrast maps, hence, we can see the differences with the 5σ curve: local variations of the achieved contrast are averaged azimuthally.
In the text

	Fig. B.2 PACO ASDI S/N map assuming the SED of sought sources is flat. The locations of injected sources are highlighted by the red boxes.
In the text

	Fig. B.3 PCA residual map. The locations of injected sources are highlighted by red boxes.
In the text

	Fig. B.4 PCA S/N map. The locations of injected sources are highlighted by red boxes.
In the text

	Fig. B.5 TLOCI residual map. The locations of injected sources are highlighted by red boxes.
In the text

	Fig. B.6 TLOCI S/N map. The locations of injected sources are highlighted by red boxes.
In the text

	Fig. B.7 Normalized injected contrasts for IFS for the three cases considered.
In the text

	Fig. C.1 Results for HD 105, DI only (left), DI+RV (right). One epoch was available for DI.
In the text

	Fig. C.2 Results for HD 377, DI only (left), DI+RV (right). One epoch was available for DI.
In the text

	Fig. C.3 Results for HD 987, DI only (left), DI+RV (right). One epoch was available for DI.
In the text

	Fig. C.4 Results for HD 1466, DI only (left), DI+RV (right). Two epochs were available for DI.
In the text

	Fig. C.5 Results for BD-12 243, DI only (left), DI+RV (right). One epoch was available for DI.
In the text

	Fig. C.6 Results for HD 8558, DI only (left), DI+RV (right). Three epochs were available for DI.
In the text

	Fig. C.7 Results for HD 17925, DI only (left), DI+RV (right). One epoch was available for DI.
In the text

	Fig. C.8 Results for HD 43989, DI only (left), DI+RV (right). Two epochs were available for DI.
In the text

	Fig. C.9 Results for HD 44627, DI only (left), DI+RV (right). One epoch was available for DI.
In the text

	Fig. C.10 Results for HD 45270, DI only (left), DI+RV (right). One epoch was available for DI.
In the text

	Fig. C.11 Results for CD-61 1439, DI only (left), DI+RV (right). One epoch was available for DI.
In the text

	Fig. C.12 Results for HD 49855, DI only (left), DI+RV (right). One epoch was available for DI.
In the text

	Fig. C.13 Results for HD 61005, DI only (left), DI+RV (right). One epoch was available for DI.
In the text

	Fig. C.14 Results for HD 118100, DI only (left), DI+RV (right). One epoch was available for DI.
In the text

	Fig. C.15 Results for HD 1642490, DI only (left), DI+RV (right). Four epochs were available for DI.
In the text

	Fig. C.16 Results for HD 174429, DI only (left), DI+RV (right). Five epochs were available for DI.
In the text

	Fig. C.17 Results for HD 181327, DI only (left), DI+RV (right). One epoch was available for DI.
In the text

	Fig. C.18 Results for HD 189245, DI only (left), DI+RV (right). One epoch was available for DI.
In the text

	Fig. C.19 Results for HD 218860, DI only (left), DI+RV (right). One epoch was available for DI.
In the text

	Fig. C.20 Results for HD 224228, DI only (left), DI+RV (right). One epoch was available for DI.
In the text

	Fig. C.21 Results for HD 90712 (left) and HD 197890 (right). No radial velocity data were available. One epoch was available for DI.
In the text

	Fig. C.22 Results for CD-31 16041. No radial velocity data were available. Two epochs were available for DI.
In the text