# The COBREX archival survey: improved constraints on the occurrence rate of wide-orbit substellar companions.
## I. A uniform re-analysis of 400 stars from the GPIES survey

V. Squicciarini, J. Mazoyer, A. -M. Lagrange, A. Chomez, P. Delorme, O. Flasseur, F. Kiefer

► To cite this version:

V. Squicciarini, J. Mazoyer, A. -M. Lagrange, A. Chomez, P. Delorme, et al.. The COBREX archival survey: improved constraints on the occurrence rate of wide-orbit substellar companions. I. A uniform re-analysis of 400 stars from the GPIES survey. Astronomy & Astrophysics - A&A, 2025, 693 (A54), 10.1051/0004-6361/202452310. [hal-04778163](https://hal.science/hal-04778163v1)

HAL Id: hal-04778163

https://hal.science/hal-04778163v1

Submitted on 12 Nov 2024

**HAL** is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire **HAL**, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# The COBREX archival survey: improved constraints on the occurrence rate of wide-orbit substellar companions

## I. A uniform re-analysis of 400 stars from the GPIES survey

V. Squicciarini<fnref n="1" />, J. Mazoyer<fnref n="1" />, A.-M. Lagrange<fnref n="1" />,<fnref n="2" />, A. Chomez<fnref n="1" />,<fnref n="2" />, P. Delorme<fnref n="2" />, O. Flasseur<fnref n="3" />, and F. Kiefer<fnref n="1" />

November 8, 2024

## ABSTRACT

**Context.** Direct imaging (DI) campaigns are uniquely suited to probing the outer regions around young stars and looking for giant exoplanet and brown dwarf companions, hence providing key complementary information to radial velocity (RV) and transit searches for the purpose of demographic studies. However, the critical 5-20 au region, where most giant planets are thought to form, remains poorly explored, lying in-between RV and DI capabilities.

**Aims.** Significant gains in detection performances can be attained at no instrumental cost by means of advanced post-processing techniques. In the context of the COBREX project, we have assembled the largest collection of archival DI observations to date in order to undertake a large and uniform re-analysis. In particular, this paper details the re-analysis of 400 stars from the GPIES survey operated at GPI@Gemini South.

**Methods.** Following the pre-reduction of raw frames, GPI data cubes were processed by means of the PACO algorithm. Candidates were identified and vetted based on multi-epoch proper motion analysis – whenever possible – and by means of a suitable color-magnitude diagram. The conversion of detection limits into detectability maps allowed for an estimate of unbiased occurrence frequencies of giant planets and brown dwarfs.

**Results.** Deeper detection limits were derived compared to the literature, with up to a twofold gain in minimum detectable mass compared to previous GPI-based publications. Although no new substellar companion was confirmed, we identified two interesting planet candidates awaiting follow-up observations. We derive an occurrence rate of $1.7_{-0.7}^{+0.9}\%$ for $5 M_{\text{Jup}} < m < 13 M_{\text{Jup}}$ planets in $10 \text{ au} < a < 100 \text{ au}$, that raises to $2.2_{-0.8}^{+1.0}\%$ when including substellar objects up to $80 M_{\text{Jup}}$. Our results are in line with the literature, but come with lower uncertainties thanks to the enhanced detection sensitivity. We confirm the finding, hinted at by previous studies, of a larger occurrence of giant planets around BA hosts compared to FGK stars; moreover, we tentatively observe a smaller occurrence of brown dwarf companions around BA stars, although larger samples are needed to shed light on this point.

**Conclusions.** While waiting for the wealth of data expected from future instrument and facilities, valuable information can still be extracted from existing data. In this regard, a complete re-analysis of SPHERE and GPI data is expected to provide the most precise demographic constraints ever provided by imaging.

**Key words.** Planets and satellites: detection – Planets and satellites: gaseous planets – brown dwarfs – Planetary systems – Techniques: high angular resolution

## 1. Introduction

Bolstered by more than 5000 confirmed detections to date, the

exoplanet field has become mature enough to accompany the

still thriving detection-oriented endeavor with follow-up stud-

ies aimed at shedding light on key questions related to the ori-

gin, the prevalence, and the architecture of planetary systems.

By unveiling statistical trends in the measured physical, orbital

and star-related properties of the exoplanet population, exoplanet

demographics seeks to connect theory and observation, in order

to fully understanding the physical processes underlying planet

formation (Biazzo et al. 2022).

The census of known exoplanets currently spans about four magnitudes in mass and in semi-major axis<fnref n="1" />. No single detection

method is adequate to probe such a large extent of the parame- ter space: it is through the combination of the different meth-

ods, each optimized for detection inside a specific niche, that

the large-scale picture can be unveiled and reconstructed (see,

e.g., Gratton et al. 2023, 2024). However, obtaining a complete

and unbiased blend from heterogeneous ingredients is hindered

by factors such as inconsistent detection criteria, completeness

and false positive assessment, uncertainty quantification, neglect

of underlying selection or observational biases (Gaudi et al.

<footnote n="1">¹ LESIA, Observatoire de Paris, Université PSL, CNRS, Sorbonne Université, Université Paris Cité, 5 place Jules Janssen, 92195 Meudon, France
e-mail: vito.squicciarini@obspm.fr</footnote>

<footnote n="2">² Univ. Grenoble Alpes, CNRS-INSU, Institut de Planétologie et d'Astrophysique de Grenoble (IPAG) UMR 5274, Grenoble, F-38041, France</footnote>

<footnote n="3">³ Centre de Recherche Astrophysique de Lyon (CRAL) UMR 5574, CNRS, Univ. de Lyon, Univ. Claude Bernard Lyon 1, ENS de Lyon, F-69230 Saint-Genis-Laval, France</footnote>

<footnote n="1">¹ Empirically estimated based on the Extrasolar Planet Encyclopedia: http://www.exoplanet.eu/.</footnote>
2021). Whenever two different methods can be simultaneously employed, their complementarity allows better characterizing individual objects (see, e.g., Gandolfi et al. 2017; Bonnefoy et al. 2018; Bourrier et al. 2018; Lacedelli et al. 2021; Kuzuhara et al. 2022; Philipot et al. 2023, Lagrange et al. under review) and strengthening the statistical trends emerging in each of the methods (Rogers 2015; Santerne et al. 2016). In the cases where different techniques probe instead different separations within the same system, the joint analysis opens up exquisite dynamical and formation studies (see, e.g., Covino et al. 2013; Bryan et al. 2016; Zhu & Wu 2018).

Radial velocity (RV) surveys have provided invaluable constraints on the physical and orbital properties of giant planets up to ~ 5 au (Wolthoff et al. 2022; Rosenthal et al. 2024). Yet, the reliability of RV trends for larger separations has been questioned (Lagrange et al. 2023), and the predicted yields for direct imaging (DI) surveys based on extrapolations of RV results have been shown to be too optimistic (see, e.g., Cumming et al. 2008; Dulz et al. 2020). On the other hand, direct imaging (DI) is mostly sensitive to young giant planets in wide ($a \ge 20$ au) orbits, providing access to the scarcely studied outskirts of planetary systems. Starting from 2004 (Chauvin et al. 2004), direct imaging has discovered ~ 30 planets ($M < 13 M_{Jup}$) (Zurlo 2024), including iconic systems like the disk-enshrouded PDS 70 (Keppler et al. 2018), the ~ 20-Myr-old $\beta$ Pictoris (Lagrange et al. 2009), 51 Eridani (Macintosh et al. 2015), AF Leporis (Mesa et al. 2023; De Rosa et al. 2023; Franson et al. 2023), and the four-planet HR 8799 (Marois et al. 2008). These detections are the main outcome of large blind surveys targeting tens (e.g. MASSIVE, Lannier et al. 2016; SEEDS, Uyama et al. 2017; LEECH, Stone et al. 2018) or hundreds of stars (e.g. NICI-PCF, Liu et al. 2010; IDPS, Galicher et al. 2016; ISPY-NACO, Launhardt et al. 2020). The forefront of DI surveys, enabled by the exquisite performances of imagers and integral field spectrographs, coupled with extreme AO systems mounted on 8-m-class telescopes, is currently represented by the 400-star SHINE (Chauvin et al. 2017) and the 600-star GPIES (Nielsen et al. 2019) surveys.

By constraining the overall frequency and the properties of wide-separation giant planets, DI studies are expected to enable a thorough comparison with concurrent formation models (see, e.g., Bowler 2016; Vigan et al. 2021); orbital properties, for instance, shed light upon their formation and dynamical evolution (Bowler et al. 2020); the dependence of frequency on stellar mass provides clues about the initial state of the disk and the formation mechanisms at play (Nielsen et al. 2019; Janson et al. 2021). However, despite years of extensive searches, it is still not clear whether the main formation channel for the observed wide-orbit population be core accretion (CA; Pollack et al. 1996; Mordasini et al. 2009), the bottom-up process responsible for the formation of planets in the Solar System, or rather a top-down star-like scenario like gravitational instability (GI; Boss 1997; Vorobyov 2013). While an interplay between the two scenarios is deemed to be favored by empirical parametric models (Reggiani et al. 2016; Vigan et al. 2021) and direct comparison with synthetic planet populations (Vigan et al. 2021) alike, understanding in an unambiguous way how each known companion was formed is still beyond reach. The large uncertainties still existing in the interpretation of the observed picture can be attributed, at least partially, to the fact that the critical 5-20 au region, where most giant planets are thought to form, remains poorly explored being exactly in-between current RV and DI capabilities.

Under given observing conditions, the final performances attainable by a high-contrast imaging observation are dictated

both by instrumental (e.g., the telescope, the science instrument, the performance of adaptive optics and coronagraphs) and post-processing components (the algorithms applied to science images to decrease the level of systematic and random noise) (Galicher & Mazoyer 2024). Depending on observing conditions, stellar brightness and angular separation, state-of-the-art instruments such as the Spectro-Polarimetric High-Contrast Exoplanet Research (SPHERE; Beuzit et al. 2019) and the Gemini Planet Imager (GPI; Macintosh et al. 2014) typically achieve raw planet-to-star contrasts as low as $10^{-3} - 10^{-5}$ (Poyneer et al. 2016; Courtney-Barrer et al. 2023). On the instrumental side, 30-m-class telescopes and space-borne coronagraphic instruments are expected to bring about a major leap forward for the field in the next decade (see, e.g., Kasdin et al. 2020; Kasper et al. 2021), whereas upgrades of existing instruments such as SPHERE+ (Boccaletti et al. 2022) and GPI 2.0 (Chilcote et al. 2018) are going to represent the forefront in the medium term; on the reduction side, advanced post-processing algorithms have been already shown to increase the contrast by as much as two orders of magnitudes compared to pre-reduced data. Therefore, the developments of more powerful reduction techniques can greatly increase detection capabilities working on observations that already exist (see, e.g., Currie et al. 2023).

In the framework of the COupling data and techniques for BReakthroughs in EXoplanetary systems exploration (CO-BREX) project, we collected more than a thousand archival SPHERE and GPI observations, assembling the largest exoplanetary direct imaging survey to date, with the aim of re-reducing them in a uniform and self-consistent way. The results of the full re-reduction of the SHINE survey are illustrated in Chomez et al. (2024). In this work, we present the re-reduction of 400 stars coming from GPIES. Despite being the largest DI observational campaign to date, just two new substellar objects were discovered during the survey: one planet (51 Eri b, Macintosh et al. 2015) and one brown dwarf (HR 2562 B, Konopacky et al. 2016). A statistical analysis of the first 300 stars was performed by Nielsen et al. (2019, hereafter N19). By combining the two surveys, it will be possible to obtain the tightest constraints to date on the occurrence of wide-orbit giant planets, hence providing an ideal test-bed to scrutinize planet formation models.

This paper is organized as follows: after laying out the selection criteria for the sample and the corresponding observations (Section 2.2), and uniformly deriving the stellar parameters of interest (Section 2.3), we describe in detail the process of data reduction (Section 2.4). Section 3 presents the results of the analysis, namely companion candidates and completeness maps. The derived occurrence rates are presented and discussed in Section 3.4. A thorough comparison with the literature is the subject of Section 4. Finally, in Section 5, we summarize the results of this work.

## 2. Data

### 2.1. Raw data collection

The observations considered in this work were collected between 2013 and 2020 by means of GPI at the Gemini South telescope. GPI is an integral-field spectrograph (IFS) with low spectral resolution (~50; Maire et al. 2014), operating in the wavelength range [0.97-2.40] $\mu$m. As the vast majority of GPI observations were gathered in the H band ([1.5-1.8] $\mu$m) (Ruffio et al. 2017) – other bands being mostly used for characterization purposes – we decided to restrict our query to H-band observations.
We downloaded therefore all H-band raw frames from GPI that are publicly available on the Gemini archive<fnref n="2" /> (~30000 frames)<fnref n="3" />. We neglected observational sequences with 5 or less frames (typically corresponding to ≤ 5-min exposure times), and observations of stars only taken for calibration purposes (easily identifiable through their program ID). The reason behind this choice is twofold: on the one hand, the known multiplicity of these stars is expected to detrimentally affect the performances attainable by post-processing; on the other hand, selection criteria for these stars are different from those of science targets, thus inducing a bias when interested in statistical considerations.

The preliminary sample obtained in this way (hereafter GPI database) is composed of 852 sequences for 655 stars. Most of the observations (715/852) within the database were collected in the course of GPIES, and additional 10 sequences are describable as follow-up observations of interesting stars from the campaign. Intertwined to GPIES observations, the remaining 127 sequences were gathered over the lifetime of the instrument as part of other scientific programs. The following Section 2.2 will elucidate how the final stellar sample was assembled, and the criteria that a sequence had to meet in order to be included in the corresponding sample of observations.

## 2.2. Sample definition

Like any other direct imaging search to date, GPIES was constructed by looking for young stars in the solar neighborhood; the reason lies in the fact that recently formed exoplanets and brown dwarfs are brighter and hotter than mature objects of the same mass due to residual formation heat, yielding a significantly more favorable planet-to-star contrast in near infrared bands. In particular, the stellar sample was assembled by merging lists of members of young moving groups from the literature (de Zeeuw et al. 1999; Zuckerman et al. 2001, 2011) with close (<100 pc) stars selected for large X-ray emission. Echelle spectra were obtained for ~2000 stars to further identify additional young stars based on lithium abundance and chromospheric activity (see Section 2.3 for details). After removing apparent binaries with angular separation $\epsilon$ [0.02", 3"] and $\Delta mag < 5$ (both before and during the campaign), and accounting for some new association members, a sample of 602 stars was finally obtained (Nielsen et al. 2019).

Due to the decommissioning of GPI – currently undertaking major upgrades to become GPI 2.0 (Chilcote et al. 2020) – in early 2020, the GPIES survey was never completed (~ 10% of the stars lack observations). It is thus vital to ascertain whether the sample of observed stars be an unbiased extraction of the full sample. An automated target-picker was employed to suggest the best targets for every telescope night, as a function of both observing conditions and stellar age/distance (McBride et al. 2011); however, we verified through a Kolmogorov-Smirnov test ($\alpha$ = 0.05) that the age and distance distributions of the first 300 stars (those from N19) are compatible with those of the full sample of observed stars ($p_{age} \approx 1$, $p_{dist} = 0.051$). Additionally, stars with known companions were not prioritized in their first-epoch observation (Nielsen et al. 2019). We can thus confidently maintain that the available GPIES observations are not affected by

selection biases, making the sample suitable for statistical studies.

The definition of the final sample was based on a combination of observational constraints and physical constraints on stellar properties. As regards the former aspect, a minimum amount of parallactic angle rotation $\Delta PA$ ~ 10° is required to enable efficiently using angular differential imaging (ADI) during post-processing (Marois et al. 2006): we conservatively adopt a minimum rotation of 12° in order to exploit angular differential imaging (ADI) during post-processing<fnref n="4" />. A single observation with extremely bad seeing was removed. With respect to the latter, we only retained stars for which youth can be established with reasonable confidence (see Section 2.3).

Given our ignorance about the selection criteria adopted for stars from non-GPIES programs, we decided not to consider them for the purpose of this paper. However, we retained non-GPIES observations of GPIES stars as valuable follow-up epochs for promising point-source candidates. The final sample employed throughout this work consists of 400 stars (515 sequences).

## 2.3. Stellar parameters

The knowledge of stellar ages is pivotal to a meaningful interpretation of direct imaging campaigns, as a large degeneracy exists between age and mass – let alone additional parameters like metallicity or a planet’s formation history – for substellar objects for which only photometric data are available (see, e.g., Spiegel & Burrows 2012).

Our primary age diagnostics is provided by kinematic membership to young associations and moving groups (hereafter YMGs). Starting from Gaia DR3 (hereafter Gaia; Gaia Collaboration et al. 2023) data, we used BANYAN $\Sigma$ (Gagné et al. 2018b) to classify a star as a member of a YMG if the associated membership probability $p > 90\%$. A second indicator was represented by the ages obtained by N19: we stress that the underlying data are not public and the derived ages, not equipped with error bars, are only available for the 300 stars presented in that study. Finally, the ages for 21 additional stars – that are not members of YMG nor targets of N19 – could be recovered after cross-matching our sample with SHINE (Desidera et al. 2024), whose thorough analysis builds upon a manifold variety of indicators (isochrones, YMG membership, activity, lithium abundance). For ages based on N19, that come with no associated uncertainty, we adopt a constant fractional uncertainty of 25%, empirically tuned to match the typical fractional uncertainty for SHINE stars.

Individual stellar parameters for each star were obtained by means of MADYS<fnref n="5" /> (v1.2, Squicciarini & Bonavita 2022), a tool for (sub)stellar parameter determination based on the comparison between photometric measurements and isochrone grids derived from theoretical (sub)stellar models. Assuming the ages described above, photometry from Gaia DR3 and 2MASS (Skrutskie et al. 2006) – corrected by extinction by integrating the 3D map by Leike et al. (2020) along the line of sight – was compared to the last version of non-rotating, solar-metallicity PARSEC isochrones (Nguyen et al. 2022).

The distance, age and spectral type of the 400 stars are shown in Figure 1. The full collection of the derived properties is sum-

<footnote n="2">² https://www.cadc-ccda.hia-iha.nrccnrc.gc.ca/en/gemini/ and https://archive.gemini.edu/</footnote>

<footnote n="3">³ We could not retrieve science commissioning data from 2014 from the old Gemini website (Macintosh et al. 2014). However, the targeted stars were only known hosts of exoplanets or disks (HR 8799, HR 4796, HD95806) that were later re-observed during GPIES.</footnote>

<footnote n="4">⁴ The distribution of $\Delta PA$ across the database shows a bump for $\Delta PA > 12$. Whereas just 12 sequences have $10° < \Delta PA \le 12°$, 49 observations are found inside the $12° < \Delta PA \le 14°$ bin.</footnote>

<footnote n="5">⁵ https://github.com/vsquicciarini/madys</footnote>
<figure><img src="image_2.png" /><figcaption>Fig. 1: Age of the final stellar sample as a function of distance. The color scale labels different spectral types. Kernel density estimates for the distribution of the two properties are provided on top and to the right of the main plot.</figcaption></figure>

259 marized in Table A, while the ages adopted for YMGs are pro- 260 vided in Table A.2.

2.4. Data reduction

262 2.4.1. Preliminary steps

263 The pre-reduction of raw GPI data is performed in two steps: the
264 goal of the first step is to build a 3D data cube (x, y, λ) starting
265 from the 2D image acquired in the detector plane; the collection
266 of the derived data cubes is then stacked into a final 4D data cube
267 (x, y, λ, t).

268 As routinely done for GPI data, we employed the version
269 1.6.0 of the GPI Data Reduction Pipeline (hereafter DRP; Perrin
270 et al. 2014, 2016; Wang et al. 2018) to handle the first part of
271 the pre-reduction. After subtracting a dark frame, bad pixels in
272 the image are substituted by interpolated values. The mapping
273 of the ~ 37000 small spectra created on the detector by GPI’s
274 lenslet array to the corresponding spaxels of the 3D data cube is
275 determined by means of master wavelength calibrations based on
276 Xenon or Argon lamps, conveniently corrected for mechanical
277 offsets induced by flexure (Wolff et al. 2014). The signal of each
278 spectrum can be thus extracted and stored into the corresponding
279 spaxel (Draper et al. 2014). Subsequent steps correct for small
280 distortion effects of the field of view and for halos induced by
281 residual atmospheric turbulence.

282 Thanks to a square grid embedded within the pattern of
283 the apodizer, a diffraction pattern of satellite spots (hereafter
284 satspots) – attenuated images of the star – is created in the image.
285 The four first-order satspots, symmetrically situated at ~ 20λ/D
286 from the star, serve three different purposes: 1) to recenter the
287 frames, by locating the position of the occulted star; 2) to cali-
288 brate the flux level of each pixel in the science image; 3) to build
289 a model of the off-axis point spread function (PSF)<fnref n="6" /> (Wang et al.
290 2014). The final operations of the DRP deal with the astrometric
291 and photometric characterization of the satellite spots, fitted by
292 a Gaussian PSF template.

