© The Authors 2025

1 Introduction and context

The first exoplanet detection around a solar-type star (51 Pegasi b, Mayor & Queloz 1995) was followed up in the same year by the first detection of a substellar companion via direct imaging (DI; Nakajima et al. 1995), opening up a new pathway toward exoplanet detection. Rapid instrumental progress and the advent of the first 8 m class telescopes allowed for the first planetary-mass object around a brown dwarf to be imaged (2MASS J1207334-393254 b; Chauvin et al. 2004) and around a star (AB Pic b; Chauvin et al. 2005), followed by the discoveries of HR 8799 bcde and β Pictoris b (Marois et al. 2008; Lagrange et al. 2009). These companions challenged the existing formation models, being much more massive than Jupiter (≥5 M_Jup) and located much further away, at separations between ten and several hundreds of astronomical units (au). Core accretion models (Pollack et al. 1996), which have been canonically applied to explain the birth of our solar system, were not able to easily explain the formation of such objects, hinting towards other complementary formation mechanisms such as gravitational instability inside the protoplanetary disk (Cameron 1978).

While radial velocity (RV) and transit techniques have been able to detect thousands of exoplanets at small separations, their ability to detect even giant planets is limited to semi-major axis less than 3 au with transit techniques, and less than 8–10 au (see e.g., Lagrange et al. 2023) with RV techniques. This has left high-contrast imaging as the only method used to reliably probe the outer region of planetary systems. Microlensing surveys have started to yield an increasing number of planetary-mass companions in the last decade (see e.g., Fig. 1 of Mroz & Poleski 2023), reaching ∼200 exoplanets mainly located in the 1−10 au range. This technique is uniquely sensitive to sub-Jovian exoplanets orbiting at separations of a few astronomical units, a population that is still out of reach of direct imaging. However, this method relies on rare events, thus requiring decade-long surveys such as the Optical Gravitational Lensing Experiment (OGLE; Udalski et al. 1992) and is still not able to provide strong statistical results for massive Jovian planets (Mroz & Poleski 2023). The future Roman space telescope will significantly increase the number of planets detected by microlensing, potentially up to a factor 10 (Penny et al. 2019; Johnson et al. 2020). As for astrometry, a small number of planetary candidates has been proposed using Gaia third data release (DR3, Gaia Collaboration 2022), although their parameters from this preliminary Gaia data may be affected by selection effects (stars properties and/or orbital configurations; see e.g. the discussions in Sozzetti et al. 2023). We note that even the preliminary results already available from DR3 may be used to boost the yield from DI and improve our knowledge of systems found by other techniques (see e.g., Winn 2022; Philipot et al. 2023a,b; Winterhalder et al. 2024).

DI faces two main challenges: (1) the very small angular separations where companions are located with regard to their host star and (2) the very high contrasts between planetary and stellar fluxes. Such limitations limit DI observations to close (typically less than 150 pc) young (less than 500 Myr) stars where the planets will be (angularly) separated from their parent star and warm enough to ensure a contrast detectable with current instruments.

Over the past two decades, high-contrast imaging surveys looking for substellar companions were conducted (for a more detailed overview, see Chauvin 2018; Currie et al. 2023a, and references therein). Starting with the first instruments using adaptive optics (AO) systems, typically targeting a few dozen stars (see e.g. Chauvin et al. 2002, 2003, 2006; Biller et al. 2007), surveys progressively increased in size and sensitivity, thanks to progress in instrumentation and image processing. Third-generation surveys, targeting ∼100 stars using coronagraphy and more advanced AO systems were able to probe the planetary regime, between 5 and 13 M_Jup, at a few dozen au (see e.g., Janson et al. 2013; Desidera et al. 2015; Chauvin et al. 2015; Vigan et al. 2017; Galicher et al. 2016). The 4th and current generation surveys consider hundreds of stars: 400 for the SpHere INfrared survey for Exoplanets (SHINE; Chauvin et al. 2017a), 600 for the Gemini Planet Imager (GPI) Exoplanets Survey (GPIES, the first half of which was published in Nielsen et al. 2019), and use dedicated planet-finder instruments equipped with extreme AO devices coupled with coronagraphy, such as the Spectro Polarimeter High-contrast Exoplanet REsearch (SPHERE; Beuzit et al. 2019) and GPI (Macintosh et al. 2006). Those surveys can efficiently probe the planetary regime beyond typically around 10 au (Vigan et al. 2021; Nielsen et al. 2019). Smaller surveys taking advantage of those facilities, more focused on very specific stellar populations, are also ongoing: for instance, the B-star Exoplanet Abundance Study (BEAST; Janson et al. 2021b) in the Scorpius-Centaurus association (hereafter Sco-Cen). With the advent of Gaia (Gaia Collaboration 2018, 2023), informed surveys are now conducted (compared to the blind surveys of the past) to try to maximize the detection efficiency of planetary companions (typically on the order of 1% for blind searches; see e.g., Currie 2020, for a survey using this approach) with three recent planetary mass companion discoveries: AF Lep b (Franson et al. 2023; Mesa et al. 2023; De Rosa et al. 2023), HIP 99770 b (Currie et al. 2023b), and Eps Ind b (Matthews et al. 2024).

The SHINE survey was conducted over 200 nights between February 2015 and September 2021, using the GTO time allocated to the SPHERE consortium. This survey was designed to find and characterize new exoplanets or brown dwarfs (BD), with the aim to explore planetary architectures, determine the frequency of giant planets beyond 10 au, and study the dependency of this frequency with the spectral type of the host star. A last aim was to study the link between disks and planets to constrain planetary systems formation processes. The first series of papers presented a sample of the first 150 stars observed during the first half of the SHINE survey (hereafter, F150) were published in 2021 with the sample characteristics presented in Desidera et al. (2021) and the performances, data processing, and source identifications in Langlois et al. (2021). From these observations, a first statistical analysis was conducted (Vigan et al. 2021), studying the frequency of giant planets, its dependency with regards to the stellar spectral types, and comparing the results to formation models. The first statistical results obtained using the F150 sample computed occurrence rates for planetary-mass companions around FGK stars of 0.7% [0.3–2.9] 68% confidence interval (CI; Vigan et al. 2021), fully compatible with the GPIES survey and its upper limit at 6.7% (95% CI interval; Nielsen et al. 2019). For higher mass stars (A and B type), SHINE and GPIES are also compatible with an occurrence rate of 8.6% [4.1–15.9] and 8.9% [5.3–13.9] (68% CI interval for both), respectively. These results still suffer from large error bars and could be improved by using additional data.

The present paper is part of a similar series with a paper describing the sample (Squicciarini et al., in prep.) and another one studying the demographics and formation models (Chauvin et al., in prep.). It adds 250 stars to the previously analyzed F150 sample to reach a total number of 400 stars, representing over 650 observations. This will be, alongside the GPIES survey (Nielsen et al. 2019) when fully released, the largest high-contrast imaging exoplanet-hunting survey ever conducted. In addition to an increased number of targets, the present analysis uses an improved post-processing algorithm compared to the one used in Langlois et al. (2021) (SPECAL Galicher et al. 2018): the PAtch-COvariance algorithm (PACO; Flasseur et al. 2018, 2020b,a)

The SHINE survey discovered two planets: HIP 65426 b (Chauvin et al. 2017b) and PDS 70 b (Keppler et al. 2018, not part of the statistical sample presented here) and one brown dwarf (HIP 64892 B; Cheetham et al. 2018) as well as four new disk detections (Lagrange et al. 2016a; Feldt et al. 2017; Sissa et al. 2018; Perrot et al. 2019), although the target selection was not optimized for it. While the SHINE survey focused on stars for which no stellar companions in the range of separation 0.1-6 arcsec at the epoch of its start (2014), a number of additional stellar companions in this separation range were discovered on SHINE data around 76 of the observed stars. Their observation is described in Bonavita et al. (2022). These stars were not considered further in our analysis.

