Preparing an unsupervised massive analysis of SPHERE high contrast data with the PACO algorithm

We aim at searching for exoplanets on the whole ESO/VLT-SPHERE archive with improved and unsupervised data analysis algorithm that could allow to detect massive giant planets at 5 au. To prepare, test and optimize our approach, we gathered a sample of twenty four solar-type stars observed with SPHERE using angular and spectral differential imaging modes. We use PACO, a new generation algorithm recently developed, that has been shown to outperform classical methods. We also improve the SPHERE pre-reduction pipeline, and optimize the outputs of PACO to enhance the detection performance. We develop custom built spectral prior libraries to optimize the detection capability of the ASDI mode for both IRDIS and IFS. Compared to previous works conducted with more classical algorithms than PACO, the contrast limits we derived are more reliable and significantly better, especially at short angular separations where a gain by a factor ten is obtained between 0.2 and 0.5 arcsec. Under good observing conditions, planets down to 5 MJup, orbiting at 5 au could be detected around stars within 60 parsec. We identified two exoplanet candidates that require follow-up to test for common proper motion. In this work, we demonstrated on a small sample the benefits of PACO in terms of achievable contrast and of control of the confidence levels. Besides, we have developed custom tools to take full benefits of this algorithm and to quantity the total error budget on the estimated astrometry and photometry. This work paves the way towards an end-to-end, homogeneous, and unsupervised massive re-reduction of archival direct imaging surveys in the quest of new exoJupiters.


Introduction
Since the discovery of the first, giant planets around solar type stars almost 30 years ago, thousand of planets have been discovered, with masses down to a few Earth-mass. Yet, solar system analogues 1 have not been detected, and we therefore still do not know if our own planetary system is unique or not. Detecting Earth twins is still not possible, and 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, and long term stellar activity is also responsible for a long term noise. As a consequence, the orbital parameters of thestill rare -RV planets announced beyond 5 au are poorly characterized (Wittenmyer et al. (2016), Fernandes et al. (2019), Fulton et al. (2021)). These limitations prevent precise comparisons between the radial distribution of giant 1 defined as planetary systems around solar-type stars, hosting inner Earth-mass planets with at least one in the habitable zone, and outer sub-Jupiter/Jupiter masses planets. planets beyond their forming regions and predictions from population synthesis models to constrain formation scenarios (Lagrange et al, 2022, submitted).
Giant Planets (GPs) 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 protoplanetary disk has dissipated and the dynamics is no longer controlled by gas. GPs could even play a role in the development of life on Earth analogues (Horner & Jones 2010), and could 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 model individual systems.
Detecting giant planets is therefore crucial to understand planetary system formation and evolution. While RV or transit techniques are best suited to detect GPs orbiting Article number, page 1 of 35 arXiv:2305.08766v1 [astro-ph.EP] 15 May 2023 A&A proofs: manuscript no. main within typically 5 au, they are not well adapted to detect and characterize 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 to detect and characterize analogues of our solar system 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 to detect a 4 M Jup planet at 5 au (Vigan et al. 2021). These limited performances are due to 1) limitations at the instrumental and data processing levels, and 2) a geometrical effect: unless on pole-on orbits, short period planets may be missed in a single observation because of a small projected separations at the time of the observation. For instance, due to geometrical effect, a single observation allows us to explore only 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 twice the star a few years apart (Lannier et al. 2017).
In this paper, we apply PACO (patch covariance), a promising detection algorithm (Flasseur et al. (2018); Flasseur et al. (2020a); Flasseur et al. (2020b) on a small sample of stars representative of SPHERE targets, 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. We moreover define our sample so as to address the following astrophysical question: how far are we from detecting solar system young giant planets siblings?
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 finally provided in Sect. 7.

Star sample
Our sample gathers 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 already analyzed with conventional postprocessing algorithms in Langlois et al. (2021), hereafter referred to as 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, i.e., at the locations of the solar system giants 2 . 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 like AB Doradus (ABDO,150 Myr), 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 3 . Table A.2 and Table A.3 provides 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 (at least 30 degrees). The atmospheric conditions were heterogeneous, with a seeing ranging from 0.5 to 1. most of the time, and the coherence time ranging between 1 to 10 ms.

Data reduction and frame centering
The reduction pipeline from raw to centered datasets is similar to the one built for the SHINE survey and 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., submitted). As part of the observing sequence, before and after the coronagraphic 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, i.e., more than one order of magnitude faster than the current pipeline, without any loss in precision.
Note that when precise astrometric measurements of known companions are needed, the satellite spots are generated during the whole coronagraphic sequence. In such a