The second step of the pre-reduction deals with stacking up 293 and recentering frames to build the final 4D data cube. In ad- 294 dition to the data cube, three files are created: 1) a 4D PSF; 2) 295 a wavelength vector; 3) a parallactic angle vector, that indicates 296 the rotation of the field of view during the sequence. Indeed, GPI 297 observations are operated in a pupil-stabilized mode, i.e. with no 298 derotator, so as to allow using ADI-based post-processing algo- 299 rithms (see, e.g., Ruffio et al. 2017). 300

This stage is achieved using pyklip (Wang et al. 2015a). 301 Compared to the 2.6 version, we introduced slight modifica- 302 tions<fnref n="7" /> to create output files formatted in a SPHERE-like way as 303 regards the data format, the PSF, and data cube flux, and the FOV 304 orientation (East to the left). In this way, we ensured the har- 305 monization of the future SPHERE+GPI sample while smooth- 306 ing the I/O integration with the post-processing algorithm (Sec- 307 tion 2.4.2). In addition to computing the image center through 308 satspot pattern, pyklip estimates satspot flux in a more precise 309 way than the DRP, a crucial step for photometric characteriza- 310 tion purposes, and empirically recomputes the wavelength vec- 311 tor based on the satspot-to-center separation (which scales with 312 λ). We visually checked the goodness of the result for all our im- 313 ages; whenever a specific satspot was – due to intrinsic dimness 314 or systematic problems – not properly fitted, inducing centering 315 offsets in one or more frames, we used a specific option of pyklip 316 to ignore it during recentering. 317

As already mentioned, satellite spots are faint images of the target star; the flux ratio between a satspot and the star, or grid ratio, was determined by Wang et al. (2014) to be Δm = 9.4 ± 0.1 mag through on-sky observations. The ~ 10% uncertainty on the grid ratio turns out to be one of the main factors in the total error budget of GPI spectrophotometry.

We performed several tests to quantify the reliability of the wavelength solution, the image centering and the photometric calibration of the PSF. The accuracy of the DRP wavelength solution was estimated by Wolff et al. (2014) to be 0.032% in H-band, well below the 1% accuracy needed to achieve a spectral characterization uncertainty < 5%. As regards the wavelength precision, we collected for every sequence the satspot positions estimated by the DRP, computed separations from pyklip’s frame centers, ξ(λ, t, s), then averaged over the temporal axis t and the satspot axis s to obtain ˙ξ(λ). We computed the ratio η = ξ(λ)/λ for every sequence, a value that ought to be constant. The 50th, 16th and 84th percentiles of the distribution across the sequences yields < η > = 37.921<sup>+0.056</sup><sub>-0.046</sub> px/µm, corresponding to a precision of 0.15%.

As regards centering precision, propagation of random uncer-
tainties in pyklip yields a centering precision along one axis
σ<sub>c,x</sub> = 0.04 px (that is, ~ 0.06 px in 2D), comparable to the one
reported in Wang et al. (2014). However, we identified a system-
atic deviation of the satspot pattern shape from a square (which is
an underlying assumption of pyklip’s centering algorithm): the
difference between the centers computed from doublets of op-
posite satspots (x<sub>c,13</sub>, y<sub>c,13</sub>) and (x<sub>c,24</sub>, y<sub>c,24</sub>), stable over time, is
ΔS = √(x<sub>c,13</sub> - x<sub>c,24</sub>)<sup>2</sup> + (y<sub>c,13</sub> - y<sub>c,24</sub>)<sup>2</sup> = 0.39 ± 0.04 px. We
consider the true center to be distributed according to a uniform
distribution between (x<sub>c,13</sub>, y<sub>c,13</sub>) and (x<sub>c,24</sub>, y<sub>c,24</sub>). The final cen-
tering error σ<sub>c</sub> is therefore:

$$
\sigma_c = \sqrt{\left(\frac{\Delta S}{\sqrt{12}}\right)^2 + 4 \cdot (\sigma_{c,x})^2} \approx 0.15 \text{ px} \sim 2 \text{ mas.} \quad (1)
$$

<footnote n="6">⁶ Unlike SPHERE, no PSF exposure is taken before and/or after the
scientific observation in GPI.</footnote>

<footnote n="7">⁷ url of GitHub repository, to be added upon acceptance</footnote>
This value was consistently employed when propagating astro-
metric uncertainties of detected sources.

Finally, we adopted a platescale of 14.161 ± 0.021 mas px⁻¹
(De Rosa et al. 2020a), assuming it to be stable over time (Tran
et al. 2016). With respect to the north offset angle, we used a
time-varying value following the prescriptions indicated in Table
4 from De Rosa et al. (2020a).

2.4.2. Post processing – PACO

Pre-reduced datasets were processed in the COBREX Data
Center, an improved version of the High-Contrast Data Cen-
ter (HC-DC<sup>8</sup>, formerly SPHERE Data Center, Delorme et al.
2017). Prompted by the promising preliminary results presented
in Chomez et al. (2023), we decided to process our archive
by means of the PAtch-COvariance algorithm (PACO; Flasseur
et al. 2018) in its robust angular and spectral differential imaging
(ASDI) mode (Flasseur et al. 2020a,b).

PACO is a post-processing algorithm that employs ASDI to model the spatial and temporal fluctuations of the image background inside small patches through a combination of weighted multivariate (i.e., accounting for the spatio-spectral correlations of the speckles field) Gaussian components. Extensive testing proved that the resulting SNR map is distributed as a normalized Gaussian N(0, 1), hence naturally providing a statistically grounded detection map upon which > 5σ detections can be identified at a controlled false alarm rate (e.g., at 5σ significance level). The algorithm was shown to be photometrically accurate and robust to false positives, and to outperform reduction methods that are routinely employed for SPHERE (Chomez et al. 2023). For these reasons, PACO was appointed by the SHINE consortium as the main reduction algorithm for the final analysis of the whole survey (Chomez et al. 2024).

In addition to this, the ASDI mode of PACO uses vectors
of spectral weights (hereafter spectral priors) to maximize the
detection capability of candidate sources exhibiting physically
representative substellar spectra. As detailed in Chomez et al.
(2023), for every star we generated 20 such priors starting from
exoplanet spectra from the ExoREM library (Charnay et al.
2019) (T<sub>eff</sub> ∈ [400, 2000] K) and suitable stellar spectra from
the BT-Nextgen AGSS2009 library (Allard et al. 2011).

As in Chomez et al. (2023), extensive injection tests were
performed on GPI datasets in order to ensure the reliability of 5σ
detection limits, an output provided by PACO after the reduction
(Flasseur et al. 2020a). After randomly picking a sample of 10
sequences, 12 synthetic sources were evenly injected in each ob-
servation's FOV; the mean flux of each source was set equal to
the 5σ detection limits estimated by PACO at the corresponding
coordinates. The whole process was repeated three times, vary-
ing the input spectrum – a flat contrast spectrum, a T-type spec-
trum and an L-type spectrum<fnref n="9" /> – of injected sources, yielding a
total 360 injected sources. The median SNR of the recovered
sources is 5.1, with little variation with spectral type, confirming
the statistical reliability of the contrast and detection confidence
estimated by PACO and underlying the statistical analysis.

Figure 2 shows the final performances attained by the PACO
reduction as 5σ detection limits. It is possible to notice that the
usage of physically-motivated spectral priors does indeed en-

<figure><img src="image_3.png" /></figure>

Fig. 2: 5σ detection limits obtained with PACO. Individual curves are plotted in gray. The median curve is plotted as a light blue solid line, the dashed lines representing the 16% and 84% percentiles of the curve distribution. The orange and green solid lines indicate the median detection limits assuming a T-type and an L-type spectral prior, respectively. The gray box marks the inner working angle of the coronagraphic mask.

hance detection capabilities. However, in order to easily allow
for comparisons with reductions performed using different algo-
rithms, we conservatively adopt in the following analysis a flat
spectral prior, that is a combination of spectral channels assum-
ing that any source has the same spectral energy distribution as
its star; this is equivalent to standard SDI-based algorithms.

2.4.3. Post processing – cADI

One of the underlying assumptions behind PACO is that the
spatial and temporal fluctuations of noise inside patches are
much stronger than the additional contribution given by physi-
cal sources happening to cross the patch itself during the expo-
sure (Flasseur et al. 2018). The assumption breaks down when
a bright source, such as a stellar companion, is present. In other
words, the algorithm is optimized to detect faint sources but can
severely subtract, or even cancel, very bright sources in the de-
rived SNR map. In order to complete the census of sources at
the bright end, we developed a custom routine based on classi-
cal angular differential imaging (cADI; Marois et al. 2006) and
performed a uniform reduction of the archive. After computing
the pixel-wise median frame of the exposure sequence, the rou-
tine subtracts it from every frame, then de-rotates the frames and
sums them up both temporally and along wavelength. The 4D
PSF is stacked along the temporal axis to build a 3D (x, y, λ)
PSF, which is then fitted by a 2D Gaussian model. The reduced
map is finally normalized by the peak of the PSF model so as
to translate it in contrast units. Detections are automatically per-
formed on the derived map by computing the variance across
annuli, centered on the target star, of width equal to 1 px, and
finding the pixels of the map beyond a certain threshold level
κ (expressed in noise standard-deviation units). In a subsequent
step, a more precise characterization by fitting the PSF model
provides the astrometry and photometry of each source.

Because of the simplicity of the noise-reduction approach,
the SNR distribution in any annulus of a given contrast maps
usually shows large deviations from Gaussianity. On the one

<footnote n="8">⁸ https://sphere.osug.fr/spip.php?rubrique16&lang=en</footnote>

<footnote n="9">⁹ Input contrast spectra, created using the same procedure as spectral priors, are arbitrarily defined using the following parameters: (T<sub>eff</sub>, log g, Z/Z<sub>☉</sub>, CO) = (1000 K, 4.0, 1.0, 0.5) for the T type and (T<sub>eff</sub>, log g, Z/Z<sub>☉</sub>, CO) = (1900 K, 3.0, 10.0, 0.6) for the L type.</footnote>
<figure><img src="image_4.png" /><figcaption>Fig. 3: Comparison between PACO and cADI performances. Sources only detected by PACO are shown as blue circles, while sources only seen through cADI are indicated as red diamonds. Common sources are plotted as blue squares with a red edge.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 4: Example of proper motion diagram. The astrometric displacement of the candidate around HD 84330B between the first and the second epoch is compatible with a background source with null motion (empty star). A bound object would have been in a position close to that marked by the filled star and within the boundaries allowed by a Keplerian motion.</figcaption></figure>

side, this issue implies that high thresholds ($κ \ge 20$) had to be adopted to ensure the stability of the detection step; on the other side, the poor robustness against outliers intrinsically prevents one from precisely defining a statistically grounded detection threshold.