After its completion, 1220 SHINE companion candidates remained unconfirmed mainly identified in the IRDIS field of view (FoV) around 133 stars. This large number of unconfirmed sources is due to an incomplete follow-up of detected sources due to losses of time because of a significant amount of poor observing conditions (SHINE was conducted in visitor mode) or scheduling constraints. Because both detections and non-detections are used in the statistical analysis (see, e.g., Vigan et al. 2021), their confirmation is critical to reducing biases and error bars.

A new and smaller survey called snapSHINE was proposed to the test the nature of 1010 candidates around 102 stars. The missing 31 stars harbor companion either too faint or too close to be accessible with snapSHINE. Finally, 76 out of the 102 targets were observed. The observing strategy of snapSHINE was designed to maximize the balance between observing time (as little as possible) and the number of sources to be re-observed (as much as possible). This snapshot technique has proven to be very effective to detect wide orbit exoplanets (Bohn et al. 2020a,b, 2021), while snapSHINE surveyed similar stars, for which the candidates were observable using this approach. We should, however, acknowledge that this still leaves a number of potentially good but unluckily unconfirmed candidates, especially at short separations where snapSHINE has poor sensitivity and, moreover, where we might expect a larger prevalence of companions (as seen in Figs. 9–10).

This paper is organized as follows. In Sect. 2, we present the SHINE and snapSHINE surveys with their associated observing conditions. In Sect. 3, we provide details of the data analysis and algorithms used. The performance of the survey is presented in Sect. 4, discussed in relation to the observing conditions, along with a comparison to the previous study. In Sect. 5, we present the detected point sources as well as the classification process used to identify them. Finally, in Sect. 6, we present our conclusion and perspectives for future works.

2 Surveys, related observing strategy and data quality

We present the observational strategy in Sect. 2.1 for SHINE and in Sect. 2.2 for snapSHINE. In Sect. 2.3, we discuss the overall data quality of both surveys.

2.1 SHINE observational strategy

All SHINE observations were taken using standardized observational templates using SPHERE. Two modes were used, either IRDIFS or IRDIFS-EXT, both of which allow simultaneous observations by the Infra-Red Dual-beam Imager and Spectrograph (IRDIS; Dohlen et al. 2008) and the infrared Integral Field Spectrograph (IFS; Claudi et al. 2008; Mesa et al. 2015) in complementary wavelength bands. IRDIS was used in the dual-band imaging mode (DBI; Vigan et al. 2010), either with the DB-H23 filter pairs in the case of the IRDIFS mode with IFS using the YJ bands, or with the DB-K12 filter pairs and YJH bands for IFS in the IRDIFS-EXT configuration. The latter configuration was mostly used for observations towards Sco-Cen (a 10–20Myr association located at 100−150 pc) stars to maximize the detection capabilities of red young late L-type dusty planets (e.g., HD 965086 b; Chauvin et al. 2018; Desgrange et al. 2022).

To efficiently schedule all of the SHINE observations and to maximize the amount of FoV rotation, critical for postprocessing techniques relying on angular (and spectral) differential imaging (A(S)DI; Marois et al. 2006), the GTO team developed a dedicated optimal scheduler: SPOT. This tool was described in Lagrange et al. (2016b) and in Langlois et al. (2021). It was used to optimize the observations for the whole survey while verifying specific constraints such as a variation of the parallactic angle of at least 20^∘ and as close as possible to 30^∘ as can be seen in Fig. 2. We note that some datasets with small rotation angles are still used: for instance, two datasets of GJ 504 have a measurement of 5^∘ in terms of the parallactic rotation, but they still allow for a robust detection of the planet/brown dwarf (depending on the adopted age) GJ 504 b (see e.g., Kuzuhara et al. 2013; Bonnefoy et al. 2018) because of its large angular separation (around 2.5 as) relative to its host star.

As previously mentioned, all SHINE observation sequences were standardized, consisting of an observation of the noncoronagraphic off-axis point-spread function (PSF) of the targeted star for flux calibration purposes. These were followed by a coronagraphic frame with four replicas of the star created by the deformable mirror (Makidon et al. 2005) and used to check the centering of the star behind the coronagraphic mask (so-called waffle-frame produced by dedicated commands on the deformable mirror to create the four replicas). Once the star is centered, the science coronagraphic sequence is recorded. Finally, another waffle frame and another off-axis PSF are recorded, to estimate the astrometric and photometric changes between the beginning and the end of the observations, followed by a sky frame for calibration purposes. Some specific observations were taken with the waffle mode during the whole science sequence, mainly in the case of orbital monitoring of known companions. Indeed, such observations allow a frame-by-frame re-centering, thereby greatly reducing the astrometric uncertainties due to any potential jitter of the telescope during the science acquisition. The standardization of most of SHINE observations also allows for a very consistent total integration time as can be seen in Fig. 3, with most of our observations achieved between 3600 s (the theoretical minimum observing time for a SHINE observation) and 5700 s of total science integration time (7200 s was the maximum time allowed).

2.2 snapSHINE observational strategy

This program was optimized to detect sources at a separation greater than 1 arcsec (corresponding to the 50−300 au range) with shallower observations compared to SHINE, with a typical sensitivity greater than 2 M_jup at most. The observations are around 28 minutes long (including telescope overheads), with a low FoV rotation, typically smaller than 4 degrees. Close or very faint sources (around 15% of the total unconfirmed sources) were then out of reach of the survey. This trade-off in sensitivity was made to maximize the number of observations, while still being sensitive to wide-orbit exoplanets such as those found by Bohn et al. (2020b); Janson et al. (2021a); Squicciarini et al. (2021).

The observations and calibration templates are similar to the SHINE survey with one minor change: the waffle mode is enabled for the whole coronagraphic sequence, to get the best astrometric precision and unambiguously identify background sources via the common proper motion test (see Sect. 5.2.1). As most unconfirmed candidates were detected using the IRDIFS setup, the latest is also used for the snapSHINE survey, with IRDIS using the BB_H filter. The change from DB_H23 to BB_H filter for IRDIS was done to maximize the number of photons received (by increasing the bandwidth) given the relatively small exposure time of the snapSHINE observations.

In the present analysis, we considered all observations performed until the end of P111, even if the program has been extended to P113. A full analysis of the survey will be presented in a dedicated paper after its completion.

Fig. 1 Overall data quality for SHINE and snapSHINE. (a) Histogram of the star magnitudes (G, H, K) for the sample considered in this study. (b) Scatter plot showing the coherence time (τ0) against the seeing with the Strehl ratio color-coded. Data points without SPARTA Strehl information are colored in white. The grey box delimits the ESO 50% best observing conditions (seeing <1 as and τ0 > 3.2 ms) whereas the yellow box delimits the 30% best (seeing <0.8 as and τ0 > 4.2 ms). (c) Scatter plot showing the 1σ raw contrast at 0.5 arcsec (directly measured on the coronagraphic frame, before post-processing) with respect to the H-band (main IRDIS band used in SHINE) star apparent magnitude. The best contrast performances are achieved for the brightest stars with a high Strehl ratio. (d) Scatter plot displaying the Strehl ratio as a function of the star’s apparent magnitude in the G band (AO spectral channel). As expected, the performances of the AO system deteriorate for fainter stars.

Fig. 2 Histogram of the total amount of parallactic angle rotation during the coronagraphic sequences for the SHINE observations considered in this survey.

Fig. 3 Histogram of the total coronagraphic integration time during the sequences for the SHINE observations considered in this survey.