Analysis pipeline
The analysis pipeline is based on the PACO A(S)DI pipeline described in Flasseur et al. (2020a) and Flasseur et al. (2020b). A few important upgrades were made, though: -Improvement of PACO robustness, i.e., capability to run PACO on diverse and heterogeneous datasets while consistently providing reliable results, see Sect. 4.1; -Optimization of spectral priors for PACO ASDI, see Sect. 4.2; -Automated and improved computation of astrometric and photometric error bars for each characterized source, see Sect. 4.3; -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) needed for further analysis. These upgrades are described in the following. 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.

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: -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-tonoise ratio (S/N) in the measured off-axis PSF (e.g., in absorption bands). -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.
Besides, an engineering effort has been made to validate through numerical experiments the faithfulness of the astrophysical quantities produced by the algorithm, in particular concerning the astrometry and photometry, and 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 will not be 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.

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 as 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, and, in Fig. 3 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 implemented in the PCA and TLOCI versions of the SPHERE data center, is somewhat similar 4 to the flat prior combination (i.e., assuming that the sought exoplanets 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 5 , 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. . 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.

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/characterization algorithms. The targets properties are provided in Table A.4, and the observing and atmospheric conditions are given in 4 In practice, PACO ASDI also accounts for a confidence weight estimated locally for each spectral channel, giving more weight to the spectral channels where the variance of the estimated flux is the lowest. This information is never accounted for in classical algorithms like TLOCI and PCA. 5 It has been shown in Flasseur et al. (2020b) that the S/N of detection is only marginally degraded in the case where the prior SEDs differ significantly from the true SED of the sources. In any case, the probability of false alarms remains controlled at the prescribed detection threshold τ .   (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.

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 H2 H3 (or K1 K2 ) where H2 and H3 (respectively K1 and K2) represent the integrated fluxes in these spectral bands. Figure 4 and Fig. 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 H2 H3 (resp K1 K2 ) ratio varies from several tens even hundreds down to 1 mostly because of strong CH 4 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 8 spectral priors uniformly spread over the considered parameter space (by steps of 0.2), 5 (by steps of 0.3), and 4 (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 can not be identified on a single frame), using 5 priors (or more) does not significantly degrade the S/N (less than 1%), while using 4 (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 τ : As an illustration, the number of pixels 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 −7 ), we have: N fp 10 6 × 2.87 × 10 −7 ∼ 0.287 , 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 choose to include five spectral priors in our library Ω IRDIS ∈ R 2×5 , achieving the desired trade-off between maximizing our detection performance and lowering the number of false positives:

Priors selection for IFS
In this section, we describe how we build a library Ω IFS ∈ R 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 ∈ R 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 purpose of illustration, the method is described for any given stellar spectral type without loss of generality. The set Ω IFS is built from a set Ω ER ∈ R N ER ×L ER of substellar spectra provided by the ExoREM models. The spectral resolution L ER of each element w ER ∈ R L ER is much larger than L, with typically L ER = 500 in this study.
For a given stellar spectrum s ∈ R L ER and a substellar spectrum w ER ∈ R L ER at the same spectral resolution L ER , we obtain an (intermediate)  Algorithm 1: Pseudo-code of the selection procedure of IFS spectral priors for PACO.
Input: ExoREM set Ω ER of sub-stellar spectra. Input: ExoREM spectrum w ref ER of reference. Input: Stellar spectrum s. Input: Target number N of spectral priors. Input: Gaussian kernel g L of width L. Input: Sampling operator R L ER /L by L ER /L. Output: Library Ω IFS of spectral priors.
Step 2. Building library of spectral priors.
Ω IFS ← Ω w ∈ R L ER by dividing the two elements component-wise, i.e. w = w ER /s , ∀ ∈ 1; L ER . 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 second re-sampling the results at the targeted spectral resolution L to get a final spectral prior w.
We aim 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 Ω ER (i.e., N ER N ), and we progressively add nonredundant atoms w in Ω IFS . We first start with a single spectral prior w ref of reference in the set Ω IFS . This spectral prior is built from a sub-stellar model w ref ER 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 . We then add to the set Ω IFS the candidate spectral prior that maximized the distance, i.e., 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 repeat this procedure for various values of N and select 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, i.e. in order to keep the S/N loss significantly 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 use again eq. 1 with the number of processed pixels for IFS n pixel = 140000. 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. 8: Example of the library Ω IFS ∈ R 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.

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 ) 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 use, 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 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.

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 recentering 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 combine four terms: -A distortion error of 0.4 mas at 1 as (Maire et al. 2021), scaling linearly with the separation of a source: -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): -An error on the re-centering of the individual frames, as estimated by the re-centering procedure described in Sect. 3: -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: -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): -PACO internal error σ PA, PACO ( • ).
Those 4 terms are also quadratically combined to obtain the full error budget σ PA,tot ( • ).

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 that 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, thus 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. During the course of an observation, SPARTA collects information on the observing conditions that are then stored. As for astrometry, the PACO algorithm provides a fitting photometric error σ PACO , but additional terms are needed to estimate the global error budget: -An 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 6 . Besides, the photometry of our faint sources cannot be estimated on each frame. Hence, we use 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 SR λi (t) at the working wavelength λ 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: Then we compute the standard deviation of the SR, after removing the values for rejected frames, as follows: where SR 0 is the SR of the PSF used. In PACO, it is the average between the off-axis PSF taken before the observing coronagraphic sequence and the one after. The error associated to the extracted spectra at each wavelength λ i is therefore: 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. -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 6 https://www.eso.org/sci/facilities/develop/ao/tecno/sparta.html 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: -PACO internal error σ PACO .
Those three terms are then quadratically summed to obtain the full photometric error budget.

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 the case 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 check 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 7 . 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 process as follows: -We compute the S/N (called S/N 2 ) that the signal should have in the second epoch using the second epoch contrast 8 map. -We measure the maximum S/N in the disk area in the second epoch data (S/N max ). -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:

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  (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 . 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 pointlike 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 which are: -The optimization zone is 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; -The radial width (in radius) of the subtraction zone is set to one FWHM; -The radial to azimuthal ratio of the subtraction zone is set to 1.5; -The optimization zone area is 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 is 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/maps are 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 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 are normalized by the off-axis PSF flux. S/N maps are directly derived from the estimated flux and its associated standard-deviation (i.e., contrast at 1σ).
We present in Sect. 5.1 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. Figure 9 (respectively, Fig. 10) shows the predicted 5σ contrast as a function of the angular separation on IRDIS (respectively, IFS) obtained with the three considered algorithms for all epochs of the 24 stars considered in this study. For this comparison, we use the combined contrast of TLOCI and PCA, as well as the contrast obtained with the f lat prior with PACO (hereafter noted PACO-flat) to possibly make the most direct comparison. Although 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 √ 2, 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, i.e. 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 f lat 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.

Contrast performance
For IRDIS, at close angular separations, PACO performs better 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 3 and a factor 5, depending on the separation, compared to PCA. It is much larger compared to TLOCI by about a factor 10 for all 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.
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.

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 at quantifying this bias via numerical experiments. For that purpose, we resort to massive injections of FPSs at contrast levels predicted by PACO and we rerun 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 Sect. 5.2.1 & 5.2.2) were performed on a dataset of HD 377 (star included in this study) for both IRDIS and IFS. 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 2-D contrast maps, hence the differences with the 5σ curve: local variations of the achieved contrast are averaged azimuthally.

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), see Appendix B. For both, from the 2-D 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 9 obtained after 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) 10 . 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 results is also close to the theoretical value predicted by a Gaussian distribution, which is 1. 9 PACO does not produce residual maps since it describes the stellar component through a statistical model rather than resorting to explicit combinations and/or subtractions of images. 10 Because we use a 4σ threshold for the analysis, we did not recover injections with a S/N below 4σ. We can however still find the median S/N by computing the S/N for which half of the injected sources (the 45 higher S/N sources in that case) have a higher S/N  These results confirm that the contrast estimates produced by PACO (e.g., shown in Fig. 9) are reliable, statistically grounded, and that the detection sensitivity is im-proved compared to TLOCI and PCA, for the whole range of angular separations.

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 cases), for which the contrast injected was also computed using the 2-D contrast maps. We also consider a flat injection which 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. 6. Results on the mini-survey