Visual inspection of all the maps ensured the reliability of the detections; given the above-mentioned caveats and the neglect of a correction for signal self-subtraction, the derived photometry will only be used to characterize the stellar companions presented in Section D.

## 3. Results

### 3.1. Exoplanet candidates

The PACO reduction of our sample yielded 91 detected sources. This number does not include a few false positives that could be recognized and removed (see Section 3.2). 11 additional sources were detected through cADI. 62 sources are detected by both methods, ensuring the overlap of the respective dynamical ranges. Figure 3 shows the sources detected by the two methods. Astrometric and photometric details for all the candidates are provided in Table C.1.

Candidate companions in DI observations are always seen as unresolved point-like sources, and no information on their distance can be discerned from a single observation; in other words, it is not clear a priori if a source is physically bound to the target star or is instead a distant background star that happens to be projected close to the target. If two or more epochs are available, the differential motion between the foreground target star (and the objects bound to it) and faraway background stars can be disentangled (Figure 4).

Whenever more than one observation was available in our sample, or if additional epochs from SPHERE could be recovered, it was possible to ascertain the proper motion of the candidates: 57 sources from the PACO reduction and 2 sources only detected with cADI were ruled out as background contaminants in this way.

If only a single observation was available, or if detection limits allowed for detection in just one epoch, we adopted an alter-

native vetting criterion that exploits color-magnitude diagrams (CMDs) to identify sources showing similar colors to known imaged planets and to set them apart from background sources. It might be argued that, given the availability of contrast spectra, a spectrum-based classification could be employed: however, we argue that such a method would be highly sensitive to both random uncertainties and the ignorance about the amount of interstellar extinction to be adopted for background-star spectra. Conversely, the photometric method based on CMDs has already been shown to be highly reliable for absolute magnitudes $H \ge 15$ mag, using an unprecedented sample of ~2000 of confirmed astrophysical background sources found in SPHERE data (Chomez et al. 2024). This usage of the CMD has been introduced by the SHINE consortium (Chauvin et al. 2017) and its construction is fully detailed in Bonnefoy et al. (2018). This tool has already been used to efficiently classify some of the sources detected in the first part of the SHINE survey (Langlois et al. 2021). As a first step, the H-band spectrum of each target star was estimated by means of synthetic stellar spectra from the BT-Nextgen AGSS2009 library<fnref n="10" /> (Allard et al. 2011), adequately degraded to match the spectral resolution of GPI. The best-matching synthetic spectrum was identified as the closest in effective temperature; the latter was empirically estimated as the median value across all literature measurements found in VizieR (Ochsenbein et al. 2000). Contrast spectra from candidate sources detected with PACO<fnref n="11" /> could thus be turned into physical spectra by multiplying them by their corresponding primaries spectra. We convolved these spectra with SPHERE H2 and H3 filters to derive synthetic H2 and H3 photometry for all our candidates; in other words, GPI spectroscopy was turned into SPHERE-like photometry both to exploit the CMD vetting

<footnote n="10">¹⁰ The library, available at http://svo2.cab.inta-csic.es/theory/newov2/index.php?models=bt-nextgen-agss2009, is defined by the following astrophysical parameters: log g[cm s⁻²] = 4.5, log Z/Z⊙ = 0, alpha enhancement = 0.</footnote>

<footnote n="11">¹¹ The status of all cADI candidates but one could be confirmed through dynamical arguments; the remaining one is too bright to allow for the CMD test.</footnote>
<figure><img src="image_6.png" /><figcaption>Fig. 5: CMD of the companion candidates detected in this work. Overplotted to known substellar objects (white squares), background stars are represented as yellow stars if identified through proper motion analysis, or as blue circles if recognized via their color. Ambiguous sources are marked as red crosses. The exclusion area (gray) is defined by the two dashed lines. The two promising candidates (from left to right: C1 (HD 24072), C2 (HIP 78663)) are indicated as red dots.</figcaption></figure>

method and to enable future comparisons between the results from the two instruments. The convolution was possible thanks to the broad extent of GPI's H band, whose wavelength window covers both SPHERE narrow-band H filters.

In this way, it was possible to place every PACO candidate in a (H2-H3, H2) CMD (Figure 5). We used confirmed background objects from the SHINE survey – that offers a larger sample statistics thanks to the wide 11"x11" field of view of IRDIS (Dohlen et al. 2008) – to build an "exclusion zone", defined as the region of the CMD that encompasses all the points within $5\sigma$ from the mean colors of background sources as a function of their absolute magnitude. As in SHINE publications, the exclusion zone was set to begin at $H2 = 16$ mag, as the existence of some planets (e.g., HR 8799 b) with $H2 \sim 15$ mag and $H2-H3 \sim 0$ mag renders the method unreliable at brighter magnitudes (Langlois et al. 2021, Chomez et al. 2024). We labeled as "companion candidates" the sources lying along the T track or having additional indications hinting towards a bound nature, and as "ambiguous" the sources in the region $H2 < 16$ mag and $H2-H3 \sim 0$.

Excluding already known substellar companions, all but nine sources can be confidently ruled out as background contaminants. Seven of these are classified as ambiguous according to

our vetting scheme, and will not be hereafter discussed. The nature of the two remaining promising candidates – whose photometry and age is consistent with 5-8 $M_{Jup}$ objects – is currently unclear. The candidate around HIP 78663 is located in a position of the CMD where the colors of bound companions overlay those of background stars; however, we classify it as a promising candidate because of a tentative $\sim 3.5\sigma$ detection in the shallower second epoch possibly hinting at common proper motion. As regards the candidate around HD 24072, in addition to the hypothesis of a bound nature, the following scenarios might be envisaged to explain its position along the young-object track:

1. a free-floating planet or brown dwarf, belonging to the same association as the target and hence possessing similar colors to substellar companion while not exhibiting a large variation of the distance modulus;

2. a statistical false positive (see Section 3.2). A spectral dependence of the photocenter of a false positive might happen to mimic, during the characterization step of PACO, blue spectra similar to those of real substellar companions lying along the T track.

Spurious detections in direct imaging have previously arisen due to extended objects (proto-planetary and debris disks) that were poorly subtracted (see, e.g., Sallum et al. 2015 and confusion by Currie et al. 2019), but we exclude this possibility given the lack of infrared excess in WISE (Wright et al. 2010) bands.

We finally notice that the HD 24072 system also comprises a low-mass star, closer to the primary than the planet candidate (see Section D); under the assumption of face-on circular orbits, we empirically verified, based on the results by Musielak et al. (2005), that the candidate would be far enough from the substellar companion to be dynamically stable.

Our reanalysis redetected all substellar companions (7 planets, 3 brown dwarfs) that we expected to find on the basis of the literature (Figure 6). Some of these companions – notably, HR 8799 c, d, and e – have just one epoch in our observing sample; consistently with the decision tree described above, we would have been able to confirm them as bound objects through proper motion test, employing additional available SPHERE or GPI epochs.

We report in Table 1 details about the astrometry and the photometry of these 10 substellar companions. Given the extensive characterization of these objects already undertaken in the literature, we deem a re-derivation of masses and semi-major axes – the main input needed for the statistical analysis of Section 3.4 – to be outside of the scope of this paper; instead, we decided to recover the most accurate values from dedicated literature works.

In addition to substellar companions, we detected 6 sources whose high luminosity points towards a stellar nature. We were able to confirm 5 of them as physically bound thanks to 1) a proper motion strongly disagreeing with background stars, and 2) astrometric wobbles indicated by the Gaia astrometric solution of the corresponding primaries; the remaining binary candidate to HD 74341B, with no second epoch and too bright to employ the CMD test, awaits confirmation. Details are provided in Section D.

## 3.2. False positives 587

As mentioned in Section 3.1, the roster of PACO candidates does not include a few detections identified as false positives, induced either 1) by real astrophysical or optical features, or 2) by statistical fluctuations of the SNR map. The former category includes residuals of the first Airy ring around very bright sources
<figure><img src="image_7.png" /><figcaption>Fig. 6: PACO or cADI detection maps for the substellar companions detected in the survey (indicated by arrows). PACO maps are to be read as SNR maps, sharing a common colorbar. Individual colorbars are shown below the two cADI maps.</figcaption></figure>

<table><caption>Table 1: Detected substellar companions to stars in the sample.</caption><thead><tr><th>Name</th><th>subsample</th><th>Date</th><th>SNR</th><th>sep<br>arcsec</th><th>PA<br>deg</th><th>ΔH<br>mag</th><th>H<br>mag</th><th>H2 - H3<br>mag</th><th>mass<br>M<sub>sup</sub></th><th>sma<br>au</th><th>source</th></tr></thead><tbody><tr><td>HD 206893 b</td><td>FGK</td><td>2016-09-22</td><td>7.2</td><td>0.268 ± 0.002</td><td>61.5 ± 0.4</td><td>11.1</td><td>13.7 ± 0.1</td><td>0.5 ± 0.1</td><td>28</td><td>9.6</td><td>1,1</td></tr><tr><td>HR 8799 c</td><td>BA</td><td>2016-09-19</td><td>37.1</td><td>0.953 ± 0.002</td><td>330.8 ± 0.2</td><td>12.5</td><td>14.8 ± 0.2</td><td>0.4 ± 0.2</td><td>8.5</td><td>41</td><td>2,2</td></tr><tr><td>HR 8799 d</td><td>BA</td><td>2016-09-19</td><td>19.4</td><td>0.666 ± 0.002</td><td>223.5 ± 0.2</td><td>12.5</td><td>14.7 ± 0.2</td><td>0.3 ± 0.2</td><td>8.5</td><td>27</td><td>2,2</td></tr><tr><td>HR 8799 e</td><td>BA</td><td>2016-09-19</td><td>9.2</td><td>0.393 ± 0.003</td><td>284.1 ± 0.4</td><td>12.2</td><td>14.4 ± 0.2</td><td>0.3 ± 0.2</td><td>9.6</td><td>16</td><td>3,2</td></tr><tr><td rowspan="5">51 Eri b</td><td rowspan="5">BA</td><td>2014-12-18</td><td>14.7</td><td>0.439 ± 0.002</td><td>171.3 ± 0.3</td><td>14.3</td><td>16.7 ± 0.2</td><td>-1.2 ± 0.1</td><td rowspan="5">4.1</td><td rowspan="5">11.1</td><td rowspan="5">4,5</td></tr><tr><td>2015-01-31</td><td>5.0</td><td>0.456 ± 0.007</td><td>170.4 ± 0.8</td><td>15.1</td><td>—</td><td>—</td></tr><tr><td>2015-09-01</td><td>10.9</td><td>0.442 ± 0.003</td><td>166.8 ± 0.4</td><td>14.7</td><td>17.1 ± 0.2</td><td>-1.0 ± 0.2</td></tr><tr><td>2015-12-20</td><td>7.5</td><td>0.443 ± 0.005</td><td>166.5 ± 0.7</td><td>14.4</td><td>16.8 ± 0.2</td><td>-0.6 ± 0.2</td></tr><tr><td>2016-09-18</td><td>12.7</td><td>0.442 ± 0.003</td><td>162.0 ± 0.4</td><td>14.4</td><td>16.8 ± 0.1</td><td>-0.8 ± 0.1</td></tr><tr><td rowspan="3">β Pic b</td><td rowspan="3">BA</td><td>2015-11-06</td><td>8.3</td><td>0.421 ± 0.004</td><td>359.1 ± 0.6</td><td>9.8</td><td>11.9 ± 0.2</td><td>0.4 ± 0.1</td><td rowspan="3">11.9</td><td rowspan="3">9.93</td><td rowspan="3">6,6</td></tr><tr><td>2015-12-22</td><td>7.7</td><td>0.241 ± 0.002</td><td>213.2 ± 0.4</td><td>9.6</td><td>11.7 ± 0.2</td><td>0.4 ± 0.1</td></tr><tr><td>2016-01-21</td><td>7.5</td><td>0.226 ± 0.002</td><td>212.3 ± 0.6</td><td>10.4</td><td>12.5 ± 0.3</td><td>0.4 ± 0.3</td></tr><tr><td rowspan="4">HD 95086 b</td><td rowspan="4">BA</td><td>2013-12-11</td><td>8.9</td><td>0.635 ± 0.003</td><td>150.6 ± 0.3</td><td>13.8</td><td>15.9 ± 0.2</td><td>0.2 ± 0.2</td><td rowspan="4">2.6</td><td rowspan="4">61.7</td><td rowspan="4">7,8</td></tr><tr><td>2016-02-29</td><td>9.3</td><td>0.627 ± 0.003</td><td>148.1 ± 0.4</td><td>13.6</td><td>15.8 ± 0.2</td><td>0.2 ± 0.1</td></tr><tr><td>2016-03-06</td><td>8.1</td><td>0.626 ± 0.004</td><td>148.1 ± 0.4</td><td>14.0</td><td>16.1 ± 0.2</td><td>0.7 ± 0.2</td></tr><tr><td>2016-03-28</td><td>4.9</td><td>0.629 ± 0.006</td><td>147.5 ± 0.6</td><td>14.2</td><td>16.4 ± 0.2</td><td>0.0 ± 0.3</td></tr><tr><td rowspan="4">HR 2562 b</td><td rowspan="4">FGK</td><td>2016-01-25</td><td>17.0</td><td>0.605 ± 0.002</td><td>294.7 ± 0.3</td><td>12.8</td><td>15.3 ± 0.2</td><td>0.6 ± 0.3</td><td rowspan="4">10.28</td><td rowspan="4">21.2</td><td rowspan="4">9,9</td></tr><tr><td>2017-02-13</td><td>21.4</td><td>0.635 ± 0.002</td><td>298.3 ± 0.2</td><td>11.3</td><td>13.8 ± 0.1</td><td>0.3 ± 0.1</td></tr><tr><td>2017-11-29</td><td>12.5</td><td>0.654 ± 0.002</td><td>297.8 ± 0.3</td><td>12.8</td><td>15.2 ± 0.2</td><td>0.4 ± 0.1</td></tr><tr><td>2018-11-19</td><td>28.1</td><td>0.677 ± 0.002</td><td>297.3 ± 0.2</td><td>11.3</td><td>13.7 ± 0.3</td><td>0.3 ± 0.3</td></tr><tr><td>HD 984 B<sup>a</sup></td><td>FGK</td><td>2015-08-30</td><td>—</td><td>0.219 ± 0.002</td><td>84.0 ± 0.4</td><td>7.7</td><td>—</td><td>—</td><td>61</td><td>28</td><td>10,10</td></tr><tr><td>PZ Tel B<sup>a</sup></td><td>FGK</td><td>2015-07-30</td><td>—</td><td>0.501 ± 0.002</td><td>59.5 ± 0.2</td><td>5.9</td><td>—</td><td>—</td><td>27</td><td>27</td><td>11,11</td></tr></tbody></table>

Notes. $^a$: reduction through the custom cADI. (mass, sma) sources for planet properties: 1: Hinkley et al. (2023); 2: Zurlo et al. (2022); 3: Brandt et al. (2021); 4: Elliott et al. (2024); 5: De Rosa et al. (2020b); 6: Lacour et al. (2021); 7: Nielsen et al. (2019); 8: Rameau et al. (2016); 9: Zhang et al. (2023); 10: Franson et al. (2022); 11: Franson & Bowler (2023). The column named "subsample" indicates whether the parent star belongs to the FGK or the BA subsample (see Sec. 3.4).

<footnote>(2 cases) and disk residuals (10 cases); the latter (8 cases) is constituted by unusually bright residuals that had no counterpart in</footnote>