2.3 Overall data quality

The whole SHINE survey has been conducted in visitor mode, hence, the quality of observations is highly variable. It is possible to retrieve information on these conditions either via the ESO MASS¹-DIMM² or using the parameters related to image quality obtained by the real-time computer SPARTA (Suárez Valles et al. 2012) of the SAXO extreme AO of SPHERE (Fusco et al. 2016; Beuzit et al. 2019). The seeing and the coherence time are retrieved using the ESO MASS-DIMM measurements. We used the SPARTA Strehl ratio; however, it is not available for all observations: data are missing for the 115 observations taken between 2015-02-05 and 2015-05-19 and between 2016-06-02 and 2016-12-16. We note that the SPARTA Strehl can be slightly different from the ESO MASS-DIMM Strehl (see, e.g., the Fig. 8 of Milli et al. 2017, and discussion therein). We can get information on the observing conditions through the prism of the ESO turbulence category in a similar fashion as in Courtney-Barrer et al. (2023), shown in Fig. 1b. SHINE is conducted in visitor mode and so, we can assume that the observing conditions during the 200 nights were random. We should then expect that around 50% of the observations were in the 50% turbulence category (i.e., seeing <1 as and τ₀ > 3.2 ms). We indeed found that 48.35% (311) of the observations are within this parameter space. 66.14% (424) of the SHINE observations are within the 70% (seeing <1.15 as and τ₀ > 2.2 ms) turbulence category and 34.1% (221) belong to the 30% (seeing <0.8 as and τ₀ > 4.2 ms) category, showing that the observing quality of a given SHINE observation behaves as expected from a collection of random observing nights. This impacts our ability to confirm faint candidate sources with a second epoch if this was obtained in poorer conditions; or, on the contrary, we can achieve a deeper second epoch, unveiling new candidates undetected in the first epoch if the conditions were better.

Figures 1c, d also illustrates the degradation of the performances with increasing star magnitude. This is quantitatively shown in Fig. 1c, which provides the 1σ raw contrast, computed on each coronagraphic frame, and then the median combined over the whole cube.

For snapSHINE, the survey was conducted in service mode, with a constraint of compliance with the 70% turbulent category (seeing <1.15 as and τ₀ > 2.2 ms). About 98.7% (75) of the observations were within the specified constraint. The observing conditions and data quality are summarized in Fig. 1, alongside the same details on SHINE.

3 Data reduction and analysis

3.1 Calibrations

The standard calibration strategy described in Langlois et al. (2021) was adopted throughout the entire survey. As for astrometry, we used the long-term calibration computed by Maire et al. (2021) for all detected point sources.

3.2 Pre-processing

The pre-processing is identical to the one used for the F150, see Langlois et al. (2021) for a detailed explanation of each step performed. In a nutshell, the first steps of the pre-reduction, consisting of dark, background, flat and bad pixel corrections, were performed at the High-Contrast Data Center (HC-DC³, formerly known as SPHERE Data Center, Delorme et al. 2017) by using the SPHERE data reduction and handling pipeline (DRH; Pavlov et al. 2008) provided by ESO.

For IFS, a few additional steps are used to improve the wavelength calibration, implement a correction for cross-talk during the spectral extraction, and improve the bad pixels correction (Mesa et al. 2015; Delorme et al. 2017). The only steps that differ from the F150 pre-reduction is the recentering of the individual frames. This was already presented in Chomez et al. (2023a); Dallant et al. (2023). In brief, we fit a 2D Gaussian on each waffle pattern to estimate their subpixel location and infer the position of the star with subpixel precision, allowing thus precise recentering. It has been shown in Dallant et al. (2023) that this re-centering increases the signal-to-noise (S/N) of detected sources by a few percent compared to the previous pipeline.

3.3 ASDI post-processing

One critical part of processing such a survey is the homogeneity of the post-processing reduction and the ease of the analysis that comes afterward. During the first part of the SHINE F150 analysis (Langlois et al. 2021), classical algorithms were used such as TLOCI (Marois et al. 2014) and PCA (Soummer et al. 2012; Mesa et al. 2015), embedded into the SPECAL software (Galicher et al. 2018). Those algorithms have the same three main drawbacks: (i) they provide non-statistically reliable detection limits (especially close to the star), thus leading (in practice) to many more false alarms than theoretically expected at the prescribed detection confidence, and requiring (ii) a manual visual and time consuming source detection by an astronomer; additionally, (iii) differential imaging techniques heavily impact the shape of the off-axis PSF due to self-subtraction, introducing biases in the retrieve astro-photometry and requiring additional reduction steps (e.g., injection of fake planets) to calibrate those biases.

To re-reduce all SHINE data (see Sect. 3.4 for the algorithm use for snapSHINE), we chose to use the PAtch COvariance (PACO; Flasseur et al. 2018) algorithm as our baseline for both IFS and IRDIS in an automated pipeline detailed in Chomez et al. (2023a). More specifically, we used robust PACO ASDI algorithm (Flasseur et al. 2020a,b). The underlying models capture the statistics of the noise at a local scale of small patches through mixtures of scaled multi-variate Gaussians. The algorithm allows for a better estimation of the spatial and spectral correlation of the noise than classical methods. Chomez et al. (2023a) demonstrated on a sub-sample of the F150 that this approach improves the IRDIS contrast by one to two magnitudes (x3-x6), depending on the angular separation, compared to more classical ADI approaches such as TLOCI and PCA. The gain for IFS is more modest because IFS data are more photonlimited than IRDIS. In both cases, the detection limits provided by PACO are statistically reliable, thus allowing for a realistic false alarm rate.

The ASDI option of PACO is used to jointly process the temporal and spectral dimensions, producing spectrally combined S/N maps, and using spectral weights to maximize the detection of candidate sources having a similar spectral signature. These weights, with as many components as wavelength (L = 2 for IRDIS and L = 39 for IFS) are referred to as spectral priors. The selection, optimization process, and associated detection performances of the spectral priors are fully explained in Chomez et al. (2023a).

3.4 No-ADI post-processing

Because of the small number of frames and the small FoV rotation, we were not able to process the snapSHINE observations with PACO. Instead, we used the NoADI routine (a simple derotated median stack of all images) embedded in SPECAL (Galicher et al. 2018) for the IRDIS processing, which is well adapted given the small △PA variations, and since the sources of interest are located in the background limited region. Additionally, we processed the IFS PSFs using the algorithm described in Bonavita et al. (2022). In a nutshell, this algorithm uses the off-axis PSFs taken before and after the coronagraphic sequence and computes a subtraction between the two after spectral-collapsing, normalization, and re-aligning of PSFs. In case of very bright sources (mainly binary stars) in the IRDIS FoV that can saturate the coronagraphic frames, they could not be characterized. Thus, we also used the no-ADI algorithm on the off-axis PSFs, as they are often taken with a neutral density to avoid saturation. Four binary stars were characterized using this method.

4 Detection sensitivity and comparison with the F150 analysis

In this section, we present the performances achieved for SHINE and snapSHINE. As a first step, we present the contrast curves for both surveys in Sect. 4.1 as well as a discussion on those limits. In the second part, in Sect. 4.2 we explain how we translated the contrast curves into detection limits expressed in Jupiter masses for the same stars considered in the F150 sample (Desidera et al. 2021; Vigan et al. 2021) to quantify the improvement using PACO ASDI compared to the previous analysis (Langlois et al. 2021).

4.1 Contrast curves and performance comparison regarding observing conditions

The achieved performance details, in terms of contrast, are presented in Fig. 4 for both the IRDIFS mode and the IRDIFS-EXT mode.