Identified point sources and status
Running the PACO algorithm along with the analysis tools described in Sects. 4.2 and 4.4 over the 40 datasets of this survey allows us to identify 60 (58 IRDIS, 3 IFS) point-like sources with an S/N above the 5σ detection confidence. 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): source classified as background in the F150 or in this analysis using a proper motion analysis.  Figure 16 shows a pie-chart diagram of the sources classification. Individual information will 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 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 classify 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 classify them as candidate companions (CC) at this stage.
Finally, we find twelve false positives during 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.

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, once expressed in masses, and, whenever available in the ESO HARPS (Mayor et al. 2003) archive, radial velocity data. Indeed, the combination of direct imaging data with RV data provides a large exploration of the star's environments, from a fraction of au out to 100 au. To convert the PACO contrasts into detection limits expressed in masses, we use luminositymass 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 assume 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 analysis, thanks to im-  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.
proved 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 M Jup planets could be detected, compared to 6 M Jup 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 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.

Conclusion
We have presented here a number of upgrades made to our reduction and analysis pipeline using PACO so as to improve 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 performance are significantly improved compared to those obtained with classical algorithms such as TLOCI and PCA. Also, we have showed that PACO provides statistically meaningful S/N maps. This work paves the way to an end-to-end, homogeneous, and unsupervised massive re-reduction of archival SPHERE direct imaging data in the quest of 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.  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.