<footnote>additional GPI or SPHERE observations with better or similar detection limits.</footnote>
With respect to the latter case, we tried to estimate the number of false positives expected to arise from statistical fluctuations. We recall that the distribution of pixel intensities in PACO SNR maps is a normalized Gaussian $N(0, 1)$. In this case, a $5\sigma$ threshold corresponds to a false alarm probability $p_{5\sigma} = 2.9 \cdot 10^{-7}$. Given the number of pixels in GPI's FOV, $N_{px} = 185^2$ and the number of effectively independent spectral priors, $N_p$, that we empirically estimate as $N_p \approx 4^{12}$, using a binomial distribution, we expect ~ 20 false positives across the entire survey (see Chomez et al. 2023). This number is larger than the number of statistical false positives that could be identified through second-epoch observations; hence, we expect that some sources labeled as CMD background sources also belong to the category.

### 3.3. Completeness

The completeness of our survey was quantified in the following way. As a first step, we azimuthally averaged the 2D detection maps provided by PACO, obtaining 1D contrast curves (detection limits at $5\sigma$).

Pending a final confirmation of the nature of the two promising candidate companions, the detection limits of the corresponding observations were adjusted accordingly to ensure the statistic reliability of the corresponding observations. The same was done for the seven datasets containing ambiguous sources. In particular, the mean contrast of each candidate was employed as a floor value in the corresponding 2D $5\sigma$ map. In other words, we pretend to have shallower observations so that a source as bright as the candidate can be at most a marginal $5\sigma$ detection. These maps were then collapsed to 1D as described above.

The 1D curves obtained in this way were converted into mass limits through MADYS, adopting the stellar parameters discussed in Section 2.3. The mass-luminosity relation is based on the ATMO evolutionary models (hereafter ATMO; Phillips et al. 2020; Chabrier et al. 2023)$^{13}$. Chemical disequilibrium is expected to critically affect the atmospheric features of cool T-type and Y-type objects (see, e.g. Leggett et al. 2015, 2017; Miles et al. 2020; Baxter et al. 2021); given that 1) the corresponding temperature range is within the reach of our analysis (see Figure 5), and 2) the effect is particularly strong in H-band observations such as those under consideration, we decided to employ the grid assuming weak chemical disequilibrium (ATMO-NEQ-W) instead of chemical equilibrium (ATMO-CEQ). We explored in Appendix E the effect of this assumption, comparing the results with those obtained under chemical equilibrium and strong disequilibrium (ATMO-NEQ-S)$^{14}$. In addition to this, the impact of model selection and age uncertainty was quantified.

Starting from mass limits, the completeness could be estimated through Exo-DMC$^{15}$ (Bonavita 2020). Within each cell of a 2D grid in the (mass, sma) plane, the detectability of $N = 1000$ companions to every star, whose orbital parameters are drawn in a Monte Carlo fashion, is computed by comparison with $5\sigma$

<figure><img src="image_8.png" /><figcaption>Fig. 7: Survey completeness as a function of companion mass and semi-major axis, computed using the ATMO-NEQ-W models. Red stars indicate known substellar companions (see Table 1).</figcaption></figure>

<table><caption>Table 2: Input parameters used for Exo-DMC.</caption><thead><tr><td>parameter</td><td>description</td></tr></thead><tbody><tr><td>no. of steps in sma</td><td>500</td></tr><tr><td>no. of steps in mass</td><td>200</td></tr><tr><td>no. of draws per cell</td><td>1000</td></tr><tr><td>semi-major axis</td><td>log-uniform in [0.1, 1000] au</td></tr><tr><td>companion mass</td><td>log-uniform in [0.1, 100] M<sub>Jup</sub></td></tr><tr><td>inclination</td><td>cos i uniform in [-1, 1]</td></tr><tr><td>eccentricity</td><td>|N(0, 0.9)|</td></tr><tr><td>longitude of node</td><td>uniform in [0, 2π]</td></tr><tr><td>longitude of periastron</td><td>uniform in [0, 2π]</td></tr><tr><td>fraction of period</td><td>uniform in [0, 1]</td></tr></tbody></table>

mass limits as a function of the projected separation. We provide in Table 2 an overview of the adopted parameters.

The final map is computed by taking the average of all the individual maps. When multiple epochs for a given star are available, the largest value for the detectability is selected for every cell. The results are shown in Figure 7. We notice that the peak sensitivity of the survey is about 88%: we interpret such a low value as the combination of three factors: 1) the small field of view of the instrument; 2) the moderate distance spread across the sample; 3) the fact that, working in semi-major axis and not in projected separation, a fraction of planets with given $a$ might be undetectable because of projection effects.

We are now able to directly compare our detection capabilities with those of N19, so as to justify a posteriori the idea of a reanalysis of those archival data. In order to avoid any possible systematic difference, a new map was computed only using the observations considered therein; moreover, instead of using Exo-DMC, we closely reproduced the original method, including the same values for distances, ages and substellar evolutionary model. The comparison, shown in Figure 8, indicates that the PACO-based reanalysis allows for a significant performance gain at all separations, which can be up to twofold in terms of detectable mass at given completeness.

<footnote>$^{12}$ The correlation between SNR maps under any two spectral priors is larger than zero. By "number of effectively independent spectral priors" we mean the ratio $N_{positives,20priors}/N_{positives,1prior}$, estimated through extensive testing.</footnote>

<footnote>$^{13}$ Using the most recent version, that features a new equation of state for dense hydrogen-helium mixtures: https://noctis.erc-atmo.eu/fsdownload/zyU96xA6o/phillips2020.</footnote>

<footnote>$^{14}$ The amount of vertical mixing in disequilibrium models is parametrized through the eddy diffusion coefficient $K_{zz}$. Constraining $K_{zz}$ is a long-standing issue (see discussion in Phillips et al. 2020).</footnote>

<footnote>$^{15}$ https://github.com/mbonav/Exo_DMC</footnote>
<figure><img src="image_9.png" /><figcaption>Fig. 8: Comparison of the survey completeness between N19 (dashed lines) and this work (solid lines). Only the observations used by N19 were used to draw this plot.</figcaption></figure>

## 3.4. Planet occurrence rates

Deriving unbiased occurrence frequencies of exoplanets is one of the main goals of large blind surveys, and in turn a crucial input to draw comparisons with formation models. Provided a large enough sample, it is additionally possible to investigate the dependence of these frequencies on host properties such as mass and metallicity, highlighting the key role of the parent star in shaping its planetary system.

We begin our investigation by focusing on the occurrence frequency $f$ for the entire stellar population represented by the GPIES sample. Extracting this quantity from the fact of having observed $N$ companions given a certain survey completeness is a typical inversion problem that can be treated within a Bayesian framework.

We employ a formalism that is similar to that used in previous direct imaging studies (see, e.g. Lafrenière et al. 2007; Lannier et al. 2016). Given a certain area $A$ in the (sma, mass) plane defined by $a_{min} < a < a_{max}$ and $m_{min} < m < m_{max}$, let us define $p_i$ the mean probability to see a companion around the $i$-th star lying within $A$. Based on our completeness analysis (Section 3.3), $p_i$ can be estimated as the mean detection probability in $A$ across the entire survey, that is the mean value in $A$ of the completeness map shown in Figure 7.

The probability $p_{det,i}$ to detect a companion in $A$ around the $i$-th star is the product of the detection probability and the underlying occurrence frequency $f$: $p_{det,i} = p_i \cdot f$. The connection with the observed planet sample is mediated by $d$, a vector whose $i$-element represents the number of companions detected within $A$ around the $i$-th star.

The likelihood of the observed data as a function of the $f$ can be estimated as the product of individual Bernoulli events, one per star:

$$L(\{d_i\}|f) = \prod_{i=1}^{N} (1 - p_{det,i})^{1-d_i} \times (p_{det,i})^{d_i} \quad (2)$$

The probability density function of $f$, that is the occurrence frequency of companions in $A$ given the data can be finally estimated through Bayes' theorem:

$$P(f|\{d_i\}) = \frac{L(\{d_i\}|f)P(f)}{\int_0^1 L(\{d_i\}|f)P(f)df} \quad (3)$$

as the posterior distribution emerging from the interplay between a suitable prior distribution $P(f)$ and the likelihood $L(\{d_i\}|f)$. We adopt two distinct priors: a uniform prior and a Jeffreys prior. The uniform prior:

$$P(f) \propto 1, \qquad \forall f \in [0, 1] \quad (4)$$

despite not incorporating any observational information, is not uninformative, as it assumes much larger weights for large values of $f$ compared to what is expected from observations. Nevertheless, the simplicity of this prior makes it widely adopted in the literature: we decided to employ it in order to allow for comparison with published results.

A Jeffreys prior has the twofold advantage of being non-informative and counterbalancing the bias that favors $f \sim 0.5$. In the case of Bernoulli events, the Jeffreys prior for the parameter $f$ is simply:

$$P(f) = \frac{1}{\sqrt{f \cdot (1-f)}} \quad (5)$$

We adopt the latter prior distribution, that has the advantage of being non-informative, as our standard choice in the following analysis.

A particularly delicate point is represented by the choice of $A$: on the one hand, selecting a too narrow range would result in a critical amplification of fluctuations from small number statistics; on the other hand, including regions where $\langle\{p_i\}\rangle \sim 0$ would require a significant amount of extrapolation due to the lack of data and, consequently, induce a flattening of the posterior distribution over the prior. An additional factor to take into account is the dependence of the results on both age uncertainty and model selection, becoming more severe as the lower mass limit is decreased (Appendix E). We decided to consider, as our nominal case, a lower mass limit of $5M_{Jup}$ and a semi-major axis range $10 \text{ au} < a < 100 \text{ au}$ as a compromise between these concurrent factors; the upper mass limit will be set to either $13 M_{Jup}$ or $80 M_{Jup}$ depending on whether brown dwarf companions are considered or not. We derive occurrence frequencies of $1.7^{+0.9}_{-0.7}\%$ when $A = [5, 13] M_{Jup} \times [10, 100]$ au, and $2.2^{+1.0}_{-0.8}\%$ when $A = [5, 80] M_{Jup} \times [10, 100]$ au, with $f$ represented by the median of the posterior and the error bars defined from the [16%, 84%] percentiles. In order to allow a straightforward comparison with previous results from the literature (Section 4), we also present additional occurrences starting from different definitions of $A$ (Table 3). We notice that no significant difference is derived from the prior choice, confirming that our careful choice for $A$ did minimize the impact of the prior; the errorbars, as expected, are smaller in Jeffreys case compared to the uniform case. In addition to this, no significant deviation arises for this
<table><caption>Table 3: Occurrence rates for different definitions of the (mass, sma) range $\mathcal{A}$ and for the two choices for the prior distribution (U: uniform; J: Jeffreys). $^a$: 95% upper limit.</caption><thead><tr><th>$\mathcal{A}$<br>$M_{\text{Jup}} \times \text{au}$</th><th>SpT</th><th>$f_U$<br>%</th><th>$f_J$<br>%</th></tr></thead><tbody><tr><td>[5, 13] × [10, 100]</td><td>all</td><td>1.9<sup>+1.0</sup><sub>-0.7</sub></td><td>1.7<sup>+0.9</sup><sub>-0.7</sub></td></tr><tr><td>[5, 13] × [10, 100]</td><td>BA</td><td>4.3<sup>+2.6</sup><sub>-1.9</sub></td><td>3.8<sup>+2.4</sup><sub>-1.7</sub></td></tr><tr><td>[5, 13] × [10, 100]</td><td>FGK</td><td>1.0<sup>+1.0</sup><sub>-0.6</sub></td><td>0.7<sup>+0.9</sup><sub>-0.5</sub></td></tr><tr><td>[5, 80] × [10, 100]</td><td>all</td><td>2.4<sup>+1.0</sup><sub>-0.8</sub></td><td>2.2<sup>+1.0</sup><sub>-0.8</sub></td></tr><tr><td>[5, 80] × [10, 100]</td><td>BA</td><td>3.5<sup>+2.1</sup><sub>-1.5</sub></td><td>3.0<sup>+2.0</sup><sub>-1.4</sub></td></tr><tr><td>[5, 80] × [10, 100]</td><td>FGK</td><td>2.2<sup>+1.3</sup><sub>-0.9</sub></td><td>1.9<sup>+1.2</sup><sub>-0.9</sub></td></tr><tr><td>[2, 13] × [10, 100]</td><td>all</td><td>3.5<sup>+1.5</sup><sub>-1.2</sub></td><td>3.2<sup>+1.5</sup><sub>-1.1</sub></td></tr><tr><td>[2, 13] × [3, 100]</td><td>all</td><td>5.3<sup>+2.1</sup><sub>-1.7</sub></td><td>5.0<sup>+2.1</sup><sub>-1.6</sub></td></tr><tr><td>[2, 13] × [5, 300]</td><td>all</td><td>5.3<sup>+2.1</sup><sub>-1.7</sub></td><td>5.0<sup>+2.1</sup><sub>-1.6</sub></td></tr><tr><td>[13, 80] × [5, 100]</td><td>all</td><td>1.4<sup>+0.9</sup><sub>-0.6</sub></td><td>1.2<sup>+0.8</sup><sub>-0.6</sub></td></tr><tr><td>[13, 80] × [5, 100]</td><td>BA</td><td>&lt; 3.8<sup>a</sup></td><td>&lt; 2.7<sup>a</sup></td></tr><tr><td>[13, 80] × [5, 100]</td><td>FGK</td><td>2.3<sup>+1.4</sup><sub>-1.0</sub></td><td>2.0<sup>+1.3</sup><sub>-0.9</sub></td></tr></tbody></table>

choice of $\mathcal{A}$ as an effect of the theoretical assumptions and observational uncertainties, ensuring the robustness of our results (Appendix E).

We consider host star metallicity not to be a factor of particular concern, as the metallicity of young star-forming regions in the solar neighborhood is typically solar with limited spread (D’Orazi et al. 2011; Biazzo et al. 2012; Baratella et al. 2020; Magrini et al. 2023). Conversely, as done in Nielsen et al. (2019) and Vigan et al. (2021), we explicitly investigate the dependence of the occurrence frequency on stellar mass. We divided our sample in three bins of stellar masses, obtaining the BA subsample ($M > 1.5M_\odot$, 160 stars), the FGK subsample ($0.5 < M \le 1.5M_\odot$, 235 stars), and the M subsample ($M \le 0.5M_\odot$, 5 stars). Given its small size, the M star sample was discarded.

Both the aggregated results and the mass-dependent ones are plotted in Fig. 9. Occurrence frequencies for different values of $\mathcal{A}$ are provided for reference in Table 3; moreover, a digitized version of completeness maps is also made available<fnref n="16" /> so as to allow interested readers to extract additional results based on different definitions of $\mathcal{A}$.

## 4. Discussion

Thanks to the analysis performed in Section 3.4, it is now possible to directly compare the results emerging from our PACO re-reduction of GPI data to previous literature works.

Figure 10 presents a juxtaposition of our results with the frequencies derived from the first 300 stars of GPIES (Nielsen et al. 2019) and the first 150 stars of SHINE (SHINE F150; Vigan et al. 2021). Moreover, the results emerging from the meta-analysis of 384 imaged stars by Bowler (2016) are shown, although we stress that they are, by design, less protected against selection biases due to the heterogeneous underlying sample. A general finding is that our results are fully compatible with literature estimates, but they are typically more precise. In particular, the comparison with N19 clearly indicates that the new analysis places the tightest constraints to date on giant planet occurrence based on GPI data, with a gain in precision being a direct consequence of the large gain in completeness (Figure 8). The smaller

<figure><img src="image_10.png" /><figcaption>Fig. 9: Occurrence frequency of GP (upper panel) and GP+BD (lower panel) from the re-analysis of the 400-star GPIES sample presented in this work. Aggregated results are shown in blue, whereas results for the BA and the FGK subsample are plotted in green and orange, respectively. The colored area encompasses the [16th, 84th] of the posterior distribution.</figcaption></figure>