Compared to the F150 analysis (Langlois et al. 2021), the IRDIS performances with the IRDIFS setup are significantly improved, reaching a median contrast of 10⁻⁶ at 0.8 arcsec and 2 × 10⁻⁷ beyond 2.5 as. The gain for the IFS data is not as high, only slightly improving the detection limits. However, the SPECAL detection curves are not reliable with regards to the false alarm probability at 5 σ, when PACO is, especially at short separation. Even if the contrast improvement is moderate, the statistical reliability of the contrast curves derived by PACO, crucial for the statistical analysis, is much stronger. For the IRDIFS-EXT setup, the achieved contrast at short separations (below 1 as) is comparable to the one achieved in the IRDIFS setup. At larger separations (>1 as), the contrast in the K band is limited by the strong thermal (background) noise impacting IRDIS images, thus reducing the achieved contrast compared to the H band, where this type of noise is significantly lower.

Following Courtney-Barrer et al. (2023), we plotted the raw and post-processing contrasts using bins of stellar G magnitude and the ESO SPHERE turbulence category classification. We considered three categories: G < 6, 6 < G < 9, and G > 9. All datasets impacted by the low-wind effect (windspeed <3 m s⁻¹; Milli et al. 2018) have been removed from this study. Additionally, only datasets with a FoV rotation between 20^∘ and 40^∘, and a total integration time between 2700 and 5400 s were considered. This parameter space is representative of the majority of our observations (see Figs. 2 and 3), although it is still tight enough to mitigate the impact of those variables on the plots, allowing us to only consider the star brightness and the observing conditions⁴. Additionally, we estimated the gain provided by PACO ASDI by computing the ratio between the median raw curves and the median post-processed curves for each star brightness and observing category. As expected, the raw contrast performance (i.e., directly measured in the raw coronagraphic frames; in our case, median-averaged spectrally and temporally over the whole 4D cube) decreases with the increasing faintness of the star, as can be seen in Fig. 5. While on IRDIS, PACO improves the contrast by a factor around 100 for separations below 1 arcsec for all stars’ brightnesses, the gain is much more dependent on the star brightness for IFS: for bright stars, the gain can go up a few hundreds, achieving better performances at small angular separations than IRDIS; but for faint stars, the gain is comparable or worse than for IRDIS. As shown in Fig. 5, for the bright and intermediate stars, we clearly see the improvement in contrast due to the AO correction between 200 and 500 mas (the dip in the raw contrast curves). In those cases, the data are still speckles-limited. However, this dip disappears for faint stars, becoming a plateau. This indicates a switch in the noise regime, from a contrast-limited one to a background-limited one. Given the already worse raw contrast in that case, it might explain why IFS can underperform compared to what is traditionally expected and explain the broader dispersion compared to IRDIS in the post-processed contrast seen in Fig. 4⁵. Finally, we note the excellent performance of SPHERE coupled with PACO ASDI when observing bright stars under the best observing conditions; for example, contrasts better than 10⁻⁶ at 0.5 arcsec can be achieved.

Finally, the contrast curves derived with SPECAL for snapSHINE are presented alongside the ones obtained with PACO ASDI in IRDIFS mode in Fig. 4 for comparison. As expected, snapSHINE is much less sensitive below 2 arcsec but reaches the expected contrast in the background limited area, being only less sensitive than SHINE by about two magnitudes. Figure 6 shows as an example a SHINE and a snapSHINE image of the same target. While the three brightest sources in the FoV are recovered in the snapSHINE image, the faintest one is not.

Fig. 4 Contrast curves at 5σ obtained on the whole survey for both instruments, IRDIS in green and IFS in red. (a) IRDIFS setup. (b) IRDIFS-EXT setup. The contrast curves of the snapSHINE survey are represented in blue. The grey overlaid rectangles at small angular separations show the average inner working angle of the coronagraphic masks.

Fig. 5 Contrast comparison between the raw contrast (median average spectrally and spatially, dashed lines) and the PACO contrast, both at 5σ. Three bins of star magnitudes and observing conditions are considered. Individual observations are represented by the thin lines and the median for each observing condition by the thick lines. The gain gives the contrast improvement between the median raw curves and the median PACO curves, provided for different observing conditions and different star brightness. The average size of the coronagraph’s inner working angle is around 100 mas.

Fig. 6 Comparison between a SHINE and a snapSHINE observation of the same star. (a) S/N map of the star TYC_8634_1393_1 with the SHINE observation using PACO (τ = 5σ). (b) Residual map of the same star with the snapSHINE observation using NoADI. Four detections are present in the PACO S/N map, highlighted by the dashed squares. Only the three red ones are redetected in the snapSHINE observation (background sources). The other one in yellow is too faint to be detected by snapSHINE. A 15 × 15 pixel 2 D median filter was applied on the snapSHINE residual map and the detection was performed manually.

4.2 Comparison with the F150 early analysis

To quantify the gain using PACO, we compared the detection limits obtained in this work to the one published in Vigan et al. (2021), only using the stars included in the F150.

As in Vigan et al. (2021), we used the 1D contrast curves and convert them into detection limits in Jupiter masses using a recent update of MADYS⁶ (v1.2, Squicciarini & Bonavita 2022), with the ages extracted from Desidera et al. (2021) and the parallaxes and fluxes from SIMBAD. For the conversion, we used the COND2003 model (Baraffe et al. 2003). To convert the mass curves into detection limits, we used EXO-DMC⁷ (Bonavita 2020). The grid parameters are designed to be similar to the ones used in Vigan et al. (2021), with the semi-major axis ranging from 0.1 to 1000 au (500 points, log uniform sampling) and the masses from 0.1 to 100 M_Jup (200 points, log uniform sampling). A thousand draws were made for every cell in the grid with the inclination taken from a uniform distribution without any constraints and the eccentricity from a Gaussian law following $\mathcal{N}$(μ = 0, σ = 0.3). Figure 7 clearly shows the gain in sensitivity obtained using PACO compared to the post-processing algorithms used in Langlois et al. (2021). The sensitivity is improved at all separations: for instance, the parameter range within which we are able to identify companions around 93% of the sample (140 stars) region is greatly extended, reaching 6 M_Jup at 100 au. At 10 au, while the previous analysis was able to reach 4 M_Jup for 20% (30 stars) of the sample when using PACO ASDI, we achieve a detection limit of 2 M_Jup. We are even able to go beyond the 1 M_Jup sensitivity for 20% of the F150 sample at 70 au, performances that the SPECAL analysis could barely reach at 200 au. We then computed the subtraction between the PACO analysis and the analysis presented in Vigan et al. (2021). This map is computed by subtracting, for every cell in the grid, the probability obtained using PACO ${\mathcal{P}}_{\text{PACO}}$(sma, M_Jup) by the one obtained using SPECAL ${\mathcal{P}}_{\text{SPECAL}}$(sma, M_Jup) and provide a clear view of the area of the parameters’ space where we improve the most compared to the previous analysis.

5 Candidate companion classification and analysis

In this section, we summarize the confirmed planets or brown dwarf companions detected during the second part of the survey in Sect. 5.1. Sub-stellar companions already presented in Langlois et al. (2021) and stellar binaries studied in Bonavita et al. (2022) will not be discussed here but their astro-photometry extracted with PACO will still be included in the VizieR table linked to this article, as the whole survey has been re-processed. The source classification scheme used for the survey is presented in Sect. 5.2. Online CDS tables with the astro-photometry of all IRDIS and IFS sources as well as all IFS spectrum for sub-stellar objects are linked to this article.

5.1 Confirmed exoplanets and brown dwarfs

One additional substellar companion has been confirmed during the second part of the SHINE survey: HIP 74865 B. Already detected in 2015 (Hinkley et al. 2015) using Sparse Aperture Masking (SAM) on VLT/NACO (Rousset et al. 2003), it is a brown dwarf of about 30−60 M_Jup orbiting a 15 Myr Sco-Cen star at about 20 au of its host star. The first image of the companion has been taken in the framework of the SHINE survey (see Fig. 8). Following this confirmation, two additional SPHERE epochs were recorded using the SAM mode of SPHERE (Cheetham et al. 2016). Those three observations and an in-depth analysis of this brown dwarf will be presented in a forthcoming paper (Cantalloube et al., in prep.). The limited number of additional confirmed companions with respect to F150 should not be considered surprising because in the course of the SHINE survey, higher priority was given to observations of the best targets, where detections were more probable.

Fig. 7 Detection performance comparison between this work and the detection limits computed in Vigan et al. (2021). Left: comparison between the sensitivity achieved by Vigan et al. (2021) using classical post-processing algorithms (Langlois et al. 2021) on the F150 (dashed lines) compared to the sensitivity achieved using PACO on the same sample of stars (plain lines). The lower limit seen around 100 au at 0.5 MJup is due to the COND model, which does not provide inputs for lower masses. Right: probability increase (${\mathcal{P}}_{\text{PACO}}$(sma, MJup) – ${\mathcal{P}}_{\text{SPECAL}}$(sma, MJup)) between the PACO analysis and the one presented in Vigan et al. (2021).

Fig. 8 S/N maps for (a) IRDIS and (b) IFS of HIP 74865 B. The coronagraph is represented by the yellow circle. The brown dwarf is clearly visible southeast of it. Maps are cropped at a 500 mas radius. The color bar represents the S/N.

5.2 Classification of the detected point sources

In order to ensure the most reliable and complete statistical analysis that will come after this work, we need to classify all sources detected during the survey. This classification must be performed as thoroughly as possible by considering each and every source as a potential planetary candidate with the objective of leaving as few as possible candidates with an ambiguous status. Those ambiguous candidates are the ones impacting the statistical analysis as their nature (i.e., background contaminants or bound sources) cannot be accurately assessed. We want to stress here that there are two ways of counting the detections: one where we sum all the detected point sources and one where we consider only the physical sources. For instance, there are six HR8799 datasets in the present analysis, where all four planets (on IRDIS) are detected on every epoch. The number of physical sources is four (one per planet), but the number of detection is 24. Combining over 700 observations, 7223 (5858 SHINE + 1365 snapSHINE) point sources were detected. In the following, both numbers will be given if different. The lack of a second observation for every star as well as heterogeneous observing conditions calls for a complex identification scheme. Its various steps are described in the following subsections. As it requires a photometric and astrometric analysis of each detected point source (as done for actual exoplanets), this process is particularly time-consuming.

5.2.1 Common proper motion test

The most robust method to assess the physical link between a candidate companion and a star is the common proper motion test because it does not make any hypothesis on its spectrum. It compares the motion of a detected point source between two observations and assesses if it is consistent with a co-moving (i.e., gravitationally bounded) object or a stationary background star with a null apparent on-sky motion (or with the motion of background objects as determined by e.g. Gaia astrometry). An example of both cases is presented in Appendix A. If the time between the two considered epochs is sufficiently large (typically, at least a year for the objects in the SHINE survey), this method can robustly separate the background objects from the bounded ones. By coupling SHINE and snapSHINE, we were able to identify 1180 physical (corresponding to 4217 detections) background sources and 26 physical bound sources (corresponding to 90 detections).

5.2.2 Direct Gaia detection

As a second step, we used a tool to cross-match detected point sources in SPHERE images with the Gaia DR3 detected point sources catalog (Gaia Collaboration 2023). The tool works as follows: it queries the Gaia DR3 archive around the target and retrieves all other Gaia detections in a given radius around the target (fixed at 8 arcsec for this analysis, slightly larger than the IRDIS FoV to avoid missing candidates). Then, the RA and DEC of all retrieved sources are propagated from the Gaia observations to the SPHERE epoch and converted in separation and parallactic angle to be easily identified in the S/N map. The tool also provides various information such as the proper motion, parallax, and fluxes of identified sources if available in the Gaia archive. An example is shown in Fig. B.1. This script allows for the identification of background sources or binary systems, even if only one observation is available. However, due to the intrinsic limitations of Gaia, we cannot use this method for candidates fainter than a G magnitude of around 21 or too close to the star, typically less than 2 arcsec in projected separation (see Squicciarini et al. 2022; Viswanath et al. 2023; Chomez et al. 2023b, for examples of non-applicable cases). This tool allowed the identification of 34 additional background sources as well as two stellar binaries, one of which was already identified in Bonavita et al. (2022).

5.2.3 Color magnitude diagram

For the remaining detected point sources, due to the already mentioned limitations of the follow-up observations, we have to use color-magnitude diagrams (CMDs) to classify them. This method was used for 1926⁸ sources. The main goal of such a diagram is to try to disentangle potentially interesting sources from background contaminants to efficiently plan for follow-up observations. The design and usage of such a diagram have been described in detail in Bonnefoy et al. (2018) and its initial usage by the GTO team in Langlois et al. (2021). As mentioned in Langlois et al. (2021), due to the IRDIS FoV, a large number of confirmed background sources (e.g., confirmed by using the common proper motion test) have been detected. Those background sources, mostly expected to be field M and K stars (Parravano et al. 2011), fall within a particular area of the CMD. Based on the large size of this sample, we chose to re-define the exclusion area presented in Langlois et al. (2021), taking advantage of the sheer number of background contaminants that we characterized in detail. This empirical distribution of bona-fide astrophysical contaminants is by construction representative of the true dispersion of contaminants in the CMD, accounting for the noise measurement but also for the natural color dispersion of background sources in SPHERE data. The statistics of the color C_color(H2 − H3 or K1 − K2) of confirmed background sources was computed (mean C̄_color and standard deviation σ_color) using bins of one magnitude to ensure enough data points. Figures C.1-C.2 display the corresponding histograms with their Gaussian fit.

The exclusion area is then defined by the area between C̄_color ± n σ_color for each bin of absolute magnitude. Here, we chose to set n = 5. Those values are then linearly interpolated to obtain the final exclusion area. All sources whose position on the CMD falls in this exclusion area are then considered as background sources. This method is applied at magnitudes fainter than the L−T transition where the expected track of T and Y objects differs from M-L type objects. However, at magnitudes brighter than the L-T transition (seen at M_H2 = 14 in Fig. 9), background contaminants and sub-stellar companions overlap on the CMD. This is why we only used the background exclusion area below the L−T transition with an additional buffer of one magnitude below the absolute magnitude of HR 8799 b to avoid miss-classification of under luminous planets. The exclusion area thus defined can be seen in the CMDs inFigs. 9–10.

Additionally to the CMD, we associate, for every remaining detection where the exclusion area can be used, a Z-score, Z aimed at identifying the most promising sources. We define this score Z score by: $$Z=\frac{|C_{\mathrm{d}\mathrm{e}\mathrm{t}}-{\overline{C}}_{\mathrm{b}\mathrm{c}\mathrm{k}\mathrm{g}}|}{{\sigma }_{\mathrm{d}\mathrm{e}\mathrm{t}}},$$(1)

with σ_det as the error on the color of the source.

Although this technique is efficient with the H23 CMD, due to the increased thermal background noise in K1 and K2, the use of the K12 CMD is much more difficult. The error bars on each data point are typically larger than in H23, thus increasing the dispersion of the color K1 − K2. As can be seen in Fig. 10, adopting the same method as in Fig. 9 leads to a very wide background exclusion area, not suitable to disentangle between interesting sources and field contaminants. We choose then to mark the candidates between 3σ and 5σ as ambiguous because while within the expected dispersion range of background sources, they are also compatible with the CMD track of young T-type exoplanets. Additionally, all candidates sharing the same color and magnitude as M and L-type objects will also be marked as ambiguous as the CMD cannot provide useful information to identify the most promising ones. We also want to stress that a detection indicated as promising based on colors should not be considered as a new exoplanet detection, but rather as the most promising candidates pending confirmation following second epoch observations.