frequency of companions is an effect of the increased completeness with no new confirmed detection. As regards SHINE F150, the much larger field of view of IRDIS (11"x11", compared to the 2.7"x2.7" FOV of GPI) ensures a much more complete coverage of the semi-major axis range of interest and thus larger room for planet detection, compensating the twofold advantage of our study in terms of sample size and reduction performances: as a result, the precision of the derived occurrence rates is similar. In this regard, the full analysis of SHINE data with PACO (Chomez et al. 2024), which combines all the advantages of the two analyses, is expected to provide an invaluable contribution to demographic studies of wide-orbit exoplanets.

In view of the profound consequences with respect to planet formation scenarios, it is extremely interesting to assess the dependence of the observed occurrence rates on stellar mass. As in N19, we employed a threshold value of $M = 1.5M_\odot$ to distinguish a BA subsample and a FGK subsample. That study claimed that a significant (3.4$\sigma$) tension between BA and FGK planet rates ($f_{\text{BA}}$ and $f_{\text{FGK}}$, respectively) exists for $\mathcal{A} = [2, 13] M_{\text{Jup}} \times [3, 100]$ au $:= \mathcal{A}_1$, with planets being more

786

787

788

789

790

791

792

793

794

795

796

797

798

799

800

801

802

803

804

805

<footnote n="16">¹⁶ url of Zenodo repository, to be created upon acceptance</footnote>
<figure><img src="image_11.png" /><figcaption>Fig. 10: Comparison of the occurrence rates (left panel: $\mathcal{A} = [2, 13] M_{\text{Jup}} \times [3, 100]$ au; right panel: $\mathcal{A} = [5, 13] M_{\text{Jup}} \times [10, 100]$ au) of giant planets with previous analyses by Bowler (2016), Nielsen et al. (2019) and Vigan et al. (2021). The left and the right half of each panel are relative to BA and FGK stars, respectively. Estimates indicated by arrows are to be read as 95% upper limits, while error bars on point estimates are defined as to encompass the 68% C.I. Question marks indicate missing data points.</figcaption></figure>

common around BA hosts; while six companions were detected inside $\mathcal{A}_1$ in the BA subsample, no companion was identified in the FGK subsample. However, we expect the result to be weak-ened by the re-revaluation of the mass of HR 2562 b (Zhang et al. 2023), a companion to an FGK star now firmly placed into the planetary-mass domain. In order to verify whether this is the case, we drew, in a Monte Carlo fashion, values from the poste-rior distributions of $f_{\text{BA}}$ and $f_{\text{FGK}}$ under $\mathcal{A}_1$. The values from the latter distribution are larger than those drawn from the former in 0.1% of the cases, implying a 3.3 $\sigma$ tension between the two distributions. Hence, our analysis confirms the finding by N19.

With respect to brown dwarf companions, no statistically sig-nificant difference in the observed rates was found by N19 be-tween BA and FGK hosts. The observation is in line with the results of previous analyses showing compatible rates across a wide range of stella types (see, e.g., Nielsen et al. 2013; Bowler et al. 2015; Lannier et al. 2016; Bowler & Nielsen 2018). Based on our analysis, a tentative (1.7 $\sigma$) tension between the two rates is found for $\mathcal{A} = [13, 80] M_{\text{Jup}} \times [5, 100]$ au, with an interest-ing inversion compared to the planetary case: in other words, giant planets appear to be more common around BA host, while brown dwarf companions tentatively appear to be more common around FGK hosts.

Although the BD trend is not statistically significant, an in-teresting analogy might be drawn with the behavior of the two empirical distributions of substellar companions introduced by Vigan et al. (2021) in the context of SHINE. A planet-like and a star-like distribution of companions – both being the product of a log-normal distribution for semimajor axis and a power-law for companion-to-star mass ratios – were simultaneously fitted to the substellar companion population, divided in three bins of mass (BA, FGK, M). Similarly to our planetary rates, the me-dian values of the planet-like posterior are larger than those of the star-like posterior for BA hosts, and smaller for FGK hosts. It might be argued that a strict distinction between giant planets and brown dwarfs based on the deuterium burning limit is not adequate to capture the complexity of the different formation mechanisms (CA and GI) involved (Chabrier et al. 2014), and

that studying together the entire population is the key to identify population trends (see, e.g., Gratton et al. 2024); while this is cer-tainly true, first-order, population-wise differences in some pa-rameter, arguing for different underlying formation channels, can sometimes be discerned using rough mass boundaries (Bowler et al. 2020).

In view of the low occurrence rates, the large extent of host star masses, the interplay of different formation channels and the impact of input assumptions, we deem it necessary to defer a thorough study of the distribution of companion properties to our future joint SPHERE+GPI analysis: thanks to its larger sample size, this sample is expected to bring about much tighter con-straints on the properties of the companion population, offering in turn the possibility to compare them both to empirical distri-butions and to synthetic populations of companions produced by formation models.

## 5. Conclusions

We have presented in this work the results of a complete re-reduction of 400 stars from the GPIES survey, one of the largest planet-hunting DI endeavors to date, by means of an advanced post-processing algorithm named PACO. The key results of this work are the following:

- the detection capabilities of the survey were greatly en-hanced by means of our novel post-processing technique, reaching up to a twofold gain in terms of detectable mass at given completeness;

- out of 102 detected sources, 2 were identified as promising companion candidates awaiting follow-up confirmation;

- thanks to the deeper detection limits provided by PACO, it was possible to place some of the deepest constraints ever provided by direct imaging on the occurrence of wide-orbit giant planets. We derive an occurrence rate of $1.7^{+1.0}_{-0.7}\%$ for $5 M_{\text{Jup}} < m < 13 M_{\text{Jup}}$ planets in $10 \text{ au} < a < 100 \text{ au}$, increasing to $2.3^{+1.0}_{-0.8}\%$ when including substellar companions up to $80 M_{\text{Jup}}$;
879 – we verified that the above-mentioned results are robust against the effect of age uncertainty, model selection, and disequilibrium chemistry;

880 – as in previous studies, we observe (3.3 σ C.L.) a larger occurrence rate of giant planets around BA hosts compared to FGK stars;

885 – we tentatively (1.7 σ C.L.) identify an inversion of this trend when considering brown dwarf companions, with FGK stars possibly hosting more such companion than their BA counterparts.

In a forthcoming study, we plan to combine the archives of SPHERE and GPI data, leading to a threefold sample size compared to this work. By applying the same reduction and analysis methods presented here, it will be possible to assess a whole series of stimulating questions related to the origin, the prevalence and the properties of wide-orbit planets. In addition to this, these endeavors will enable a decisive step towards the coveted combination of demographic constraints derived through different detection techniques, delivering in turn key inputs for planet formation models suited to a wide variety of host stars.

*Acknowledgements.* This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (COBREX; grant agreement n° 885593).

SPHERE is an instrument designed and built by a consortium consisting of IPAG (Grenoble, France), MPIA (Heidelberg, Germany), LAM (Marseille, France),

LESLIA (Paris, France), Laboratoire Lagrange (Nice, France), INAF - Observatorio 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/LESLIA (Paris), and Observatoire de Lyon/CRAL,