Fig. 9 H2–H3 color-magnitude diagram with every non confirmed bound source detected during the SHINE survey. Error bars for background sources are removed for clarity purposes. The background exclusion area is indicated by the grey area and corresponds to a 5σ exclusion area.

Fig. 10 Same as Fig. 9, but for the K12 IRDIS filter and a 3σ exclusion area. The 5σ area is represented by the dotted lines. Sources between 3 and 5σ are labeled as AMBIGUOUS (see explanation in text Sect. 5.2.3).

5.2.4 False positive identification

Given the unprecedented scale of this study and the number of observations, false positives are to be expected. Identifying false positives is a non-trivial task: it requires (at least) two observations and the confirming observation (i.e., the epoch where we are trying to re-detect a previously detected source) must be of higher quality to avoid any statistical noise fluctuation that would prevent the re-detection.

Using this strategy, we were able to identify 191 false positives (100 on IRDIS, 91 on IFS). Their on-sky astrometry can be seen in Figs. D.1 and D.2, representing a polar plot (sep, PA) of all confirmed false positives for IRDIS and IFS, as well as their radial distribution. Their distributions present a slight increase at short separation in the speckle dominated area compared to a simple unbroken function of r². This excess is partly due to speckles causing extra false positives and partly because the quality of the noise estimation at these shorter separations is more affected when the FoV rotation is low, given the lack of data diversity. This latter effect is a statistical counterpart of the self-subtraction that affects sources in ADI-based algorithms.

Chomez et al. (2023a) estimated the false positive rate for both IRDIS and IFS by using observations with an inverted PA rotation to destroy any physical signal, leaving only false positives. Using these results, we can estimate the expected number of false positives: 260 for IRDIS (0.4 per observation, ∼650 observations) and 162 for IFS (0.25 per observation, ∼650 observations). As demonstrated in Chomez et al. (2023a), we can properly assess false positive rates under controlled conditions and the numbers then derived are representative of a standard SHINE observation. However, because of the variable observing conditions (see Fig. 1) and given the presence of a variety of astrophysical sources in our data, we acknowledge that the false positive rate will likely be higher, although this number is reduced by automated and human flagging of spurious detections. Given the large number of sources without a second epoch in the BCKG_CMD and PROMISING categories, it is likely that some of them, especially at low S/N, are false positives. A more detailed discussion about false positives is available in Appendix D.

Fig. 11 Pie chart summarizing the number of physical sources in each category for the SHINE F400 survey. CC stand for Confirmed Companion, BCKG_GAIA for sources directly detected by GAIA and confirmed as background stars, BCKG_PM for background stars confirmed via proper motion analysis and BCKG_CMD, sources classified as background using a color magnitude diagram.

5.2.5 Summary of identifications

The final classification is as follows: “CC” for confirmed companion (stars, brown dwarfs, and planets) via proper motion, “BCKG_PM” for background sources confirmed via proper motion, “BCKG_CMD” for background sources identified via the CMD’s exclusion area, “BCKG_GAIA” for background sources identified with Gaia, “PROMISING” for promising candidates (e.g., candidates whose positions on the CMD is deemed interesting but no confirmation via proper motion can be made, either because no second epoch are available or the archival data are not of good enough quality) and finally “AMBIGUOUS” for ambiguous cases, where the CMD cannot distinguish between background and bound sources or the common proper motion test is non-conclusive or no reliable astro-photometry can be extracted. The pie chart in Fig. 11 summarizes the classification of the physical sources detected by SHINE according to the classes described above.

This analysis leaves us with 434 sources without a conclusive status (AMBIGUOUS + PROMISING). Most of those sources are likely to be background sources, but their exact nature cannot be confirm at this stage; hence, it is possible that several of these detections are indeed real bound companions. Because the detection of any additional planet or brown dwarfs will have a strong impact on the derived statistics, mitigation strategies need to be implemented. While a detailed explanation of such strategies is beyond the scope of this paper and will be fully explained in Paper VI, an example of such process can be found in Delorme et al. (2024). In that work, the contrast map has been thresholded due to the presence of an unclassified source.

Figure 12 show an histogram of the cumulative normalized distribution of the separation of all detected sources, divided according to their classification classes. Interestingly, the cumulative distribution of the sources classified as background with the proper motion (astrometric determination) are a nearly perfect match to that of the background classified using the CMDs (photometric determination). To compare the two distributions, we used the Kolmogorov-Smirnov test (see e.g., Massey 1951) with the null hypothesis of them being drawn from the same distribution. Using an α level of 0.05 (corresponding to 95% confidence), we find that the two samples are drawn from the same distribution (p-value of 0.09). We can then be confident in our classification scheme, using the CMD for those two categories. Additionally, they match the (r/r_max)², which is the expected behavior: the number of sources visible in a FoV should increase quadratically with its radius.

As for the comparison between the AMBIGUOUS and BCKG_PM distribution, the p-value is 0.006, suggesting that the two underlying distribution are different. We note that it is possible that we are still suffering from two biases: (i) with the use of snapSHINE and due to the variable SHINE observing conditions, we are sometimes unable to confirm sources close to their stars and (ii) most of those ambiguous cases are around few stars. Given the limitation of SPHERE (both in term of sensitivity and inner working angle), it is unlikely that more than one or two sources per star are real bound companions (this would even stand as an overly optimistic bound) and most of them are likely background stars. The number of individual stars harboring an ambiguous candidate is 69, which is six times lower than the number of ambiguous candidates. However, this excess (especially at short separation) can be seen as a hint toward the presence of real physical bound companion in this class. We then carefully checked the literature for existing archival observations taken outside the framework of SHINE (mainly with NACO, GPI, and SPHERE) to vet ambiguous candidates. This excess is quite similar to the one shown in Fig. D.1 and can be interpreted as an excess due to false positives. However, the S/N of the AMBIGUOUS sources, given their relatively high absolute magnitude, is typically larger than the one of the false positives. While it is likely that some false positives are hidden in this class, is it unlikely that they are the unique cause of this excess.

As for the two last classes (PROMISING and BCKG_GAIA), the number of sources within is too small to compute any reliable statistics. We can still see the limitation in the separation of the sources detectable by Gaia, as more than 90% of them are located at a separation greater than 3.5 arcsec.

Fig. 12 Cumulative normalized distribution of the separation of the sources according to their classes. Classes colors are identical to those in Figs. 9, 10, and 11.

6 Conclusion

In this work, we analyze in a homogeneous manner, over 700 datasets using advanced post-processing techniques. This survey, alongside GPIES (Nielsen et al. 2019) and its recent reprocessing using PACO (Squicciarini et al. 2025), is the largest survey ever conducted in the search for young giant planets in the solar neighborhood. We derived improved and statistically reliable contrast curves, critical to the upcoming study on demographics and formation mechanisms.

Combining SHINE and snapSHINE, over 7000 detections were reported, each needing a precise and thorough analysis to leave as few as possible sources in an ambiguous state. This tedious, but essential, analysis is mandatory to derive meaningful statistical constraints on the population and formation models.

Such blind surveys, designed before the advent of Gaiainformed surveys, are still an essential component of planetary studies as they allow for the unbiased study of large statistical samples of stars. Although no new confirmed exoplanets have been found in this re-analysis and in the second half of the survey, due to the increased sensitivity, 24 additional promising sources have been unveiled. Those sources are either too far away from their host stars or (if bound) too light to cause any significant Gaia signal, needing further observations to confirm their nature.