and is supported by a grant from Labex OSUG@2020 (Investissements d'avenir - ANR10 LABX56).

This work is based on observations obtained at the Gemini Observatory, which is

operated by the Association of Universities for Research in Astronomy, Inc., under a cooperative agreement with the NSF on behalf of the Gemini partnership:

the National Science Foundation (United States), the National Research Council (Canada), CONICYT (Chile), the Australian Research Council (Australia),

Ministério Ciência, Tecnologia e Inovação (Brazil) and Ministerio de Ciencia,

Tecnología e Innovación Productiva (Argentina).

This research has made use of data obtained from or tools provided by the portal exoplanet.eu of The Extrasolar Planets Encyclopaedia.

This research has made use of the VizieR catalogue access tool, CDS, Strasbourg, France (DOI: 10.26093/cds/vizier). The original description of the VizieR service was published in 2000, A&AS 143, 23.

This research has made use of the SIMBAD database, CDS, Strasbourg Astronomical Observatory, France. The original description of the SIMBAD database was published in 2000, A&AS 143, 9.

This work has made use of data from the SHINE GTO survey, operated at SPHERE@VLT.

This research has made use of the SVO Filter Profile Service "Carlos Rodriguez", funded by MCIN/AEI/10.13039/501100011033/through grant PID2020-112949GB-I00.

We thank Schuyler Grace Wolff for her help to calculate GPI calibration files using GPI DRP.

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.

Software: numpy (Harris et al. 2020), astropy (Astropy Collaboration et al. 2013), astroquery (Ginsburg et al. 2019), madys (Squicciarini & Bonavita 2022), GaiaPMEX (Kiefer et al. 2024), pyklip (Wang et al. 2015b), Exo-DMC (Bonavita 2020), PACO (Flasseur et al. 2018).

## References

948

Allard, F., Homeier, D., & Freytag, B. 2011, in Astronomical Society of the Pacific Conference Series, Vol. 448, 16th Cambridge Workshop on Cool Stars, Stellar Systems, and the Sun, ed. C. Johns-Krull, M. K. Browning, & A. A. West, 91

Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33

Baraffe, I., Chabrier, G., Barman, T. S., Allard, F., & Hauschildt, P. H. 2003, A&A, 402, 701

Baratella, M., D'Orazi, V., Carraro, G., et al. 2020, A&A, 634, A34

Baxter, C., Désert, J.-M., Tsai, S.-M., et al. 2021, A&A, 648, A127

Bell, C. P. M., Mamajek, E. E., & Naylor, T. 2015, MNRAS, 454, 593

Beuzit, J. L., Vigan, A., Mouillet, D., et al. 2019, A&A, 631, A155

Biazzo, K., Bozza, V., Mancini, L., & Sozzetti, A. 2022, in Astrophysics and Space Science Library, Vol. 466, Demographics of Exoplanetary Systems, Lecture Notes of the 3rd Advanced School on Exoplanetary Science

Biazzo, K., D'Orazi, V., Desidera, S., et al. 2012, MNRAS, 427, 2905

Boccaletti, A., Chauvin, G., Wildi, F., et al. 2022, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 12184, Ground-based and Airborne Instrumentation for Astronomy IX, ed. C. J. Evans, J. J. Bryant, & K. Motohara, 121841S

Bonavita, M. 2020, Exo-DMC: Exoplanet Detection Map Calculator, Astrophysics Source Code Library, record ascl:2010.008

Bonavita, M., Gratton, R., Desidera, S., et al. 2022, A&A, 663, A144

Bonnefoy, M., Perraut, K., Lagrange, A. M., et al. 2018, A&A, 618, A63

Boss, A. P. 1997, Science, 276, 1836

Bourrier, V., Dumusque, X., Dorn, C., et al. 2018, A&A, 619, A1

Bowler, B. P. 2016, PASP, 128, 102001

Bowler, B. P., Blunt, S. C., & Nielsen, E. L. 2020, AJ, 159, 63

Bowler, B. P., Liu, M. C., Shkolnik, E. L., & Tamura, M. 2015, ApJS, 216, 7

Bowler, B. P. & Nielsen, E. L. 2018, in Handbook of Exoplanets, ed. H. J. Deeg & J. A. Belmonte, 155

Brandt, G. M., Brandt, T. D., Dupuy, T. J., Michalik, D., & Marleau, G.-D. 2021, ApJ, 915, L16

Bryan, M. L., Knutson, H. A., Howard, A. W., et al. 2016, ApJ, 821, 89

Chabrier, G., Baraffe, I., Phillips, M., & Debras, F. 2023, A&A, 671, A119

Chabrier, G., Johansen, A., Janson, M., & Rafikov, R. 2014, in Protostars and Planets VI, ed. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning, 619-642

Charnay, B., Bézard, B., & Baudino, J.-L. 2019, in EPSC-DPS Joint Meeting 2019, Vol. 2019, EPSC-DPS2019-1450

Chauvin, G., Desidera, S., Lagrange, A. M., et al. 2017, in SF2A-2017: Proceedings of the Annual meeting of the French Society of Astronomy and Astrophysics, ed. C. Reylé, P. Di Matteo, F. Herpin, E. Lagadec, A. Lançon, Z. Meliani, & F. Royer, Di

Chauvin, G., Lagrange, A. M., Dumas, C., et al. 2004, A&A, 425, L29

Chilcote, J., Konopacky, Q., De Rosa, R. J., et al. 2020, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 11447, Ground-based and Airborne Instrumentation for Astronomy VIII, ed. C. J. Evans, J. J. Bryant, & K. Motohara, 114471S

Chilcote, J. K., Bailey, V. P., De Rosa, R., et al. 2018, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 10702, Ground-based and Airborne Instrumentation for Astronomy VII, ed. C. J. Evans, L. Simard, & H. Takami, 1070244

Chomez, A., Delorme, P., & Lagrange, A.-M. 2024, A&A, submitted

Chomez, A., Lagrange, A. M., Delorme, P., et al. 2023, A&A, 675, A205

Courtney-Barrer, B., De Rosa, R., Kokotanekova, R., et al. 2023, A&A, 680, A34

Covino, E., Esposito, M., Barbieri, M., et al. 2013, A&A, 554, A28

Cumming, A., Butler, R. P., Marcy, G. W., et al. 2008, PASP, 120, 531

Currie, T., Biller, B., Lagrange, A., et al. 2023, in Astronomical Society of the Pacific Conference Series, Vol. 534, Protostars and Planets VII, ed. S. Inutsuka, Y. Aikawa, T. Muto, K. Tomida, & M. Tamura, 799

Currie, T., Marois, C., Cieza, L., et al. 2019, ApJ, 877, L3

De Rosa, R. J., Nguyen, M. M., Chilcote, J., et al. 2020a, Journal of Astronomical Telescopes, Instruments, and Systems, 6, 015006

De Rosa, R. J., Nielsen, E. L., Wahhaj, Z., et al. 2023, A&A, 672, A94

De Rosa, R. J., Nielsen, E. L., Wang, J. J., et al. 2020b, AJ, 159, 1

de Zeeuw, P. T., Hoogerwerf, R., de Bruijne, J. H. J., Brown, A. G. A., & Blaauw, A. 1999, AJ, 117, 354

Delorme, P., Meunier, N., Albert, D., et al. 2017, in SF2A-2017: Proceedings of the Annual meeting of the French Society of Astronomy and Astrophysics, ed. C. Reylé, P. Di Matteo, F. Herpin, E. Lagadec, A. Lançon, Z. Meliani, & F. Royer, Di

Desidera, S., Chauvin, G., Bonavita, M., et al. 2021, A&A, 651, A70

Desidera, S., Squicciarini, V., & Brandner, W. 2024, A&A, in prep.

949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023

948

949

950

951

952

953

954

955

956

957

958

959

960

961

962

963

964

965

966

967

968

969

970

971

972

973

974

975

976

977

978

979

980

981

982

983

984

985

986

987

988

989

990

991

992

993

994

995

996

997

998

999

1000

1001

1002

1003

1004

1005

1006

1007

1008

1009

1010

1011

1012

1013

1014

1015

1016

1017

1018

1019

1020

1021

1022

1023
<table>
  <tbody>
    <tr>
      <td>1024</td>
      <td>Dohlen, K., Langlois, M., Saisse, M., et al. 2008, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 7014, Ground-based and Airborne Instrumentation for Astronomy II, ed. I. S. McLean & M. M. Casali, 70143L</td>
    </tr>
    <tr>
      <td>1025</td>
      <td>D'Orazi, V., Biazzo, K., & Randich, S. 2011, A&A, 526, A103</td>
    </tr>
    <tr>
      <td>1029</td>
      <td>Draper, Z. H., Marois, C., Wolff, S., et al. 2014, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 9147, Ground-based and Airborne Instrumentation for Astronomy V, ed. S. K. Ramsay, I. S. McLean, & H. Takami, 91474Z</td>
    </tr>
    <tr>
      <td>1033</td>
      <td>Dulz, S. D., Plavchan, P., Crepp, J. R., et al. 2020, ApJ, 893, 122</td>
    </tr>
    <tr>
      <td>1034</td>
      <td>Elliott, A., Boyajian, T., Ellis, T., et al. 2024, PASA, 41, e043</td>
    </tr>
    <tr>
      <td>1035</td>
      <td>Flasseur, O., Denis, L., Thiébaut, É., & Langlois, M. 2018, A&A, 618, A138</td>
    </tr>
    <tr>
      <td>1036</td>
      <td>Flasseur, O., Denis, L., Thiébaut, É., & Langlois, M. 2020a, A&A, 637, A9</td>
    </tr>
    <tr>
      <td>1037</td>
      <td>Flasseur, O., Denis, L., Thiébaut, É., & Langlois, M. 2020b, A&A, 634, A2</td>
    </tr>
    <tr>
      <td>1038</td>
      <td>Franson, K. & Bowler, B. P. 2023, AJ, 165, 246</td>
    </tr>
    <tr>
      <td>1039</td>
      <td>Franson, K., Bowler, B. P., Brandt, T. D., et al. 2022, AJ, 163, 50</td>
    </tr>
    <tr>
      <td>1040</td>
      <td>Franson, K., Bowler, B. P., Zhou, Y., et al. 2023, ApJ, 950, L19</td>
    </tr>
    <tr>
      <td>1041</td>
      <td>Gagné, J., Faherty, J. K., & Mamajek, E. E. 2018a, ApJ, 865, 136</td>
    </tr>
    <tr>
      <td>1042</td>
      <td>Gagné, J., Mamajek, E. E., Malo, L., et al. 2018b, ApJ, 856, 23</td>
    </tr>
    <tr>
      <td>1043</td>
      <td>Gaia Collaboration, Vallenari, A., Brown, A. G. A., et al. 2023, A&A, 674, A1</td>
    </tr>
    <tr>
      <td>1044</td>
      <td>Galicher, R., Marois, C., Macintosh, B., et al. 2016, A&A, 594, A63</td>
    </tr>
    <tr>
      <td>1045</td>
      <td>Galicher, R. & Mazoyer, J. 2024, Comptes Rendus Physique, 24, 133</td>
    </tr>
    <tr>
      <td>1046</td>
      <td>Gandolfi, D., Barragán, O., Hatzes, A. P., et al. 2017, AJ, 154, 123</td>
    </tr>
    <tr>
      <td>1047</td>
      <td>Gaudi, B. S., Meyer, M., & Christiansen, J. 2021, in ExoFrontiers; Big Questions in Exoplanetary Science, ed. N. Madhusudhan (IOP Publishing Ltd), 2–1</td>
    </tr>
    <tr>
      <td>1049</td>
      <td>Ginsburg, A., Sipőcz, B. M., Brasseur, C. E., et al. 2019, AJ, 157, 98</td>
    </tr>
    <tr>
      <td>1050</td>
      <td>Gratton, R., Bonavita, M., Mesa, D., et al. 2024, A&A, 685, A119</td>
    </tr>
    <tr>
      <td>1051</td>
      <td>Gratton, R., Mesa, D., Bonavita, M., et al. 2023, Nature Communications, 14, 6232</td>
    </tr>
    <tr>
      <td>1052</td>
      <td>Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357</td>
    </tr>
    <tr>
      <td>1053</td>
      <td>Hinkley, S., Lacour, S., Marleau, G. D., et al. 2023, A&A, 671, L5</td>
    </tr>
    <tr>
      <td>1054</td>
      <td>Janson, M., Squicciarini, V., Delorme, P., et al. 2021, A&A, 646, A164</td>
    </tr>
    <tr>
      <td>1055</td>
      <td>Kasdin, N. J., Bailey, V. P., Mennesson, B., et al. 2020, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 11443, Space Telescopes and Instrumentation 2020: Optical, Infrared, and Millimeter Wave, ed. M. Lystrup & M. D. Perrin, 114431U</td>
    </tr>
    <tr>
      <td>1056</td>
      <td>Kasper, M., Cerpa Urra, N., Pathak, P., et al. 2021, The Messenger, 182, 38</td>
    </tr>
    <tr>
      <td>1057</td>
      <td>Keppler, M., Benisty, M., Müller, A., et al. 2018, A&A, 617, A44</td>
    </tr>
    <tr>
      <td>1058</td>
      <td>Kervella, P., Arenou, F., Mignard, F., & Thévenin, F. 2019, A&A, 623, A72</td>
    </tr>
    <tr>
      <td>1059</td>
      <td>Kervella, P., Arenou, F., & Thévenin, F. 2022, A&A, 657, A7</td>
    </tr>
    <tr>
      <td>1060</td>
      <td>Kiefer, F., Lagrange, A.-M., Rubini, P., & Philipot, F. 2024, arXiv e-prints, arXiv:2409.16992</td>
    </tr>
    <tr>
      <td>1061</td>
      <td>Konopacky, Q. M., Rameau, J., Duchêne, G., et al. 2016, ApJ, 829, L4</td>
    </tr>
    <tr>
      <td>1062</td>
      <td>Kuzuhara, M., Currie, T., Takarada, T., et al. 2022, ApJ, 934, L18</td>
    </tr>
    <tr>
      <td>1063</td>
      <td>Lacedelli, G., Malavolta, L., Borsato, L., et al. 2021, MNRAS, 501, 4148</td>
    </tr>
    <tr>
      <td>1064</td>
      <td>Lacour, S., Wang, J. J., Rodet, L., et al. 2021, A&A, 654, L2</td>
    </tr>
    <tr>
      <td>1065</td>
      <td>Lafrenière, D., Doyon, R., Marois, C., et al. 2007, ApJ, 670, 1367</td>
    </tr>
    <tr>
      <td>1066</td>
      <td>Lagrange, A. M., Gratadour, D., Chauvin, G., et al. 2009, A&A, 493, L21</td>
    </tr>
    <tr>
      <td>1067</td>
      <td>Lagrange, A. M., Philipot, F., Rubini, P., et al. 2023, A&A, 677, A71</td>
    </tr>
    <tr>
      <td>1068</td>
      <td>Langlois, M., Gratton, R., Lagrange, A. M., et al. 2021, A&A, 651, A71</td>
    </tr>
    <tr>
      <td>1069</td>
      <td>Lanner, J., Delorme, P., Lagrange, A. M., et al. 2016, A&A, 596, A83</td>
    </tr>
    <tr>
      <td>1070</td>
      <td>Launhardt, R., Henning, T., Quirrenbach, A., et al. 2020, A&A, 635, A162</td>
    </tr>
    <tr>
      <td>1071</td>
      <td>Leggett, S. K., Morley, C. V., Marley, M. S., & Saumon, D. 2015, ApJ, 799, 37</td>
    </tr>
    <tr>
      <td>1072</td>
      <td>Leggett, S. K., Tremblin, P., Esplin, T. L., Luhman, K. L., & Morley, C. V. 2017, ApJ, 842, 118</td>
    </tr>
    <tr>
      <td>1073</td>
      <td>Leike, R. H., Glatzle, M., & Enßlin, T. A. 2020, A&A, 639, A138</td>
    </tr>
    <tr>
      <td>1074</td>
      <td>Lindegren, L., Klioner, S. A., Hernández, J., et al. 2021, A&A, 649, A2</td>
    </tr>
    <tr>
      <td>1075</td>
      <td>Liu, M. C., Wahhaj, Z., Biller, B. A., et al. 2010, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 7736, Adaptive Optics Systems II, ed. B. L. Ellerbroek, M. Hart, N. Hubin, & P. L. Wizinowich, 77361K</td>
    </tr>
    <tr>
      <td>1076</td>
      <td>Macintosh, B., Graham, J. R., Barman, T., et al. 2015, Science, 350, 64</td>
    </tr>
    <tr>
      <td>1077</td>
      <td>Macintosh, B., Graham, J. R., Ingraham, P., et al. 2014, Proceedings of the National Academy of Science, 111, 12661</td>
    </tr>
    <tr>
      <td>1078</td>
      <td>Macintosh, B. A., Anthony, A., Atwood, J., et al. 2014, in Proceedings of the SPIE, Vol. 9148 (SPIE), 152–165</td>
    </tr>
    <tr>
      <td>1079</td>
      <td>Magrini, L., Viscasillas Vázquez, C., Spina, L., et al. 2023, A&A, 669, A119</td>
    </tr>
    <tr>
      <td>1080</td>
      <td>Maire, J., Ingraham, P. J., De Rosa, R. J., et al. 2014, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 9147, Ground-based and Airborne Instrumentation for Astronomy V, ed. S. K. Ramsay, I. S. McLean, & H. Takami, 914785</td>
    </tr>
    <tr>
      <td>1081</td>
      <td>Marley, M. S., Saumon, D., Visscher, C., et al. 2021, ApJ, 920, 85</td>
    </tr>
    <tr>
      <td>1082</td>
      <td>Marois, C., Lafrenière, D., Doyon, R., Macintosh, B., & Nadeau, D. 2006, ApJ, 641, 556</td>
    </tr>
    <tr>
      <td>1083</td>
      <td>Marois, C., Macintosh, B., Barman, T., et al. 2008, Science, 322, 1348</td>
    </tr>
    <tr>
      <td>1084</td>
      <td>McBride, J., Graham, J. R., Macintosh, B., et al. 2011, PASP, 123, 692</td>
    </tr>
    <tr>
      <td>1085</td>
      <td>Mesa, D., Gratton, R., Kervella, P., et al. 2023, A&A, 672, A93</td>
    </tr>
    <tr>
      <td>1086</td>
      <td>Miles, B. E., Skemer, A. J. I., Morley, C. V., et al. 2020, AJ, 160, 63</td>
    </tr>
    <tr>
      <td>1087</td>
      <td>Mordasini, C., Alibert, Y., & Benz, W. 2009, A&A, 501, 1139</td>
    </tr>
    <tr>
      <td>1088</td>
      <td>Murphy, S. J., Lawson, W. A., & Bessell, M. S. 2013, MNRAS, 435, 1325</td>
    </tr>
  </tbody>
</table>

<table>
  <tbody>
    <tr>
      <td>1104</td>
      <td>Musielak, Z. E., Cuntz, M., Marshall, E. A., & Stuit, T. D. 2005, A&A, 434, 355</td>
    </tr>
    <tr>
      <td>1105</td>
      <td>Nguyen, C. T., Costa, G., Girardi, L., et al. 2022, A&A, 665, A126</td>
    </tr>
    <tr>
      <td>1106</td>
      <td>Nielsen, E. L., De Rosa, R. J., Macintosh, B., et al. 2019, AJ, 158, 13</td>
    </tr>
    <tr>
      <td>1107</td>
      <td>Nielsen, E. L., Liu, M. C., Wahhaj, Z., et al. 2013, ApJ, 776, 4</td>
    </tr>
    <tr>
      <td>1108</td>
      <td>Ochsenbein, F., Bauer, P., & Marcout, J. 2000, A&AS, 143, 23</td>
    </tr>
    <tr>
      <td>1109</td>
      <td>Pecaut, M. J. & Mamajek, E. E. 2016, MNRAS, 461, 794</td>
    </tr>
    <tr>
      <td>1110</td>
      <td>Perrin, M. D., Ingraham, P., Follette, K. B., et al. 2016, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 9908, Ground-based and Airborne Instrumentation for Astronomy VI, ed. C. J. Evans, L. Simard, & H. Takami, 990837</td>
    </tr>
    <tr>
      <td>1111</td>
      <td>Perrin, M. D., Maire, J., Ingraham, P., et al. 2014, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 9147, Ground-based and Airborne Instrumentation for Astronomy V, ed. S. K. Ramsay, I. S. McLean, & H. Takami, 91473J</td>
    </tr>
    <tr>
      <td>1112</td>
      <td>Philipot, F., Lagrange, A. M., Rubini, P., Kiefer, F., & Chomez, A. 2023, A&A, 670, A65</td>
    </tr>
    <tr>
      <td>1113</td>
      <td>Phillips, M. W., Tremblin, P., Baraffe, I., et al. 2020, A&A, 637, A38</td>
    </tr>
    <tr>
      <td>1114</td>
      <td>Pollack, J. B., Hubickyj, O., Bodenheimer, P., et al. 1996, Icarus, 124, 62</td>
    </tr>
    <tr>
      <td>1115</td>
      <td>Poyneer, L. A., Palmer, D. W. Macintosh, B., et al. 2016, Appl. Opt., 55, 323</td>
    </tr>
    <tr>
      <td>1116</td>
      <td>Rameau, J., Chauvin, G., Lagrange, A. M., et al. 2013, A&A, 553, A60</td>
    </tr>
    <tr>
      <td>1117</td>
      <td>Rameau, J., Nielsen, E. L., De Rosa, R. J., et al. 2016, ApJ, 822, L29</td>
    </tr>
    <tr>
      <td>1118</td>
      <td>Reggiani, M., Meyer, M. R., Chauvin, G., et al. 2016, A&A, 586, A147</td>
    </tr>
    <tr>
      <td>1119</td>
      <td>Rogers, L. A. 2015, ApJ, 801, 41</td>
    </tr>
    <tr>
      <td>1120</td>
      <td>Rosenthal, L. J., Howard, A. W., Knutson, H. A., & Fulton, B. J. 2024, ApJS, 270, 1</td>
    </tr>
    <tr>
      <td>1121</td>
      <td>Ruffio, J.-B., Macintosh, B., Wang, J. J., et al. 2017, ApJ, 842, 14</td>
    </tr>
    <tr>
      <td>1122</td>
      <td>Sallum, S., Follette, K. B., Eisner, J. A., et al. 2015, Nature, 527, 342</td>
    </tr>
    <tr>
      <td>1123</td>
      <td>Santerne, A., Moutou, C., Tsantaki, M., et al. 2016, A&A, 587, A64</td>
    </tr>
    <tr>
      <td>1124</td>
      <td>Skrutskie, M. F., Cutri, R. M., Stiening, R., et al. 2006, AJ, 131, 1163</td>
    </tr>
    <tr>
      <td>1125</td>
      <td>Spiegel, D. S. & Burrows, A. 2012, ApJ, 745, 174</td>
    </tr>
    <tr>
      <td>1126</td>
      <td>Squicciarini, V. & Bonavita, M. 2022, A&A, 666, A15</td>
    </tr>
    <tr>
      <td>1127</td>
      <td>Stone, J. M., Skemer, A. J., Hinz, P. M., et al. 2018, AJ, 156, 286</td>
    </tr>
    <tr>
      <td>1128</td>
      <td>Tran, D., Konopacky, Q., & GPIES Collaboration. 2016, in American Astronomical Society Meeting Abstracts, Vol. 228, American Astronomical Society Meeting Abstracts #228, 120.01</td>
    </tr>
    <tr>
      <td>1129</td>
      <td>Uyama, T., Hashimoto, J., Kuzuhara, M., et al. 2017, AJ, 153, 106</td>
    </tr>
    <tr>
      <td>1130</td>
      <td>Vigan, A., Fontanive, C., Meyer, M., et al. 2021, A&A, 651, A72</td>
    </tr>
    <tr>
      <td>1131</td>
      <td>Vorobyov, E. I. 2013, A&A, 552, A129</td>
    </tr>
    <tr>
      <td>1132</td>
      <td>Waisberg, I., Klein, Y., & Katz, B. 2023a, Research Notes of the American Astronomical Society, 7, 125</td>
    </tr>
    <tr>
      <td>1133</td>
      <td>Waisberg, I., Klein, Y., & Katz, B. 2023b, Research Notes of the American Astronomical Society, 7, 78</td>
    </tr>
    <tr>
      <td>1134</td>
      <td>Wang, J. J., Perrin, M. D., Savransky, D., et al. 2018, Journal of Astronomical Telescopes, Instruments, and Systems, 4, 018002</td>
    </tr>
    <tr>
      <td>1135</td>
      <td>Wang, J. J., Rajan, A., Graham, J. R., et al. 2014, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 9147, Ground-based and Airborne Instrumentation for Astronomy V, ed. S. K. Ramsay, I. S. McLean, & H. Takami, 914755</td>
    </tr>
    <tr>
      <td>1136</td>
      <td>Wang, J. J., Ruffio, J.-B., De Rosa, R. J., et al. 2015a, pyKLIP: PSF Subtraction for Exoplanets and Disks, Astrophysics Source Code Library, record ascl:1506.001</td>
    </tr>
    <tr>
      <td>1137</td>
      <td>Wang, J. J., Ruffio, J.-B., De Rosa, R. J., et al. 2015b, pyKLIP: PSF Subtraction for Exoplanets and Disks, Astrophysics Source Code Library, record ascl:1506.001</td>
    </tr>
    <tr>
      <td>1138</td>
      <td>Wolff, S. G., Perrin, M. D., Maire, J., et al. 2014, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 9147, Ground-based and Airborne Instrumentation for Astronomy V, ed. S. K. Ramsay, I. S. McLean, & H. Takami, 91477H</td>
    </tr>
    <tr>
      <td>1139</td>
      <td>Wolthoff, V., Reffert, S., Quirrenbach, A., et al. 2022, A&A, 661, A63</td>
    </tr>
    <tr>
      <td>1140</td>
      <td>Wright, E. L., Eisenhardt, P. R. M., Mainzer, A. K., et al. 2010, AJ, 140, 1868</td>
    </tr>
    <tr>
      <td>1141</td>
      <td>Zhang, S. Y., Duchêne, G., De Rosa, R. J., et al. 2023, AJ, 165, 219</td>
    </tr>
    <tr>
      <td>1142</td>
      <td>Zhu, W. & Wu, Y. 2018, AJ, 156, 92</td>
    </tr>
    <tr>
      <td>1143</td>
      <td>Zuckerman, B. 2019, ApJ, 870, 27</td>
    </tr>
    <tr>
      <td>1144</td>
      <td>Zuckerman, B., Bessell, M. S., Song, I., & Kim, S. 2006, ApJ, 649, L115</td>
    </tr>
    <tr>
      <td>1145</td>
      <td>Zuckerman, B., Rhee, J. H., Song, I., & Bessell, M. S. 2011, ApJ, 732, 61</td>
    </tr>
    <tr>
      <td>1146</td>
      <td>Zuckerman, B., Song, I., Bessell, M. S., & Webb, R. A. 2001, ApJ, 562, L87</td>
    </tr>
    <tr>
      <td>1147</td>
      <td>Zurlo, A. 2024, arXiv e-prints, arXiv:2404.05797</td>
    </tr>
    <tr>
      <td>1148</td>
      <td>Zurlo, A., Goździewski, K., Lazzoni, C., et al. 2022, A&A, 666, A133</td>
    </tr>
    <tr>
      <td>1149</td>
      <td>Zurlo, A., Goździewski, K., Lazzoni, C., et al. 2022, A&A, 666, A133</td>
    </tr>
  </tbody>
</table>
## Appendix A: The stellar sample

<table><caption>Table A.1: Stellar properties for the sample considered in this work. The full table is available in electronic form at the CDS.</caption><thead><tr><th>star name</th><th>ra<sup>a</sup></th><th>dec<sup>a</sup></th><th>parallax<sup>a</sup></th><th>SpT<sup>b</sup></th><th>Gmag<sup>a</sup></th><th>Hmag<sup>c</sup></th><th>E(B-V)</th><th>YMG</th><th>p<sub>memb</sub></th><th>age</th><th>age ref.<sup>d</sup></th><th>mass</th></tr><tr><th></th><th>hms</th><th>dms</th><th>mas</th><th></th><th>mag</th><th>mag</th><th>mag</th><th></th><th></th><th>Myr</th><th></th><th>M<sub>☉</sub></th></tr></thead><tbody><tr><td>HD 104467</td><td>12 01 39.1168</td><td>-78 59 16.915</td><td>10.18 ± 0.12</td><td>G3V(e)</td><td>8.44</td><td>6.97</td><td>0.02</td><td>EPSC</td><td>1.00</td><td>3.7<sup>+4.6</sup><sub>-1.4</sub></td><td>B</td><td>1.74<sup>+0.09</sup><sub>-0.24</sub></td></tr><tr><td>HD 105874A</td><td>12 11 14.8135</td><td>-52 13 03.187</td><td>8.07 ± 0.99</td><td></td><td>7.88</td><td></td><td>0.04</td><td>LCC</td><td>0.99</td><td>15 ± 3</td><td>B</td><td>1.69 ± 0.08</td></tr><tr><td>HD 118991A</td><td>13 41 44.7704</td><td>-54 33 33.934</td><td>11.28 ± 0.10</td><td>B8.5Vn</td><td>5.24</td><td>5.45</td><td>0.02</td><td>FIELD</td><td>0.72</td><td>15.0 ± 3.8</td><td>N</td><td>3.32 ± 0.17</td></tr><tr><td>HD 129926B</td><td>14 46 00.5907</td><td>-25 26 39.973</td><td>32.62 ± 0.02</td><td>G1V</td><td>6.95</td><td>5.72</td><td>0.01</td><td>FIELD</td><td>0.71</td><td>500 ± 120</td><td>N</td><td>1.08 ± 0.05</td></tr><tr><td>HD 131399A</td><td>14 54 25.3089</td><td>-34 08 34.038</td><td>10.20 ± 0.70</td><td>A1V</td><td>7.07</td><td></td><td>0.04</td><td>UCL</td><td>0.99</td><td>16 ± 2</td><td>B</td><td>1.94 ± 0.10</td></tr><tr><td>HD 137919A</td><td>15 30 21.31</td><td>-41 55 08.33</td><td>7.93 ± 0.72</td><td></td><td></td><td>6.46</td><td>0.03</td><td>UCL</td><td>0.99</td><td>16 ± 2</td><td>B</td><td>3.54 ± 0.18</td></tr><tr><td>HD 141943</td><td>15 53 27.2916</td><td>-42 16 00.71</td><td>16.63 ± 0.02</td><td>G2</td><td>7.79</td><td>6.41</td><td>0.01</td><td>FIELD</td><td>0.48</td><td>16 ± 4</td><td>N</td><td>1.22<sup>+0.09</sup><sub>-0.06</sub></td></tr><tr><td>HD 147553A</td><td>16 23 56.7146</td><td>-33 11 57.828</td><td>7.23 ± 0.04</td><td>B9.5V(n)</td><td>7.00</td><td>7.01</td><td>0.04</td><td>UCL</td><td>0.95</td><td>16 ± 2</td><td>B</td><td>2.51 ± 0.13</td></tr><tr><td>HD 16699A</td><td>2 38 44.2802</td><td>-52 57 03.053</td><td>17.27 ± 0.02</td><td>F8V</td><td>7.75</td><td>6.70</td><td>0.00</td><td>ARG</td><td>0.98</td><td>45 ± 5</td><td>B</td><td>1.22 ± 0.06</td></tr><tr><td>HD 16699B</td><td>2 38 45.0461</td><td>-52 57 08.451</td><td>16.69 ± 0.21</td><td>G8V</td><td>8.24</td><td>6.63</td><td>0.00</td><td>ARG</td><td>0.95</td><td>45 ± 5</td><td>B</td><td>1.15 ± 0.06</td></tr></tbody></table>

Notes. Data taken from: <sup>a</sup>: Gaia DR3; <sup>b</sup>: Simbad; <sup>c</sup>: 2MASS; <sup>d</sup>: (B)ANYAN, (N)ielsen+19, (S)HINE. Details about the derivation of E(B-V), YMG membership, ages and masses are provided in Section 2.3.

<table><caption>Table A.2: Adopted ages for the YMG of interest.</caption><thead><tr><th>YMG</th><th>Age</th><th>Source</th></tr></thead><tbody><tr><td>AB Doradus</td><td>149<sup>+31</sup><sub>-49</sub></td><td>1</td></tr><tr><td>Argus</td><td>45 ± 5</td><td>2</td></tr><tr><td>β Pic MG</td><td>24 ± 5</td><td>1</td></tr><tr><td>Carina</td><td>45<sup>+11</sup><sub>-7</sub></td><td>3</td></tr><tr><td>Carina-Near</td><td>200 ± 50</td><td>4</td></tr><tr><td>Columba</td><td>42<sup>+6</sup><sub>-4</sub></td><td>3</td></tr><tr><td>ε Cha</td><td>3.7<sup>+4.6</sup><sub>-1.4</sub></td><td>5</td></tr><tr><td>Lower Centaurus-Crux</td><td>15 ± 3</td><td>6</td></tr><tr><td>Tucana Horologium Ass.</td><td>45<sup>+5</sup><sub>-4</sub></td><td>3</td></tr><tr><td>TW Hya Ass.</td><td>10 ± 3</td><td>3</td></tr><tr><td>Upper Centaurus-Lupus</td><td>16 ± 2</td><td>6</td></tr><tr><td>Upper Scorpius</td><td>10 ± 3</td><td>6</td></tr><tr><td>Volans-Carina Ass.</td><td>87<sup>+5</sup><sub>-9</sub></td><td>7</td></tr></tbody></table>

Notes. Sources: 1: Desidera et al. (2021); 2: Zuckerman (2019); 3: Bell et al. (2015); 4: Zuckerman et al. (2006), assuming a relative 25% error; 5: Murphy et al. (2013); 6: Pecaut & Mamajek (2016); 7: Gagné et al. (2018a).
## 1175 Appendix B: Observation logs

<table><caption>Table B.1: Observing log for the observations considered in this work. The full table is available in electronic form at the CDS.</caption><thead><tr><th rowspan="2">star name</th><th rowspan="2">GPIES name</th><th rowspan="2">obs. night</th><th colspan="2">τ<sub>0</sub><sup>a</sup></th><th rowspan="2">seeing<sup>b</sup></th><th rowspan="2">airmass</th><th rowspan="2">int. time<sup>c</sup></th><th rowspan="2">ΔPA</th><th rowspan="2">program</th></tr><tr><th>ms</th><th>arcsec</th><th>s</th><th>deg</th></tr></thead><tbody><tr><td>HIP 2472</td><td>HIP2472</td><td>2013-11-13</td><td></td><td>0.680</td><td>1.077</td><td>30x1x49.46</td><td>15.7</td><td>GS-ENG-GPI-COM</td></tr><tr><td>HIP 53524</td><td>HD95086</td><td>2013-12-10</td><td></td><td>0.330</td><td>1.356</td><td>21x1x119.29</td><td>15.0</td><td>GS-ENG-GPI-COM</td></tr><tr><td>HIP 64995</td><td>HD 115600</td><td>2014-04-22</td><td></td><td>0.365</td><td>1.157</td><td>58x1x49.46</td><td>32.5</td><td>GS-2014A-SV-403</td></tr><tr><td>HIP 11964</td><td>CC Eri</td><td>2014-11-08</td><td></td><td>0.865</td><td>1.057</td><td>36x1x59.65</td><td>31.9</td><td>GS-2014B-Q-500</td></tr><tr><td>HIP 12964</td><td>HR 826</td><td>2014-11-08</td><td></td><td>0.930</td><td>1.015</td><td>36x1x59.65</td><td>69.6</td><td>GS-2014B-Q-500</td></tr><tr><td>HIP 560</td><td>HR 9</td><td>2014-11-08</td><td></td><td>0.890</td><td>1.031</td><td>35x1x59.65</td><td>44.1</td><td>GS-2014B-Q-500</td></tr><tr><td>HIP 19893</td><td>gam Dor</td><td>2014-11-09</td><td></td><td>0.635</td><td>1.073</td><td>38x1x59.65</td><td>28.7</td><td>GS-2014B-Q-500</td></tr><tr><td>HIP 12413</td><td>HR 789</td><td>2014-11-09</td><td></td><td>0.785</td><td>1.039</td><td>40x1x59.65</td><td>43.5</td><td>GS-2014B-Q-500</td></tr><tr><td>HIP 490</td><td>HD 105</td><td>2014-11-09</td><td></td><td>0.615</td><td>1.026</td><td>40x1x59.65</td><td>31.2</td><td>GS-2014B-Q-500</td></tr><tr><td>HIP 25283</td><td>HD 35650</td><td>2014-11-09</td><td></td><td>0.650</td><td>1.023</td><td>40x1x59.65</td><td>24.9</td><td>GS-2014B-Q-500</td></tr></tbody></table>

Notes. $^a$: coming from MASS measurements, not available before April 2015 and after April 2017. $^b$: average between MASS and DIMM measurements; stale MASS measurements (non-zero values repeated over long-daily to monthly-periods of time; see Poyneer et al. 2016) were identified and removed. No MASS/DIMM values were available after April 2017. $^c$: int. time = number of frames × number of co-added images × Detector Integration Time per frame. ΔPA represents the parallactic rotation over the sequence.

1176

1177

## 1178 Appendix C: Companion candidates

<table><caption>Table C.1: Companion candidates detected in this work. The table is available in electronic form at the CDS.</caption><thead><tr><th rowspan="2">star name</th><th rowspan="2">obs. night</th><th rowspan="2">SNR</th><th colspan="2">arcsec</th><th colspan="2">contrast</th><th colspan="2">H2</th><th rowspan="2">algo.</th><th rowspan="2">classification</th></tr><tr><th>deg</th><th>mag</th><th>mag</th><th>mag</th></tr></thead><tbody><tr><td>HD 104467</td><td>2018-03-26</td><td>5.8</td><td>0.368 ± 0.005</td><td>7.5 ± 0.9</td><td>13.3</td><td>15.4 ± 0.1</td><td>0.3 ± 0.2</td><td>PACO</td><td>ambiguous</td></tr><tr><td>HD 118991A</td><td>2015-04-04</td><td>45.9</td><td>1.162 ± 0.002</td><td>217.7 ± 0.2</td><td>10.8</td><td>11.6 ± 0.1</td><td>0.3 ± 0.2</td><td>PACO</td><td>pm bkg</td></tr><tr><td>HD 131399A</td><td>2017-02-15</td><td>9.6</td><td>0.801 ± 0.003</td><td>193.9 ± 0.3</td><td>13.5</td><td>15.4 ± 0.1</td><td>0.4 ± 0.2</td><td>PACO</td><td>pm bkg</td></tr><tr><td>HD 131399A</td><td>2017-04-20</td><td>10.8</td><td>0.802 ± 0.003</td><td>194.0 ± 0.3</td><td>13.2</td><td>15.1 ± 0.1</td><td>0.4 ± 0.3</td><td>PACO</td><td>pm bkg</td></tr><tr><td>HD 24072</td><td>2018-11-22</td><td>—</td><td>0.193 ± 0.002</td><td>16.5 ± 0.4</td><td>6.9</td><td>—</td><td>—</td><td>cADI</td><td>star comp</td></tr><tr><td>HD 24072</td><td>2018-11-22</td><td>5.8</td><td>0.466 ± 0.004</td><td>337.1 ± 0.5</td><td>14.1</td><td>14.3 ± 0.7</td><td>-2.0<sup>+0.1</sup><sub>-0.3</sub></td><td>PACO</td><td>interesting</td></tr><tr><td>HD 36869</td><td>2016-12-17</td><td>17.2</td><td>0.783 ± 0.002</td><td>212.6 ± 0.2</td><td>12.7</td><td>15.9 ± 0.1</td><td>0.1 ± 0.1</td><td>PACO</td><td>pm bkg</td></tr><tr><td>HD 36869</td><td>2016-12-17</td><td>15.0</td><td>0.432 ± 0.002</td><td>107.5 ± 0.3</td><td>12.1</td><td>15.3 ± 0.1</td><td>0.0 ± 0.1</td><td>PACO</td><td>pm bkg</td></tr><tr><td>HD 74341B</td><td>2015-12-20</td><td>10.4</td><td>1.454 ± 0.007</td><td>342.3 ± 0.4</td><td>13.8</td><td>16.3 ± 0.1</td><td>0.0 ± 1.2</td><td>PACO</td><td>cmd bkg</td></tr><tr><td>HD 84330B</td><td>2015-12-18</td><td>20.7</td><td>1.123 ± 0.003</td><td>247.5 ± 0.3</td><td>12.8</td><td>15.7 ± 0.1</td><td>0.0 ± 0.1</td><td>PACO</td><td>pm bkg</td></tr><tr><td>HD 84330B</td><td>2016-03-18</td><td>11.8</td><td>1.072 ± 0.003</td><td>246.0 ± 0.3</td><td>13.1</td><td>16.0 ± 0.1</td><td>0.1 ± 0.4</td><td>PACO</td><td>pm bkg</td></tr></tbody></table>

1179

1180

Notes. H2: absolute SPHERE H2 magnitude. H2-H3: SPHERE H2-H3 color. Classification: cmd bkg = background star via CMD; pm bkg: background star via proper motion analysis; sub comp: substellar companion; star comp: stellar companion. $^a$: unconfirmed, see Section D.
1181

## Appendix D: Stellar companions

1182

1183

1184

1185

1186

1187

1188

1189

1190

1191

1192

1193

1194

1195

1196

1197

1198

1199

1200

1201

1202

1203

1204

1205

1206

1207

1208

1209

1210

1211

1212

1213

1214

1215

1216

1217

1218

1219

1220

1221

1222

1223

1224

1225

1226

1227

1228

1229

1230

1231

1232

1233

1234

1235

1236

1237

1238

1239

1240

discrepancy to the fact that these sources lie at the edge of the coronagraph, where the transmission is much lower than elsewhere across the field of view.

1241

1242

1243

In addition to the brown dwarfs HD 984 B and PZ Tel B, the ADI reduction identified 7 bright companion candidates (two of them detected twice in two different epochs). Proper motion analysis allowed us to identify one of them (namely, the one seen next to HIP 61087) as a background object and one as a bound companion (around HIP 74696). For the remaining objects, for which only one observation was available in our sample, we searched for archival detections in the literature. It turns out that all the candidates but one (around HD 74341B) had been already imaged in the course of past campaigns, but just one (around HIP 38160) had already been confirmed as a comoving object through follow-up observations (Rameau et al. 2013). Therefore, we performed the proper motion analysis for all the systems, using the astrometric measurements reported in Table D.1.

The proper motion test confirmed that the 5 sources with two epochs exhibit a significantly different motion compared to static background objects, with large displacements related to orbital motion (Figure D.1).

In order to clarify the status of HD 74341 B, to further corroborate the bound nature of the other objects and, finally, to test the reliability of the derived photometric masses, we ran the GaiaPMEX tool (Kiefer et al. 2024) to see if the astrometry of the primary from Gaia and/or Hipparcos showed hints of wobbles indicative of the presence of an unseen companion. GaiaPMEX comes equipped with a model of the Renormalized Unit Weight Error (ruwe; see Lindegren et al. 2021) and of the Gaia-Hipparcos proper motion anomaly (PMa; see Kervella et al. 2019, 2022) distribution expected for a single star as a function of stellar magnitude and colors. The evaluation of whether astrometric information is consistent with an unseen companion is performed in the following way. After defining a log-uniform grid of companion masses $M_c \in [0.1, 3000] M_{Jup}$ and semi-major axes $a \in [0.01, 1000]$ au with $30 \times 30$ bins, the program draws, within each bin, 100 ($\log M_c$, $\log a$)-doublets from a uniform distribution. As regards the other orbital parameters, they are randomly extracted from the distributions described in Table D.2. We employ stellar parallaxes from Gaia DR3, while stellar masses are recovered from our analysis described in Section 2.3.

At each node of the mass-sma grid, distributions of the ruwe and/or PMa are determined given the target and its hypothetical orbiting companion, and compared to the actual ruwe and/or PMa; the derivation of confidence regions for possible companion masses and semi-major axes can be finally obtained through Bayesian inversion.

Clear astrometric detections were found for all the targets but HD 74341B (Figure D.2), due to the absence of the star in the Hipparcos catalog; due to the much shorter timespan of the astrometric measurements underlying the ruwe (~ 3 yr, compared with the ~ 24 yr of the PMa), the sensitivity of the ruwe at the relatively large separation of the companion candidate is virtually null.

We find good agreement between the photometric and the dynamical masses for three stellar companions (those around HIP 74696, HD 24072, and HIP 26369). As regards the companions to HIP 67199 and HIP 38160, which are the ones located at the shortest separations from the star, we confirm their bound nature but we find largely underestimated masses, despite the expedient to use no-ADI instead of cADI<fnref n="17" />. We attribute the

<footnote n="17">¹⁷ The no-ADI algorithm can be thought as a cADI but with no median-subtraction step. The advantage of the method is to avoid self-subtraction of signal from the source, a problem becoming more severe</footnote>

<footnote>at shorter separations: on the other hand, this is obtained at the price of much poorer detection limits.</footnote>
<table><caption>Table D.1: Stellar companions identified in the sample with their astrometric and photometric properties. In addition to GPI measurements, literature astrometry is reported too.</caption><thead><tr><th>STAR</th><th>DATE</th><th>SEP<br>mas</th><th>PA<br>deg</th><th>CONTRAST<br>mag</th><th>ruwe</th><th>PMa<br>mas yr<sup>-1</sup></th><th>z<sub>ruwe</sub></th><th>z<sub>PMa</sub></th><th>MASS<br>M<sub>Jup</sub></th><th>source</th></tr></thead><tbody><tr><td rowspan="2">HIP 67199</td><td>2015-04-04</td><td>114 ± 6</td><td>354.3 ± 2.9</td><td>8.3<sup>a,b</sup></td><td rowspan="2">1.02</td><td rowspan="2">6.10 ± 0.03</td><td rowspan="2">0.2</td><td rowspan="2">37.4</td><td rowspan="2">&gt;20<sup>c</sup></td><td rowspan="2">TW<br>W23b</td></tr><tr><td>2019-03-07</td><td>147.28 ± 0.17</td><td>51.41 ± 0.06</td><td>—</td></tr><tr><td>HD 74341B</td><td>2015-12-20</td><td>744 ± 2</td><td>75.9 ± 0.3</td><td>4.3</td><td>0.74</td><td>—</td><td>1.5</td><td>—</td><td>530 ± 45</td><td>TW</td></tr><tr><td rowspan="2">HIP 26369</td><td>2018-01-06</td><td>155 ± 1</td><td>222.5 ± 0.6</td><td>3.9</td><td rowspan="2">3.80</td><td rowspan="2">10.40 ± 0.22</td><td rowspan="2">28.1</td><td rowspan="2">21.8</td><td rowspan="2">125 ± 35</td><td rowspan="2">TW<br>B22</td></tr><tr><td>2017-01-16</td><td>284.71 ± 0.3</td><td>201.78 ± 0.06</td><td>—</td></tr><tr><td rowspan="2">HD 24072</td><td>2018-11-22</td><td>193 ± 2</td><td>16.5 ± 0.4</td><td>5.3</td><td rowspan="2">1.97</td><td rowspan="2">—</td><td rowspan="2">12.2</td><td rowspan="2">—</td><td rowspan="2">355<sup>+95</sup><sub>-100</sub></td><td rowspan="2">TW<br>B22</td></tr><tr><td>2017-12-02</td><td>124.97 ± 0.77</td><td>7.81 ± 0.35</td><td>—</td></tr><tr><td rowspan="2">HIP 38160</td><td>2015-04-08</td><td>128 ± 6</td><td>283.4 ± 2.5</td><td>5.6<sup>b</sup></td><td rowspan="2">1.34</td><td rowspan="2">30.65 ± 0.09</td><td rowspan="2">2.7</td><td rowspan="2">92.3</td><td rowspan="2">&gt;240<sup>c</sup></td><td rowspan="2">TW<br>R13</td></tr><tr><td>2009-11-25</td><td>141 ± 13</td><td>117.08 ± 2.28</td><td>—</td></tr><tr><td rowspan="3">HIP 74696</td><td>2015-07-29</td><td>156 ± 2</td><td>357.8 ± 0.6</td><td>5.0</td><td rowspan="3">0.87</td><td rowspan="3">5.10 ± 0.03</td><td rowspan="3">1.3</td><td rowspan="3">28.9</td><td rowspan="3">380<sup>+52</sup><sub>-50</sub></td><td rowspan="3">TW<br>W23a</td></tr><tr><td>2019-08-11</td><td>139 ± 6</td><td>25.3 ± 2.3</td><td>5.3</td></tr><tr><td>2023-04-19</td><td>118.94 ± 0.18</td><td>56.94 ± 0.07</td><td>—</td></tr></tbody></table>

Notes. $^a$: no-ADI contrast; $^b$: upper limit; $^c$: lower limit. W23a: Waisberg et al. (2023a); W23b: Waisberg et al. (2023b); B22: Bonavita et al. (2022); R13: Rameau et al. (2013); TW: this work.

<figure><img src="image_12.png" /></figure>

Fig. D.1: Proper motion test for the five stellar companion candidates with multiple epochs. As in Figure 4, a filled star indicates the displacement expected for a bound source with no relative motion to the star, whereas an empty star marks the location of a static background source. Second epochs are labeled by a '2', third epochs by a '3'. The large deviations from the filled star are likely caused by orbital motion.
<figure><img src="image_13.png" /></figure>

Fig. D.2: Upper panel: flux maps (in contrast units) showing the stellar companion candidates detected with cADI. Lower panel: GaiaPMEX (sma, mass) maps, with contours outlining the area corresponding to the 68% and 95% confidence level. Photometric masses (dots) or lower limits (arrows) are overplotted for comparison. The HD 74341B map should be interpreted as a nondetection, the white area being incompatible with the absence of a signal.
<table><caption>Table D.2: Physical and orbital parameters used in GaiaPMEX.</caption><thead><tr><th>Parameter</th><th>type</th><th>bounds or law</th></tr></thead><tbody><tr><td>log M<sub>c</sub></td><td>uniform</td><td>log M<sub>c</sub> ± Δ log M<sub>c</sub></td></tr><tr><td>log a</td><td>uniform</td><td>log a ± Δ log a</td></tr><tr><td>e</td><td>uniform</td><td>0–0.9</td></tr><tr><td>ω</td><td>uniform</td><td>0–π</td></tr><tr><td>Ω</td><td>uniform</td><td>0–2π</td></tr><tr><td>φ</td><td>uniform</td><td>0–1</td></tr><tr><td>I<sub>c</sub></td><td>sin I<sub>c</sub> uniform</td><td>0–π/2</td></tr><tr><td>ω̄</td><td>normal</td><td>N(ω̄, σ<sup>2</sup><sub>ω̄</sub>)</td></tr><tr><td>M<sub>⋆</sub></td><td>normal</td><td>N(M<sub>⋆</sub>, σ<sup>2</sup><sub>M<sub>⋆</sub></sub>)</td></tr></tbody></table>

Notes. *e*: eccentricity; *ω*: periastron longitude; *Ω*: longitude of ascending node; *φ*: phase; *I<sub>c</sub>*: inclination; *ω̄*: parallax; *M<sub>⋆</sub>*: stellar mass.
# Appendix E: Effect of input assumptions on occurrence rates

We explored the dependence of the results derived in Section 3 on several input assumptions: the uncertainty on stellar age, the choice of the substellar evolutionary model, the degree of dis-equilibrium chemistry of planet atmospheres. In principle, all of them are expected to induce systematic deviations in the luminosity-mass relation, possibly impacting the reliability of the derived occurrence rates.

As a first step, we evaluated the impact of model selection by repeating the computations from Table 3 using the AMES-Cond models (Baraffe et al. 2003) and the Sonora Bobcat models (Marley et al. 2021). AMES-Cond models ignore the effect of dust opacity and are therefore more appropriate for objects with $T_{\text{eff}} \lesssim 1300$ K compared to fully dusty models such as the AMES-Dusty models (Baraffe et al. 2003). The derived completeness maps are shown on the left side of Figure E.1; the differences between completeness values are plotted on the right side. In this regard, we stress that, given that it is the mean detection probability across the (mass, sma) area $A$ that enters into Eq. 2, absolute differences are a more accurate proxy than relative differences when evaluating the impact of completeness maps variations on the derived frequency posteriors. Inspection of Figure E.1 clearly indicates that the discrepancies are the widest in the mass range [1, 5] $M_{\text{Jup}}$, and rapidly decrease at larger masses: this can be seen as a consequence of the stronger cooling rate of less massive objects, combined with the larger theoretical uncertainties at lower masses.

As a consequence of this observation, we expect the lower mass value selected to define $A$, $M_{\text{low}}$ to have a large impact on the accuracy of the results. We quantified this effect by computing occurrence rates under the three models for $M_{\text{low}} = [1, 2, 3, 4, 5] M_{\text{Jup}}$. As expected, the problem exacerbates for lower values of $M_{\text{low}}$ (Figure E.2). This test justifies our conservative choice for $A$ (Sec. 3.4).

Afterwards, we investigated the dependence of the results on the assumption of chemical equilibrium: in particular, we used two suites of ATMO models that assume 1) chemical equilibrium (ATMO-CEQ) or 2) strong chemical disequilibrium (ATMO-NEQ-S), that is, a different relation for the vertical mixing coefficient (Phillips et al. 2020). Given a certain H-band magnitude, the fractional mass difference, computed as a function of age and ATMO-NEQ-W mass ($m \in [1, 10] M_{\text{Jup}}$), can be as large as 30% compared to the chemical equilibrium case. The variation is larger at lower masses and larger ages, that is, at lower effective temperatures. The derived completeness maps, analogous to Figure E.1, are shown in Figure E.3. Given our careful choice of mass boundaries, it is possible to argue that completeness values within $A$ are dominated by projection effects rather than detection limits. Hence, we expect the occurrence rates to be fully consistent with those of the original analysis.