Instrumental improvements such as SAXO+ (an upgrade of the AO system of SPHERE; Boccaletti et al. 2020) and GPI 2.0 (Chilcote et al. 2022) or improved observing techniques such as “dark hole” on SPHERE (Potier et al. 2020, 2022) will allow us to increase the sensitivity at small angular separations, unveiling potential new sub-stellar companions. Coupling such surveys as SHINE with other indirect techniques such as radial velocity and absolute astrometry (mainly with Gaia DR4, scheduled in 2026) will be precious to study the radial distribution of giant exoplanets from a fraction to hundreds of astronomical units.

Data availability

Tables with the astro-photometry of all detected sources by IRDIS and IFS, as well as the IFS spectrum of all sub-stellar companions, are available at the CDS via anonymous ftp to cdsarc.cds.unistra.fr (130.79.128.5) or via https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/697/A99.

Acknowledgements

This project has 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 High Contrast Data Centre, jointly operated by OSUG/IPAG (Grenoble), PYTHEAS/LAM/CeSAM (Marseille), OCA/Lagrange (Nice), Observatoire de Paris/LESIA (Paris), and Observatoire de Lyon/CRAL, and is supported by a grant from Labex OSUG@2020 (Investissements d’avenir – ANR10 LABX56). A.Z. acknowledges support from the ANID – Millennium Science Initiative Program – Center Code NCN2021_080 and NCN2024_001. T.B. acknowledges financial support from the FONDECYT postdoctorado project number 3230470. We are grateful to the anonymous referee for the insightful comments provided during the peer-review, which largely contributed to raising the quality of the manuscript.

Appendix A Example of proper motion tests

Example of two source identifications using the common proper motion test:

• HD 95086 b (see e.g. Chauvin et al. 2018; Desgrange et al. 2022) in Fig. A.1

• A background source in Fig. A.2.

This method uses the proper motion of the targeted star to disentangle co-moving sources (i.e., gravitationally bounded) from background sources, usually remote stars with little to no on-sky proper motion.

Fig. A.1 Common proper motion plot for the planet HD 95086 b. The position of the planet at the first epoch is represented by the black cross, the expected position at the second epoch, if background, by the red hollow cross, and the measured position at the second epoch by the plain red cross (errorbars are included but due to the high S/N and the scale, they are not visible). The two red crosses are widely separated, displaying a motion non-compatible with a background source. Additionally, the motion of the planet is compatible with an orbital motion around the star.

Fig. A.2 Same as Fig. A.1, but for a background source presents in the IRDIS FoV. The two red crosses are compatible with each other, displaying a motion compatible with a background source.

Appendix B Example of a Gaia query to identify sources in the SPHERE FoV

Below, we present an S/N map of an observation around the star TYC_5164_0567_1, where three point-like sources were detected using PACO, labeled (1), (2), and (4). The Gaia query around this star (with a radius set at 8 arcsec to avoid missing sources) returns three detections indicated by the three green arrows. The properties of those detected stars are shown in Table B.1. The Gaia detection labeled (3) is outside the FoV of SPHERE, hence, it is not visible in the image. The detection labeled (1) has a similar parallax and proper motion as the primary, hinting to a bound nature, confirmed by Bonavita et al. (2022). Detection (2) is much further away than the primary, it is a background star. Finally, detection (4) is not seen by Gaia as its contrast is lower than the detection capabilities of Gaia, limited to sources with typically a G magnitude greater than 21.

Fig. B.1 S/N map of the star TYC_5164_0567_1 with three Gaia detected point-like sources labeled as (1), (2), and (4). The star is hidden behind the coronagraphic mask represented by the yellow circle at the center of the image. The three green arrows show direct Gaia detections.

Appendix C Color distribution of confirmed background sources

Histograms used to define the background exclusion area used in the CMDs in Fig. 9 and 10 are shown below for the IRDIFS configuration in Fig. C.1 and for the IRDIFS-EXT in Fig. C.2.

Fig. C.1 Histogram of the colors of confirmed background sources via proper motion by bin of magnitudes used to define the background exclusion area for the H23 CMD in Fig. 9.

Fig. C.2 Same as C.1, for the K12 CMD in Fig. 10.

Appendix D False positives identification and distribution

The identification of false positives is an important subject for any massive analysis such as this one. In summary, there are two complementary approaches to assess the statistical relevance of the claimed false alarm rate:

1. Comparing the empirical occurrence of detected FWHM-wide sources, identified as false positives, to their expected theoretical occurrence rate at the prescribed detection threshold.

2. Comparing the empirical distribution of the noise considered in the detection criterion, in the absence of astrophysical sources, to the expected theoretical distribution (i.e., a Gaussian with zero mean and unit variance). In case of PACO the relevant value is the estimated SNR per pixel.

Approach (1) is feasible when two epochs are available. It allows the clear identification of false positives, particularly if the second epoch is of higher quality, theoretically enabling the re-detection of sources with improved S/N. This eliminates random noise excursions that could lower the S/N and prevent the re-detection of a true astrophysical source. Where possible, we applied this method systematically on the survey, comparing the radial (angular) occurrence of false positives to the theoretical distribution (expected to follow a r² law if the noise in our SNR maps follows the same normal distribution in the speckle-noise dominated area and in the background-noise dominated area), see Sect. 5.2.4 and Figs. D.1-D.2. Our findings indicate a reasonable match between the two, especially given the heterogeneity of the survey, with a slight discrepancy at shorter angular separations within the speckle-dominated regime for IRDIS data. If we consider a pure thresholded analysis with PACO under 1" for IRDIS, we were able to identify 17 false positives in this area. This number is higher than the prediction in r² if the noise distribution is identical over the whole image. However, 8 of those 17 false positives are visible speckles identified thanks to the manual review of sources, leaving only 9 false positives, a number compatible with the expected one. The peculiar noise regime in the speckle-dominated area is therefore still responsible for a small increase in false positives, but this is partially compensated by reviewing detected sources.

Approach (2) is significantly more challenging, both to implement and then to interpret.

• Considering jointly the distribution of the detection criterion on the whole survey by stacking all S/N maps introduces biases, even with detected sources masked, from pathological cases (e.g., disks, smeared sources, or instrumental/reduction artifacts such as remanence leftovers). Such cases occur regularly across the hundreds of datasets we analyzed and result in false detections that are not solely representative of the noise distribution, as they arise from some form of signal rather than pure noise excursions. Even if not common, they can easily outnumber the small expected number of noise-induced excursions at four or five σ.

• Even stacking a carefully vetted subsample of good quality datasets will not directly yield the expected noise distribution because of the contribution to the histogram of astrophysical sources below the detection threshold and also above through stochastic amplification of a weak astrophysical source by random positive noise excursion. The latter is inherently unquantifiable, but the large number of background sources in our data above 5σ, combined with an expected increase of such astrophysical contamination at fainter fluxes (probing a larger volume of the Milky Way), ensure that these sources will bias the derived distribution.

As a result, the only reliable way to fully study the empirical false positive rate on the whole survey would be to reprocess all reductions using inverted parallactic angles. This would mitigates astrophysical signals by preventing constructive summation after derotation and stacking. However, systematically applying this method requires a large amount of computational resources and falls beyond the scope of this work given the tests already performed to ground the statistical relevance of our conclusions. Nonetheless, the following points should be emphasized:

• Comparisons of the empirical distribution of the detection criterion (with reversed parallactic rotation) to the theoretical distribution have already been conducted in several case-by-case studies (under good to median observing conditions) to validate the statistical relevance of PACO, see e.g. Figure 3 in Flasseur et al. (2018), Figures 2 and 4 in Flasseur et al. (2020b), and Figure 4 in Flasseur et al. (2020a). These studies consistently found a good match in all the cases analyzed.