Finally, we provide similar completeness maps to quantify the dependence on age uncertainty: Figure E.4 shows the variation of the maps when assuming lower and upper values for stellar ages.

All the occurrence rates derived in this Section are visually compared in Figure E.5. It is evident that any doublet of estimates is compatible within the errors, making the estimates presented in this work robust against systematic effects.
<figure><img src="image_14.png" /><figcaption>Fig. E.1: Effect of model selection on survey completeness: maps assuming the Ames-COND model (top row) and the Sonora model (bottom row). Left panels show completeness maps, while right panels indicate the difference relative to the map used for the analysis. The green dashed box indicates our nominal choice of *A*.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Fig. E.2: Trend of the uncertainties in the derived occurrence rates with $M_{\text{low}}$. Each model – shown as a triplet (nominal ages, lower ages, upper ages) – is plotted in a different color. Horizontal offsets have been applied to each line for the sake of visualization.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>Fig. E.3: Effect of non-equilibrium chemistry on survey completeness: maps using the ATMO models assuming equilibrium chemistry (top row) and strong disequilibrium chemistry (bottom row). Left panels show completeness maps, while right panels indicate the difference relative to the map used for the analysis. The green dashed box indicates our nominal choice of *A*.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>Fig. E.4: Effect of age uncertainty on survey completeness: maps assuming lower (top row) and upper (bottom row) values for stellar ages. Left panels show completeness maps, while right panels indicate the difference relative to the map used for the analysis. The green dashed box indicates our nominal choice of $\mathcal{A}$.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Fig. E.5: Effect of model selection and age uncertainty on planet occurrence ($\mathcal{A}$ = [5, 13] $M_{Jup}$ × [10, 100] au): results for the entire sample (red squares), the BA subsample (blue diamonds), the FGK subsample (green circles) using the standard ATMO-NEQ-weak model (ATMO-NW); the same model with lower (ATMO-NW, L) and upper (ATMO-NW, U) ages; the ATMO model with no (ATMO-C) and strong (ATMO-NS) disequilibrium chemistry; the AMES-Cond model (AMES-COND) and the Sonora Bobcat model (SONORA). A Jeffreys prior is assumed (see Section 3.4).</figcaption></figure>