• In addition, Chomez et al. (2023a) applied the same approach of inverted parallactic angles to estimate the false positive rate under good to median observing conditions. By extrapolating these results to the dataset analyzed in SHINE F400, we estimate 260 false positives for IRDIS and 162 for IFS. We acknowledge that this extrapolated value serves as a lower bound, as it accounts only for random noise excursions and not for the other sources of false positives discussed earlier. However, our final detections are not merely the result of applying a 5 σ threshold to the S/N maps. PACO incorporates several internal mechanisms to filter out pathological pixels and sources, and a manual review process during the analysis further eliminates additional spurious detections.

To summarize, we can reliably assess false positive rates under controlled observing conditions, and these estimates are representative of a standard SHINE observation. However, due to the variability in observing conditions and the presence of diverse astrophysical sources in our data, we acknowledge that the actual false positive rate on the whole survey is likely higher, and remains extremely challenging to evaluate more precisely. Nevertheless, this rate is mitigated by PACO’s automated mechanisms and the additional manual flagging of spurious detections during the analysis.

Fig. D.1 Polar plot showing the identified false positives on IRDIS and their radial distribution. The distribution is cut below 150 mas to ensure that the spectral leverage is enough to identify speckles.

Fig. D.2 Polar plot showing the identified false positives on IFS and their radial distribution. The distribution is cut below 150 mas to ensure that the spectral leverage is enough to identify speckles.

References

All Tables

All Figures

Fig. 1 Overall data quality for SHINE and snapSHINE. (a) Histogram of the star magnitudes (G, H, K) for the sample considered in this study. (b) Scatter plot showing the coherence time (τ0) against the seeing with the Strehl ratio color-coded. Data points without SPARTA Strehl information are colored in white. The grey box delimits the ESO 50% best observing conditions (seeing <1 as and τ0 > 3.2 ms) whereas the yellow box delimits the 30% best (seeing <0.8 as and τ0 > 4.2 ms). (c) Scatter plot showing the 1σ raw contrast at 0.5 arcsec (directly measured on the coronagraphic frame, before post-processing) with respect to the H-band (main IRDIS band used in SHINE) star apparent magnitude. The best contrast performances are achieved for the brightest stars with a high Strehl ratio. (d) Scatter plot displaying the Strehl ratio as a function of the star’s apparent magnitude in the G band (AO spectral channel). As expected, the performances of the AO system deteriorate for fainter stars.

In the text

	Fig. 2 Histogram of the total amount of parallactic angle rotation during the coronagraphic sequences for the SHINE observations considered in this survey.
In the text

	Fig. 3 Histogram of the total coronagraphic integration time during the sequences for the SHINE observations considered in this survey.
In the text

	Fig. 4 Contrast curves at 5σ obtained on the whole survey for both instruments, IRDIS in green and IFS in red. (a) IRDIFS setup. (b) IRDIFS-EXT setup. The contrast curves of the snapSHINE survey are represented in blue. The grey overlaid rectangles at small angular separations show the average inner working angle of the coronagraphic masks.
In the text

Fig. 5 Contrast comparison between the raw contrast (median average spectrally and spatially, dashed lines) and the PACO contrast, both at 5σ. Three bins of star magnitudes and observing conditions are considered. Individual observations are represented by the thin lines and the median for each observing condition by the thick lines. The gain gives the contrast improvement between the median raw curves and the median PACO curves, provided for different observing conditions and different star brightness. The average size of the coronagraph’s inner working angle is around 100 mas.

In the text

Fig. 6 Comparison between a SHINE and a snapSHINE observation of the same star. (a) S/N map of the star TYC_8634_1393_1 with the SHINE observation using PACO (τ = 5σ). (b) Residual map of the same star with the snapSHINE observation using NoADI. Four detections are present in the PACO S/N map, highlighted by the dashed squares. Only the three red ones are redetected in the snapSHINE observation (background sources). The other one in yellow is too faint to be detected by snapSHINE. A 15 × 15 pixel 2 D median filter was applied on the snapSHINE residual map and the detection was performed manually.

In the text

Fig. 7 Detection performance comparison between this work and the detection limits computed in Vigan et al. (2021). Left: comparison between the sensitivity achieved by Vigan et al. (2021) using classical post-processing algorithms (Langlois et al. 2021) on the F150 (dashed lines) compared to the sensitivity achieved using PACO on the same sample of stars (plain lines). The lower limit seen around 100 au at 0.5 MJup is due to the COND model, which does not provide inputs for lower masses. Right: probability increase (${\mathcal{P}}_{\text{PACO}}$(sma, MJup) – ${\mathcal{P}}_{\text{SPECAL}}$(sma, MJup)) between the PACO analysis and the one presented in Vigan et al. (2021).

In the text

	Fig. 8 S/N maps for (a) IRDIS and (b) IFS of HIP 74865 B. The coronagraph is represented by the yellow circle. The brown dwarf is clearly visible southeast of it. Maps are cropped at a 500 mas radius. The color bar represents the S/N.
In the text

	Fig. 9 H2–H3 color-magnitude diagram with every non confirmed bound source detected during the SHINE survey. Error bars for background sources are removed for clarity purposes. The background exclusion area is indicated by the grey area and corresponds to a 5σ exclusion area.
In the text

	Fig. 10 Same as Fig. 9, but for the K12 IRDIS filter and a 3σ exclusion area. The 5σ area is represented by the dotted lines. Sources between 3 and 5σ are labeled as AMBIGUOUS (see explanation in text Sect. 5.2.3).
In the text

	Fig. 11 Pie chart summarizing the number of physical sources in each category for the SHINE F400 survey. CC stand for Confirmed Companion, BCKG_GAIA for sources directly detected by GAIA and confirmed as background stars, BCKG_PM for background stars confirmed via proper motion analysis and BCKG_CMD, sources classified as background using a color magnitude diagram.
In the text

	Fig. 12 Cumulative normalized distribution of the separation of the sources according to their classes. Classes colors are identical to those in Figs. 9, 10, and 11.
In the text

Fig. A.1 Common proper motion plot for the planet HD 95086 b. The position of the planet at the first epoch is represented by the black cross, the expected position at the second epoch, if background, by the red hollow cross, and the measured position at the second epoch by the plain red cross (errorbars are included but due to the high S/N and the scale, they are not visible). The two red crosses are widely separated, displaying a motion non-compatible with a background source. Additionally, the motion of the planet is compatible with an orbital motion around the star.

In the text

	Fig. A.2 Same as Fig. A.1, but for a background source presents in the IRDIS FoV. The two red crosses are compatible with each other, displaying a motion compatible with a background source.
In the text

	Fig. B.1 S/N map of the star TYC_5164_0567_1 with three Gaia detected point-like sources labeled as (1), (2), and (4). The star is hidden behind the coronagraphic mask represented by the yellow circle at the center of the image. The three green arrows show direct Gaia detections.
In the text

	Fig. C.1 Histogram of the colors of confirmed background sources via proper motion by bin of magnitudes used to define the background exclusion area for the H23 CMD in Fig. 9.
In the text

	Fig. C.2 Same as C.1, for the K12 CMD in Fig. 10.
In the text

	Fig. D.1 Polar plot showing the identified false positives on IRDIS and their radial distribution. The distribution is cut below 150 mas to ensure that the spectral leverage is enough to identify speckles.
In the text

	Fig. D.2 Polar plot showing the identified false positives on IFS and their radial distribution. The distribution is cut below 150 mas to ensure that the spectral leverage is enough to identify speckles.
In the text