© A. Vigan et al. 2021

1 Introduction

In the past 20 yr, large-scale direct-imaging surveys for exoplanets have discovered approximately 60 substellar and planetary-mass companions around young nearby stars (see, e.g., Wagner et al. 2019). Early surveys were relatively modest in size, with samples of 50–100 stars, while ongoing surveys target 500–600 stars. The largest projects to date are the SpHere INfrared Exoplanets (SHINE) project conducted with SPHERE (Chauvin et al. 2017) and the Gemini Planet Imager (GPI) Exoplanet Survey (GPIES; Macintosh et al. 2015). SHINE and GPIES have added three new exoplanet detections to the growing cohort of directly imaged objects (Macintosh et al. 2015; Chauvin et al. 2017; Keppler et al. 2018) and several additional higher mass brown dwarfs (Konopacky et al. 2016; Cheetham et al. 2018). However, new exoplanet detections are just one goal of large-scale direct-imaging surveys. These surveys also provide key spectral and orbital characterization data for known exoplanets (e.g., De Rosa et al. 2016; Samland et al. 2017; Chauvin et al. 2018; Wang et al. 2018; Müller et al. 2018; Cheetham et al. 2019; Lagrange et al. 2019; Maire et al. 2019), and statistical constraints on the distribution of such objects at star–planet separations > 20 au (e.g., Kasper et al. 2007; Nielsen & Close 2010; Heinze et al. 2010; Janson et al. 2011; Vigan et al. 2012, 2017; Biller et al. 2013; Rameau et al. 2013; Brandt et al. 2014; Galicher et al. 2016; Lannier et al. 2016; Stone et al. 2018; Baron et al. 2019; Nielsen et al. 2019).

In particular, the current generation of surveys strongly constrains the distribution of wide (> 10 au) giant exoplanets and substellar companions to young stars (Reggiani et al. 2016; Nielsen et al. 2019). The distribution of gas giant planets and brown dwarf companions as a function of mass and orbital separation can provide insight into formation mechanisms because different formation channels (e.g., planet formation in a disk versus brown dwarf binary formation in a protostellar core) may dominate in different circumstances (mass ratio of companion to host, orbital separation, and total system mass). Contrast limits from these surveys, and the cohort of detected objects, can be used with a Bayesian approach to constrain the fraction of systems hosting planetary or substellar companions and to explore distribution functions of their architectures (semimajor axis, mass, or eccentricity distributions). Diverse demographic models can be tested: (a) parametric models based on a wide range of point estimates of frequency over fixed ranges of mass and orbital separation (e.g., Reggiani et al. 2016; Meyer et al. 2018), extrapolated to the mass and separation ranges probed by direct-imaging surveys; and (b) population synthesis models (e.g., Mordasini et al. 2009; Forgan & Rice 2013), which use numerical estimates based on physical theories of various formation mechanisms to predict a population of exoplanets, which can be compared to our observations.

Imaging surveys have yielded significantly fewer exoplanet detections than expected using extrapolations of radial velocity (RV) planet populations to larger semimajor axes (e.g., Cumming et al. 2008). These extrapolations predicted dozens of detections with optimistic assumptions. While this is disappointing from the perspective of detection, these results constrain the distribution of giant exoplanets and brown dwarf companions at separations > 10 au from host stars. SHINE achieves typical contrasts of 10⁻⁵ –10⁻⁶ at separations of 0.4–0.5′′ (Langlois et al. 2021, hereafter Paper II). Based on evolutionary models used for mass–luminosity conversion and on the ages and distances of targets in our sample (Desidera et al. 2021, hereafter Paper I), we expect that the SHINE survey will be sensitive to 1–75 Jupiter mass (M_Jup) companions at separations 5–300 au. We note, however, that as for all direct-imaging surveys, the sensitivity to the lowest masses improves for larger semimajor axes and is expected to reach a minimum only at a few dozen astronomical units. Predictions of what SHINE will see depend on the planet mass function, the orbital distribution, any correlations between the two, and perhaps on host star properties. Because only a small number of companions are detected (typically a few in a given large-scale survey), we must simplify the models to a few free parameters, preferably based on measured populations of substellar companions and extrasolar planets obtained by other methods (e.g., Meyer et al. 2018).

Several formation mechanisms can lead to the formation of 1–75 M_Jup companions that are detected in these surveys. In addition to formation channels for very low-mass binary companions (e.g., Kratter et al. 2010; Offner et al. 2010), companions can be formed in multiple modes as planets in circumstellar disks as well. The core accretion (CA) scenario is a bottom-up framework where a solid core of a few Earth masses forms first (Mizuno 1980; Pollack et al. 1996; Alibert et al. 2004), and then the rapid accretion of gas builds up gas giant planets (Piso et al. 2015b; Venturini et al. 2015). In contrast, the disk instability, or gravitational instability (GI), is a top-down binary star-like framework where planets form very quickly in the outer parts of disks from clumps that detach from the rest of the disk, become gravitationally bound, and contract into a giant planet (Boss 1998; Vorobyov 2013). Multiple theoretical approaches provide simulated populations of planets that formed through these mechanisms, which can then be compared to planets thatare detected through direct imaging. The comparison of direct-imaging observations with theoretical predictions was pioneered by Janson et al. (2011) and Rameau et al. (2013). Vigan et al. (2017) were then the first to compare observations to the outputs of population synthesis models based on the GI scenario. The authors found that these models can describe a fraction of the population, wide-orbit giant companions. With the improved sensitivity in mass and semimajor axis of new surveys, it becomes realistic to compare observations to predictions of both CA models (e.g., Mordasini et al. 2017) and GI models (e.g., Forgan & Rice 2013; Forgan et al. 2015).

In this paper we present a first statistical analysis of the properties of the population of 1–75 M_Jup companions at orbital separations from 5–300 au based on the first 150 stars observed in the SHINE survey. In Sect. 2 we present the target sample considered in our analysis, the detections that are taken into account, how the detection limits were derived and converted into mass, and finally, the survey sensitivity derived from the observations. In Sect. 3 we introduce the exoplanet population models to which we compare our observations, and in Sect. 4 we present the simulation tools we used for the comparison. In Sect. 5 we present all of our results, and finally, in Sect. 6 we discuss them in a broader context and compare the SHINE results to previously published surveys.

2 Statistical sample and detection limits

This section provides information regarding the sample of targets we considered (Sect. 2.1), the observations and data analysis (Sect. 2.2), how the planetary candidates were treated (Sect. 2.3), the statistical weight attributed to the detections (Sect. 2.4), and finally, the mass conversion of the detection limits (Sect. 2.5). The detailed properties of the statistical sample are separately treated in a companion paper (Paper I), and the observations and data analysis are discussed in a second companion paper (Paper II).

2.1 Statistical sample

The SHINE survey is being performed by the SPHERE consortium and exploits 200 nights of guaranteed time of observation (GTO). The main goal of SHINE is to observe a sample of 500 stars out of a larger sample of 800 nearby young stars to search for new substellar companions. The sample is oversized with respect to the available telescope time by a factor of approximately two on the basis of the adopted observing strategy, which assumes observations across meridian passage in order to achieve the maximum field-of-view (FoV) rotation for optimal angular differential imaging. This requires some flexibility in the target list in order to optimize the scheduling. Moreover, a sample of at least a few hundred objects is mandatory to achieve robust inference on the frequency of planets because the expected frequency of substellar companions is likely low (e.g., Vigan et al. 2017).

The sample includes a broad range of stellar masses to explore the effect of this parameter on planet frequency. The stellar masses in the sample range from ~3.0 M_⊙ to 0.3–0.5 M_⊙, and the faint end is determined by the working limit of the adaptive optics system of SPHERE (Sauvage et al. 2016; Beuzit et al. 2019). The 800 targets were divided into four priority bins, called P1, P2, P3, and P4 in decreasing order of priority, which were used to schedule the observations. Roughly a dozen targets of special interest were added to the sample for scientific reasons (presence of known substellar objects, disks, RV planets, etc.) and classified as P0 (highest priority). Some of these special objects were originally drawn from the 800-star sample built for the statistical analysis, but some were also added at a later stage. Both the selection of the targets from a wide database of nearby young stars and the priority ranking were based on simulations of planet detectability with SPHERE that was performed before the start of the survey (Paper I). These simulations adopted two different expected distributions and dependences on the stellar mass for the planet population in order to avoid biasing our results by relying on a single assumed distribution.

In addition to obvious magnitude and declination limits, the selection criteria exclude known spectroscopic and close visual binaries within the radial FoV of the SPHERE/IRDIS science camera (6′′). No selection was made in favor of stars with known disks or IR excess.

A total of 150 targets with first-epoch observations until February 2017 were included in the present analysis, and second-epoch observations extended until 2019. At this stage, the sample was not complete in any aspect, considering that the scheduling was optimized for the whole survey, but is fully representative of the whole sample. The sample includes 53 BA stars, 77 FGK stars, and 20 M stars. The median stellar age of the sample used in this early statistical analysis is 45 Myr (90% between 11 and 450 Myr), the median stellar mass is 1.15 M_⊙ (90% between 0.57 and 2.37 M_⊙), and the median distance is 48 pc (90% limits of 11 and 137 pc).

Most of the stars in the present sample belong to young moving groups. The age range considered for group members corresponds to the typically accepted spread of the mean age of the group. Age spread within the groups is not included, although mild kinematic outliers are considered individually, and their age uncertainties are typically larger than those of bona fide members. Therefore we expect that the effect is negligible for most groups, which show no clear evidence of a large age spread, and that it may be present only for targets in the Scorpius-Centaurus OB2 association, which are known to have an age spread of up to 50% (Pecaut & Mamajek 2016). However, the targets in this association constitute only 20% of our present sample, therefore the overall effect should be small. The general design of the survey, the sample selection, and the simulations performed for building it, and the parameters of the individual targets in this series of papers are fully described in Paper I.

2.2 Observations and data analysis

The complete observing strategy, data analysis, and detection performance for the targets in the sample are described in Paper II. All observations were performed with the SPHERE instrument at the Very Large Telescope (VLT; Beuzit et al. 2019) in either IRDIFS or IRDIFS-EXT mode, that is, the two near-IR (NIR) subsystems, IFS and IRDIS, observed in parallel. The IFS covers a 1.7′′ ×1.7′′ FoV and IRDIS covers a circular unvignetted FoV of diameter ~9′′. Some targets were observed multiple times because of known companions and/or the detection of (new) candidate companions. This varied the observations for each target.

All the data were downloaded at the SPHERE data center (Delorme et al. 2017a) and processed with the SpeCal software (Galicher et al. 2018) for speckle suppression, derivation of detection limits, and astro-photometry of the detected candidates. The final data products were then transferred to a private part of the DIVA+ database¹ (Vigan et al. 2017) from where they were retrieved for our analysis.

More specifically, we used the 5σ IRDIS and IFS detection limits of each observation for all the targets in the sample. These detection limits are derived based on the noisein the speckle-subtracted image, compensated for the throughput of the algorithm (calibrated with simulated planet injections), the transmission of the coronagraph (calibrated from measurements in SPHERE), and the small sample statistics (Mawet et al. 2014). More details are provided in Galicher et al. (2018) and Paper II.

2.3 Planetary candidates

For 91 of the 150 targets in the sample, candidates were identified in the first-epoch observations. Valid follow-up observations were obtained for 45 targets, complemented by archival or published data for a total of 39 targets. The use of archival data enabled us to recover the position of a significant number of candidates on a very long temporal baseline and minimize the need for follow-up observations (see Paper II). Multi-epoch data were therefore obtained for 66 targets, which left only 25 targets without follow-up. In practice, follow-up observations were attempted for all of these 25 targets but could not be performed because of scheduling problems or poor weather during the observing runs.

The stellar proper motion of the targets with candidates in the sample is 81 ± 64 mas yr⁻¹, with a minimum and maximum of 13 and 454 mas yr⁻¹, respectively. Follow-up observations were only scheduled after a time span that would unambiguously enable us to distinguish between a bound companion and background source, resulting in temporal baselines of 1.94 ± 1.22 yr. The motion of the stars with candidates over their respective time baselines is 156 ± 145 mas for the SHINE data, with a minimum and maximum of 10 and 500 mas yr⁻¹, respectively (even extending beyond a few arcseconds when the archival data are considered). With a typical astrometric accuracy of a few mas, the follow-up observations were therefore distant enough in time to reliably assess the status of candidates.

We were not always able to obtain a clear confirmation of companion or background status for the 66 targets with follow-up observations, sometimes despite multi-epoch follow-up of some candidates. This is in most cases due to the non-redetection of some candidates because the observing conditions between epochs varied. Follow-up observations of all remaining candidates is foreseen in a future dedicated programme.

Of the 1454 individual candidates, 16 were confirmed as companions or were already known to be companions (Table 1), 1134 were confirmed as background using either relative astrometry or classification based on color-magnitude diagrams (see Paper II), but 304 remain unconfirmed, sometimes despite multiple observations. This count is largely dominated by one target close to the galactic plane (TYC 7879-0980-1) that has only one epoch and more than 100 candidates in the IRDIS FoV, and a handful of other targets with a few dozen candidates. Based on the distance of the targets in the sample, the projected separation of unconfirmed candidates ranges from 3 to 1300 au. In Appendix A we provide a cumulative histogram of the number of candidates as a function of projected semimajor axis to illustrate that a small number of targets largely contributes to the total number of undefined candidates.

The statistics of young substellar companions beyond 300 au has been well established in the past decade by numerous direct-imaging surveys that used various instruments and targeted stars with a wide range of properties (e.g., Nielsen & Close 2010; Heinze et al. 2010; Janson et al. 2011; Vigan et al. 2012, 2017; Rameau et al. 2013; Galicher et al. 2016; Lannier et al. 2016; Stone et al. 2018; Baron et al. 2019). One of the main goals of the new generation of exoplanet imagers such as SPHERE is to set the first meaningful constraints below 100 au. To further this goal, we restricted our analysis to projected semimajor axes ≤ 300 au. With this upper limit in terms of physical separation, the number of unconfirmed candidates decreases to only 187, again mostly clustered around a handful of targets.

For our targets with incomplete follow-up and unconfirmed candidates, two approaches can be followed in the statistical analysis. Either the detection limits for individual targets can be raised above the level of the brightest unconfirmed candidate, regardless of its position in the FoV, or the limit can be cut at the separation of the closest unconfirmed candidate. For this analysis we chose the latter approach in order to retain the best possible sensitivity at the closest separations.

2.4 Statistical weight of detections

In the construction of the complete SHINE sample, each individual target was attributed a scientific priority from P1 to P4 based on planet detectability simulations (Paper I). The additional very high priority bin (P0) was also created to enforce the observation of specific targets, such as those with already known companions or with companions detected in parallel with SHINE, but by other teams (e.g., GPIES or open-time programs with SPHERE). The original scientific priority of stars with detections are listed in Table 1, along with the reassignment to the P0 priority when applicable.

The reassignment of P0 priority to some targets based on a priori knowledge of the presence of companions necessarily introduces astatistical bias. Without a priori knowledge, these targets may not have been immediately observed in the course of the survey, or may even have had a very low probability of being observed (e.g., HIP 107412, Milli et al. 2017). To properly take into accountthe previously known detections in the sample, we introduced a statistical weight related to the probability that the target would have been observed if the companion had not been known before. The value of this weight is between zero and one.

The most straightforward case is for completely new detections around HIP 65426 (Chauvin et al. 2017) and HIP 64892 (Cheetham et al. 2018). These two targets were in the P1 category and were observed as part of the normal course of the survey. Each of these detections are therefore counted as full detections (statistical weight of 1.0).

Then, there are cases of targets that were known to have a companion, but were not given a higher priority based on this knowledge. Only two such targets are included in the current sample: η Tel and CD -35 2722. These objects were observed independently of the fact that there was knowledge about a substellar companions, and we can safely assume that if their companions had not been known, they would certainly have been detected. The latter statement is not a strong assumption because of the relatively low contrast and large angular separation of η Tel B (Lowrance et al. 2000) and CD -35 2722 B (Wahhaj et al. 2011). Each of these detections are therefore also counted as full detections (statistical weight of 1.0).

Finally, there are cases of targets for which the priority was boosted to P0 because of previously known companions (HIP 78530, β Pic, HR 8799, HD 95086, PZ Tel, AB Pic, and GSC 8047-0232) or because of the discovery of a companion by another team after the start of the SHINE survey (51 Eri and HIP 107412). For these stars, the assigned statistical weight is equal to the probability that a star from the same priority bin (P1 to P4) would have been observed by a fixed date, specifically, the date where the early SHINE statistical sample was frozen. Because the current sample was frozen in the course of the survey, this date also corresponds to the time where 100% of first-epoch observations were obtained for the stars in the sample. Following this analysis, detections around targets that originally were in the P1, P2, and P4 priority bins² were attributed a statistical weight of 0.60, 0.35, and 0.01 respectively. The weight values were computed numerically a posteriori, based on the definition of the sample and on the dates of all the SHINE observations. For example, a weight of 0.6 implies that 60% of the stars within the original priority class of that particular star were observed at the point at which the survey was frozen for the analysis, independently of the stellar types.

The statistical weight of each detection considered in the analysis is taken into account in the Markov chain Monte Carlo (MCMC) simulations. These are described in Sect. 4.

2.5 Mass conversion of the detection limits

To convert the detection limits obtained in luminosity space into mass detection limits, it is necessary to use a mass-luminosity relationship, L(M). Whereas for old (≳1 Gyr) systems this relationship is essentially unique for gas giants, at young ages, the value of the post-formation luminosity still remains uncertain (Marley et al. 2007; Spiegel & Burrows 2012; Marleau & Cumming 2014; Bonnefoy et al. 2014a,b). In recent years, first steps toward predicting the post-formation luminosity of planets have been taken (Berardo et al. 2017; Berardo & Cumming 2017; Cumming et al. 2018; Marleau et al. 2017, 2019b). While detailed predictions are not quite available yet, these theoretical studies suggest that warm or hot starts are more likely (see also the discussion in Sect. 5.2). This agrees with observational results that cold starts are disfavored for massive companions (see review in e.g., Nielsen et al. 2019). For lower-mass companions, the question remains open from the observational side. For example, a cold start is allowed by the data for 51 Eri b, but they do not exclude a hot start either (e.g., Rajan et al. 2017; Samland et al. 2017).

When luminosity was converted into mass, we used hot starts as the fiducial model, but also consider warm starts in Sect. 5.2 and Appendix B. Specifically, we took the Bern EXoplanet cooling tracks (BEX) coupled with the COND atmospheric models (Allard et al. 2001) and assumed hot-start or warm-start initial conditions (Marleau et al. 2019a). Extending the fits of Marleau et al. (2019a, Eqs. (1b) and (1c), respectively), we took the post-formation (i.e., initial) luminosity L_pf   of the BEX-hot and BEX-warm tracks as a function of planet mass M_p as $$\begin{array}{l}{\tilde{M}}_n\equiv \frac{M_p}{n M_{\text{Jup}}},\\ L_{\text{pf}}{}^{\text{BEX-hot}}=2.62\times {10}^{-5}{\tilde{M}}_1^{1.4} L_{\odot },\\ L_{\text{pf}}{}^{\text{BEX-warm}}=\{\begin{array}{ll}4.35\times {10}^{-6}{\tilde{M}}_1^{0.5} L_{\odot },\hfill & M_p⩽10 M_{\text{Jup}}\hfill \\ 1.39\times {10}^{-5}{\tilde{M}}_{10}^7 L_{\odot },\hfill & 10<M_p/M_{\text{Jup}}⩽20\hfill \\ 1.74\times {10}^{-3}{\tilde{M}}_{20}^{1.4} L_{\odot },\hfill & 20 M_{\text{Jup}}⩽M_p.\hfill \end{array}\end{array}$$

Finally, we also considered the COND-2003 cooling tracks (Baraffe et al. 2003), which have even higher L_pf   (M_p) that is not based on any formation model.

2.6 Survey sensitivity

In order to constrain the statistical properties of our observed sample, we first converted the observed detection limits into the same parameter space as the models, that is, from projected separation to semimajor axis and from detection contrast to companion mass, so as to determine the completeness of the survey in terms of semimajor axis a and companion mass M_p. For each star, we defined a grid of semimajor axis and mass values uniformly distributed in log space, with 500 values ranging from 0.1 to 10 000 au in a, and 200 values between 0.1 to 100 M_Jup in M_p. For each cell in the grid, we generated 10⁴ companions with arguments of periastron and orbital phases drawn from uniform distributions, taking into account the orbital velocities along the orbit (i.e., considering the fact that an eccentric companion spends more time near apastron). We used a uniform distribution in inclination in order to simulate random orientations of orbits in space. For the eccentricity distribution, we considered the recent results derived by Bowler et al. (2020) for directly imaged exoplanets and brown dwarf companions. For this parameter we adopted a Beta distribution with parameters [α = 0.95, β = 1.30], which corresponds to the best fit to the full sample of wide substellar companions studied in Bowler et al. (2020).

For each simulated companion, we then computed the corresponding projected separation from the drawn orbital elements and the semimajor axis a of that grid point. We finally determined whether the companion is detectable in our observations by verifying that the mass value M_p of that cell lies above the contrast curve converted into mass of the considered star at the obtained projected separation (see Sect. 2.5), and that this projected separation value lies within the FoV for that star. The fraction of detectable companions in each grid cell provides the fractional completeness as a function of mass and semimajor axis for each star in our sample. Summing all derived completeness maps and dividing by the number of targets, we obtained the average 2D completeness of the survey. This task was repeated using the mass limits obtained with the various evolutionary models described in Sect. 2.5, and considering the nominal, minimum, and maximum ages of the stellar primaries. This provided a separate completeness map for each specific analysis to be performed.

Using the completeness maps for each of the targets in the sample, we computed the depth of search of the complete survey, which provides the number of stars around which the survey is sensitive for a given substellar companion mass and semimajor axis. This metric is useful for estimating the statistical strength of the results presented later. The depth of search for the 150 stars of our sample is presented in Fig. 1, based on the nominal stellar ages and the BEX-COND-hot models (see Sect. 2.5). The core of the sensitivity (>100 stars) reaches 7–9 au for objects >10 M_Jup. At lower masses, the sensitivity to the lowest masses around at least 100 targets is reached at ~100 au with a mass of ~3 M_Jup. Sensitivity to 1 M_Jup planets is only reached around ~30 stars at separations of 100–200 au.

The mean completeness map for the whole sample provides the average sensitivity of the survey, that is, the average probability of detecting an object of given mass and semimajor axis. This metric enables a direct comparison of the SHINE survey to surveys performed using the previous generations of instruments. In Fig. 2 we compare the sensitivity of SHINE with that of the NaCo-LP survey (Chauvin et al. 2015; Vigan et al. 2017), which were both computed using detection limits converted into mass using the COND-2003 evolutionary tracks (Baraffe et al. 2003). While the two surveys do not share strictly identical samples, they both target a large pool of relatively young nearby stars, so that the probability of detection in the mass versus semimajor axis space averaged over all targets is a good metric for comparison. Clearly, the new generation of instruments such as SPHERE provides a significant boost in sensitivity for 1–10 M_Jup planets in the 5–50 au range. However, the core of the sensitivity (probability >50%) still remains beyond 10 au, even for the most massive substellar companions (10–300 au for companions >10 M_Jup).

We also plot in Fig. 2 an estimate of the range of H₂ O and CO snow lines for the stars in the SHINE sample. The snow lines are estimated based on a parametric disk temperatureprofile as derived from the composition of Solar System bodies (Lewis 1974) and on observations of a large sample of protoplanetary disks (Andrews & Williams 2005, 2007a,b). The average evaporation temperatures for H₂ O and CO have been reported in Öberg et al. (2011), specifically, they are 135 and 20 K, respectively. Because protoplanetary disk physics and chemistry are complex, these estimates of the snow lines locations are approximate, but they enable a first-order comparison of the sensitivity of SHINE in locations that are important for giant planet formation. It is interesting to note that SHINE has some sensitivity to massive objects at the level of the water snow-line, which might constitute a turnover point in the giant planet occurrence rate (Fernandes et al. 2019), although the core of the sensitivity is shifted toward larger orbital separations. If the water snow-line is indeed a turnover point, the low detection rate of new planetary companions in the SHINE and GPIES surveys might qualitatively indicate that this turnover might apply to low masses where SHINE (this work) and GPIES (Nielsen et al. 2019) have little sensitivity.

Fig. 1 Depth of search of the SHINE survey for the 150 stars in the sample. The black and white contour lines give the numbers of stars around which the survey is sensitive to substellar companions as a function of mass and semimajor axis. The mass conversion of the detection limits is based on the nominal stellar ages and on the BEX-COND-hot evolutionary models (Marleau et al. 2019a). The colored circles represent the detected substellar companions in the sample. The color indicates the spectral type of the primary star (BA, FGK, or M). The size of the symbol is proportional to the weight of the detection in the statistical analysis (see Sect. 2.4 and Table 1 for details).

Fig. 2 Comparison of the sensitivities of the NaCo-LP (Vigan et al. 2017; dashed red lines) and SHINE (this work; solid black lines) surveys (with the current sample), based on the average probability of detecting a companion as a function of its mass and semimajor axis. The analysis is based on detection limits that were converted using the COND-2003 evolutionary tracks for both surveys. The contours for the NaCo-LP are not labeled but are the same as for SHINE, and correspond to equal levels of detection probability. The range of semimajor axes spanning the H2 O and CO snow lines for the stars in the sample are overplotted (see Sect. 2.6 for details).

3 Exoplanet population modeling

We here compare our observations to two different types of exoplanet population models. The first type is a parametric model based on inputs from both theoretical and observational work (Sect. 3.1), which aims at being a better representation than the simple power-law distributions in mass and semimajor axis used previously (e.g., Lafrenière et al. 2007; Kasper et al. 2007; Nielsen & Close 2010; Vigan et al. 2012). Although relatively straightforward, this remains a simplified parametric approach to describing the giant exoplanet population. The second type of model is based on exoplanet population synthesis models, which by definition rely on very detailed (although often simplified) physical modeling of the planet formation, interactions, and evolution (Sect. 3.2). The parametric and population model types both include a top-down binary star-like formation component and a bottom-up planet-like formation component in an attempt to capture different formation pathways for the observed detections in the SHINE sample.

3.1 Parametric model

We compared our results to a parametric model that was developed to explain a wide range of observations. Details of this model are presented in Meyer et al. (in prep.). We provide here an overview of the key features of the model for our needs in the context of the SHINE survey.

For each bin of stellar spectral type (BA, FGK, and M), the model comprises two parts that represent two different populations of substellar companions: one is a planet-like population, and the other is a binary star-like population. For each of these two parts, we considered different distributions of objects as a function of mass and semimajor axis: a distribution of planets as a function of orbital separation (a) and a planet mass function (b) for the planet-like population, and an orbital distribution of low-mass binary companions (c) and a companion mass ratio distribution (d) for the binary star-like population. The planet part of the model (a and b) and the binary star part of the model (c and d) require a different normalization. In principle, all parameters of the model can be fit to the data. However, because our survey includes a limited number of observations, we only fit the normalization of the planet part and binary part separately (two free parameters). These normalization factors represent the amplitudes, or the relative frequencies, of having a very low-mass binary-like companion or a planet-like companion. Combined, the resulting fit represents the total probability for a star to have one or more substellar companions.

For part(a), the orbital distribution of gas giant planets, we assumed a Gaussian distribution in log a, a being the semimajor axis, with fixed mean and sigma. These properties likely depend on host star mass, and based on results to date, we adopted a log-normal distribution with mean loga = 0.45 and σ = 0.52 for M stars (Meyer et al. 2018; Fernandes et al. 2019), loga = 0.58 and σ = 0.69 for FGK stars, and loga = 0.79 and σ = 0.77 for BA stars (Meyer et al., in prep.). For part (b), the planet mass function, we assumed a power law where the frequency f depends on the ratio of the companion to the host star mass, q = M_p∕M_⋆, that is, f ∝ q^β, with β = −1.31 for all stellar types (Cumming et al. 2008; Wagner et al. 2019). We furthermore assumed that the planet mass function does not depend on orbital separation. The amplitude factor associated with the product of these two functions, f_PPL∕LN, is the first of our two fit variables. The planet part of our phenomenological model, which combines (a) and (b), is abbreviated PPL/LN (planet power-law, log-normal) from here on.

For part(c) we assumeda log-normal surface density of binary companions, as measured for stellar masses (e.g., Raghavan et al. 2010 for FGK stars and Winters et al. 2019 for M dwarfs) with mean log a = 1.30 and σ = 1.16 for M dwarfs, loga = 1.70 and σ = 1.68 for FGK stars, and loga = 2.59 and σ = 0.79 for BA stars (De Rosa et al. 2014). For part (d) we assumed a universal companion mass ratio distribution, which is roughly flat with the mass ratio (power-law slope of 0.25; Reggiani & Meyer 2013). We assumed that the companion mass ratio distribution extends to the minimum mass for fragmentation (cf. Reggiani et al. 2016) and that the companion mass ratio distribution does not depend on orbital separation. The other amplitude factor associated with the product of these two functions, f_BDB, is our second fit variable. The binary part of our phenomenological model, which combines (c) and (d), is abbreviated BDB (brown dwarf binary) from here on.

An illustrative comparison of the output populations with the survey sensitivity around FGK stars is provided in Fig. 3. The BDB and PPL/LN parts of the model are clearly visible: the density of planetary companions (PPL/LN) is highest at low masses, with a peak at a few astronomical units and a density decreasing toward higher masses and larger orbital separations, while the density of binary companions (BDB) is highest for higher masses and then slowly decreases toward planetary masses.

In our analysis, we fit only the relative frequencies f_BDB and f_PPL∕LN for the parametric model and for each bin of stellar spectral type (BA, FGK, and M). We also computed the total frequency for the sum of the planetary and binary parts of the model, f_{BDB +PPL∕LN}.

3.2 Population model

The population model consists of two different population synthesis models based on the GI scenario and the CA scenario, which are described in Sects. 3.2.1 and 3.2.2, respectively. Combined, they comprise the full population model. These models are currently computed only for solar-mass stars, therefore we compare them only to the observations of FGK stars in the sample (see Sects. 3.2.3 and 5). Comparison with higher and lower mass stars will be the subject of future work.

3.2.1 GI population

The synthetic GI populations are based on those first presented by Forgan & Rice (2013) and then updated by Forgan et al. (2018). These models involved running, in advance, a suite of 1D disk models that smoothly proceed from an epoch in which the GI dominates their evolution (Rice & Armitage 2009) to an epoch in which it is dominated by an alternative angular momentum transfer mechanism, such as the magnetorotational instability (Balbus & Hawley 1991). These models also include photoevaporation, which plays an important role in disk dispersal (Owen et al. 2011). The outer radius of each disk was taken to be 100 au, which optimizes the likelihood of the disk to undergo fragmentation, after which dynamical interactions can then sculpt the semimajor axis distribution.The disk-to-star mass ratios varied from 0.125 to 0.375, and the host star masses varied from 0.8 to 1.2 M_⊙.

To generate the synthetic populations, a disk model was selected and fragments were then placed in this disk. The innermost fragment was placed at the smallest radius where fragmentation is possible, typically beyond ~50 au (Rafikov 2005; Clarke 2009), and the subsequent fragments were then placed at separations that were initially a random number of Hill radii (uniform distribution between 1.5 and 3 Hill radii). The fragment masses were set by the local Jeans mass, their radii were set using the assumption that they are equivalent to the initial radii of star-forming cores, and their initial temperatures were set to be the virial temperature (Nayakshin 2010).

The fragments then followed a tidal downsizing process where they contracted and cooled, and evolved through disk migration and n-body interactions.Grains within the fragment can grow and sediment, potentially forming a solid core. When the radius of an embryo exceeds its Hill radius, it can be tidally disrupted, potentially allowing for the emergence of a terrestrial-mass protoplanetary core. Each system was evolved for a duration of 1 Myr to ascertain the effect of object–object scattering on the planetary orbital parameters (Forgan et al. 2015). Although each system was evolved for a time that is shorter than the observed ages of the objects to which we would like to compare to disk fragment models, this relatively short simulation time was used partly to reduce computational expense and partly because systems that produce scattering events express this instability within a few ten thousand years (Chambers et al. 1996; Chatterjee et al. 2008).

This process was repeated many times to produce a large population of planetary systems that have formed via GI. These systems were then used as input for the SHINE simulations for comparison with our observational results. The relative frequency of systems with at least one companion associated with the GI model of formation is noted f_GI from here on.

Fig. 3 Comparison of the depth of search of the SHINE survey for the 77 FGK stars in the sample with a population of 20 000 draws from our parametric model presented in Sect. 3.1. The contour lines give the numbers of stars around which the survey is sensitive to substellar companions as a function of mass and semimajor axis. The PPL/LN part of the model is represented with shades of red (low density of companions) to yellow (high density of companions), and the BDB part of the model is represented with shades of white (low density of companions) to blue (high density of companions). Only the detections around FGK stars are plotted.

3.2.2 CA population

The synthetic CA populations were obtained using the new Bern generation 3 model of planetary formation and evolution described in Emsenhuber et al. (2020a), which corresponds to an update of the model presented in Mordasini (2018). This model in turn has evolved out of earlier versions of the Bern model described in Alibert et al. (2004), Mordasini et al. (2012), and Benz et al. (2014). The model self-consistently evolves a 1D gas disk, the dynamical state of the solids, the accretion by the protoplanets, gas-driven migration of the protoplanets, the interiors of the planets, and their dynamical interactions. The specific population we used is population NG76 from the new-generation planetary population synthesis (NGPPS) series.

For the gas disk, the model assumes that it is viscously evolving (Lynden-Bell & Pringle 1974) and the macroscopic viscosity is given by the standard α parameterization (Shakura & Sunyaev 1973). The vertical structure was computed using a vertically integrated approach (Nakamoto & Nakagawa 1994) that includes the effect of stellar irradiation. We included additional sink terms for the accretion by the planets, and both internal and external photoevaporation, following Clarke et al. (2001) and Matsuyama et al. (2003), respectively.

The model assumes planetesimal accretion in the oligarchic regime (Ida & Makino 1993; Ohtsuki et al. 2002; Thommes et al. 2003). The model solves the internal structure equations (Bodenheimer & Pollack 1986) for the gas envelope. In the initial (or attached) phase, the envelope is in equilibrium with the surrounding disk gas, and accretion is governed by the ability of the planet to radiate the gravitational energy released from the accretion of both solids and gas. When the accretion rate exceeds the supply from the disk, the envelope is no longer in equilibrium with the disk and contracts (Bodenheimer et al. 2000). Planets undergo gas-driven migration, and the dynamical interactions are followed by means of an n-body simulation.

After 20 Myr, the model transitions into the evolution stage, where the planets are followed individually up to 10 Gyr. In this stage, the model computes the thermodynamical evolution of the envelope, atmospheric escape, and tidal migration, but the gravitational interactions with other planets in the system are not considered.

To obtain a synthetic population, we followed the procedure outlined in Mordasini et al. (2009) and Emsenhuber et al. (2020b). The distributions for the disk masses follow Tychoniec et al. (2018), and we used the relationship described by Andrews et al. (2010) to determine the characteristic radius that defines the radial distribution of the gas. The inner edge of the disk is based on the work of Venuti et al. (2017), with a log-normal distribution in period with a mean of 4.7 d. The dust-to-gas ratio was obtained as described in Mordasini et al. (2009) from the observed stellar [Fe/H], but we used the primordial solar metallicity as a reference (Lodders 2003) without an enhancement factor. The initial slope of the surface density of solids is steeper than the slope of the gas disk, following Ansdell et al. (2018).

The population used here consists of 1000 systems with 1 M_⊙ stars. Each disk started with 100 planetary embryos of lunar mass (10⁻² M_⊕), whose initial positions were randomly selected between the inner edge of the disk up to 40 au, with a uniform probability in the logarithm of the semimajor axis.

The generated systems were then used as input for the SHINE simulations. The relative frequency of systems in which at least one companion is associated with the CA mode of formation is noted f_CA from here on.

Fig. 4 Comparison of the depth of search of the SHINE survey for the 77 FGK stars in the sample with the population synthesis models based on the CA and GI formation scenarios presented in Sects. 3.2.1 and 3.2.2, respectively. The contour lines give the numbers of stars around which the survey is sensitive to substellar companions as a function of mass and semimajor axis. The CA companions are represented with shades of red (low density of companions) to yellow (high density of companions), and the GI companions are represented with shades of white (low density of companions) to blue (high density of companions). The apparent lower density of CA objects arises because the vast majority of the CA population is located outside the range of mass and semimajor axis considered in this plot. Only the detections around FGK stars are plotted.

3.2.3 Full population model

The two population synthesis models described above were combined to form the full population model. Because the population synthesis models were computed only for solar-mass stars, we restricted our analysis with this model to the 77 FGK stars that are part of the present SHINE sample. In our analysis we fit the relative frequencies f_CA and f_GI that are associated with the CA and GI parts of the model, respectively, and the total frequency for the sum of the two parts, f_GI+CA.

An illustrative comparison of the output populations with the survey sensitivity is provided in Fig. 4. Similarly to what has been described in Vigan et al. (2017) for the GI population, a large cluster of massive objects (>10 M_Jup) is located at separations of 50–100 au where the SHINE survey is the most sensitive. In contrast, the CA population only shows a rather small population of 1–30 M_Jup objects that are scattered at separations ranging from a few up to a few dozen astronomical units.

4 Statistical tools

We used a statistical tool based on the MCMC sampling method described in Fontanive et al. (2018, 2019) to constrain the companion fractions of our observed sample. The tool was built using the emcee (Foreman-Mackey et al. 2013) python implementation of the affine-invariant ensemble sampler for MCMC (Goodman & Weare 2010). The code was adapted to use two separate exoplanet population models, each made of two parts: the parametric model presented in Sect. 3.1, and the population model described in Sect. 3.2. In all simulations, the shapes of the underlying companion distributions in mass and semimajor axis were fixed to those of the models, leaving as only MCMC parameters the relative companion frequencies of the two model populations considered. The companion fractions f₁ and f₂ of populations 1 and 2, respectively, are defined over fixed semimajor axis and companion mass ranges, [a_min, a_max] and [M_p,min, M_p,max]. We sought the posterior distributions of f₁ and f₂ given our observed data, where 1 and 2 designate the two parts of our models, either BDB and PPL/LN for the parametric model, or GI and CA for the population model.

In order to take the uncertainties on the measured masses and semimajor axes of the detected planets and brown dwarfs around the observed targets into account, we followed the method of sampling approximation to the marginalized likelihood from Hogg et al. (2010). This offers a powerful approach in the framework of Bayesian statistics to inform a population-level likelihood using the posterior distributions of individual systems. At each step in the MCMC, K = 10³ sets of semimajor axes and masses are generated for each of the N_comp detected companions. Values are drawn from Gaussian distributions centered on the measured masses and semimajor axes, with Gaussian widths set to the uncertainties of the measurements (Table 1). When no uncertainties are available, the measured value is always chosen. When the 1σ interval is asymmetric around the most likely value, we defined a two-piece Gaussian (see, e.g., Wallis 2014) from which values were randomly selected. When the drawn values are between [a_min, a_max] and [M_p,min, M_p,max], a companion was counted towards the detections in that region of the parameter space. For each iteration k, we thus obtained a number N_sys,k of systems with at least one companion in the probed range, which might vary when the drawn parameters occasionally fell outside the ranges of interest. We note that N_sys,k may be smaller than the total number of detected companions when multiple planets or brown dwarfs are found around the same star.

For all iterations, we started by estimating the total number of companions expected to be detected in our observations. This was done by drawing simulated companions between [a_min, a_max] and [M_p,min, M_p,max] from the two model distributions, and injecting them into the combined completeness map defined in Sect. 2.6, using only the targets considered in a specific analysis (e.g., the BA, FGK, or M stars) with the selected evolutionary models and stellar ages. For the parametric model, N₁ = N₂ = 10⁴ companions were drawn from the continuous separation and mass ratio distributions describing the BDB and PPL/LN populations. We used the mean stellar mass of the studied subset to convert the mass ratios of the model companions into corresponding companion masses. When we worked with the synthetic population models (FGK stars only), we injected all companions found in each model within the considered semimajor axis and mass limits, adding up to totals of N₁ and N₂ companions, respectively. The expected total number of detections λ around the observed targets is then given by $$\lambda =\left[\frac{f_1}{N_1}{\sum }_{i=1}^{N_1}{p}_i+\frac{f_2}{N_2}{\sum }_{j=1}^{N_2}{p}_j\right]\times N_{\star },$$(2)

where N_⋆ is the number of stars in the studied subsample, and p_i and p_j are the probabilities of detecting simulated companions i and j from model populations 1 and 2, given the survey sensitivity. The first term in the square brackets thus provides the fraction of detectable companions from population 1 with companion frequency f₁, and the second term the fraction from population 2 with companion frequency f₂. The sum of these two terms gives the total fraction of companions that can be detected in the survey from the injected populations. This value was then multiplied by N_⋆ to obtain the total number of companions expected to be detected for respective companion frequencies f₁ and f₂ for the two parts of the model population.

The number λ of expected substellar detections may then be compared to the observed number of systems N_sys,k using Poisson statistics, as was done in Fontanive et al. (2018), providing a value ${\mathcal{L}}_{\text{P},k}$ at each step k. Averaged over the K iterations, this provides the first part of the likelihood function, which allows us to constrain the overall companion fraction. As detailed in Sect. 2.4, some of the detections are weighted to correct for observational biases due to the presence of previously known companions. The total number N_sys,k of detected systems is thus given by the sum of the effective detection rates for the companions to retain, listed in Table 1 (counting the HR 8799 system only once).

The second part of the likelihood compares the position of the companions in the mass–semimajor axis space to the model distributions in order to scale the relative companion frequencies of the two populations. This was done by defining a joint 2D probability density describing the semimajor axis–mass distributions of the combined model populations, weighting each part of the model by taking into account the relative companion fractions of each population, f₁ and f₂. Following the approach from Fontanive et al. (2018), we were then able to compute the probabilities of the detected companions being drawn from this overall model distribution. The full model probability density was convolved with the completeness map for the targets we investigated, so as to represent the distribution of companions that could be observed given the survey sensitivity. When the semimajor axis and mass log spaces are divided into bins of width 0.2 dex, the probability of observing a companion in a specific mass–semimajor axis bin is given by the volume below the probability density function delimited by the edges of that bin. For each of the N_comp detections in the considered subsample, we thus computed at each of the K steps the integral within the bin enclosing the drawn mass and semimajor axis of the companion. When the drawn values for a companion fall outside the considered parameter space, the integral was set to 0. For each detection n, the probability that this companion is drawn from the joint model distribution is given by averaging the integrals obtained for each iteration k (${\mathcal{L}}_{\text{nk}}$).

The final likelihood $\mathcal{L}$ is then computed as $$\mathcal{L}=\left({\displaystyle \prod _{n=1}^{N_{\text{comp}}}\frac{1}{K}}{\sum }_{k=1}^K{\mathcal{L}}_{nk}\right)\times \left(\frac{1}{K}{\sum }_{k=1}^K{\mathcal{L}}_{\text{P},k}\right),$$(3)

where ${\mathcal{L}}_{\text{nk}}$ is the integral computed above for the nth observed giant planet or brown dwarf companion at the kth iteration, and ${\mathcal{L}}_{\text{P},k}$ is the Poisson likelihood of the total number of detected systems at the kth step (N_sys,k) for an expected number of systems λ given by Eq. (2). The term in the left set of brackets hence gives the probability that the detected companions are drawn from the overall model distribution, and the probability for each companion is calculated as the average over the K iterations, and the final value for this term given by the product of the value for each detection. The term on the right provides the probability for detecting the number of companions that the survey yielded, taken again as the average of the Poisson likelihoods over the K iterations.

We adopted uniform priors between 0 and 1 for the two companion fractions, f₁ and f₂. The combination of the prior distributions and likelihood function according to Bayes’ theorem allows for the calculation of the posterior distribution for our two model parameters f₁ and f₂. In each simulation presented below, the code was run with 10³ walkers taking 10⁴ steps each. The initial 500 steps were discarded to remove the so-called burn-in phase, as a mean acceptance fraction was reached after some hundred steps.

5 Results

In this section we present the results from our comparison of the two models described above to our observations. We begin with results from our base parametric model (Sect. 5.1), then explore variations of the results as a function of input assumptions (Sect. 5.2), and conclude with tests of specific planet formation models using our population model (Sect. 5.3).

5.1 Frequency of substellar companions from parametric models versus host star mass

Our basic parametric model explores the companion frequency as a function of the companion mass ratio distribution and orbital distribution. For each host star, the mass range over which we are sensitive (1–75 M_Jup) results in a unique range of mass ratio q (approximately 0.0005–0.2). We explored the companion frequency from 5 to 300 au and used the age-dependent mass-luminosity conversions from the BEX-COND-hot model (see Sect. 2.5). We performed the mass-luminosity conversion for the nominal age, and the minimum and maximum ages, presented in Paper I. Our fitting explores the best-fit combinations of relative frequencies for the brown dwarf binary companion model (BDB) and the planet distributions (PPL/LN). All the results are presented in Fig. 5, and the corresponding point estimates of frequency (probability that one or more companions lie the quoted ranges) are reported in Table 2.

Figure 5 presents the probability density function (PDF) of the integrated frequency that one or more companions lies within 1–75 M_Jup and 5–300 au derived for the BDB and PPL/LN parts of our parametric model, and for the combined model (BDB + PPL/LN), as a function of spectral type. The frequency estimate for the planet contribution is significantly higher for the higher mass BA stars in our sample than for the lower mass M dwarfs. This is consistent with the idea that for 1–75 M_Jup the range of q probed for higher mass stars (0.0005–0.036 for 2 M_⊙) is lower for the planet mass function (dN_p ∝ q^−1.3) than in lower mass stars (q = 0.0015–0.228 for a 0.3 M_⊙ star). Similarly, the frequency of brown dwarf companions is much higher for low-mass M dwarfs than for higher mass BA stars. This reflects the fact that the binary brown dwarf companion mass range probed in our survey (1–75 M_Jup) is at higher q for lower mass stars than for higher mass stars (dN_BD ∝ q^0.25). This is qualitatively similar to results from the GPIES survey (Nielsen et al. 2019).

In this framework, where we combine two components in the model that represent planet-like and star-like formation pathways, it is interesting to study the degeneracy between the two components of the model. In Fig. 6 we show the degeneracy between the relative frequencies derived for the individual parts of the model for the BA, FGK, and M stars in the sample. For the BA and M stars, the observations are well explained by only a single part of the model, either PPL/LN or BDB. This is extremely clear for BA stars, where the likelihood is clustered close to zero for f_BDB and significantly elongated for f_PPL∕LN. While this is slightly less pronounced for M stars, the fact that only an upper limit can be derived for f_PPL∕LN is an indication that the contribution of the BDB part to the model is higher than that of the PPL/LN part. We note, however, that this result may be due to the small size of the M sample (20 stars) and the single detection we have in that subset (cf. Lannier et al. 2016). At higher confidence levels, similar probabilities are found in the correlation plot for roughly equal contributions from both parts of the model, and for either part being the predominant underlying population.

The result is more nuanced for the FGK stars, which appear as a transition between BA stars, dominated by planet-like formation over this range of q for companions, and M stars, dominated by binary star-like formation. FGK stars have a comparable contribution from both parts of the model in this range of q, but the total together is a lower frequency overall than either the BA or the M dwarf subsample. The contribution of the BDB and PPL/LN parts of the model is clearly inverted with respect to the BA stars: the BDB part dominates and the PPL/LN part has a small contribution, but the two parts are still required to fit the data best. While it is difficult to determine the formation scenario of individual objects, at the population level, the observed companions around FGK stars are therefore most likely explained by a combination of planet-like and star-like formation scenarios. To fully understand the transition from BA to M, it would be interesting to work in smaller bins of stellar spectral types, but the current data do not allow this because the overall number of detections is small.

Finally, our results appear to show a local minimum in the frequency of substellar companions around FGK stars. Because our sample contains only 20 M stars, we caution that this result should not be overinterpreted. The analysis of the full SHINE sample at the end of the survey will provide much stronger constraints based on a subsample of M stars that is two to three times larger than the current one.

Fig. 5 Probability density functions of the frequencies of substellar companions around BA (left), FGK (center), and M stars (right) based on the parametric model, computed for companions with masses in the range Mp = 1–75 MJup and semimajor axes in the range a = 5–300 au, and using the BEX-COND-hot evolutionary tracks for the mass conversion of the detection limits. Each plot shows the PDFs for the relative frequencies of the two components of the model (fBDB and fPPL∕LN), and for the total frequency for the full model (fBDB +PPL∕LN). The plain lines show the PDFs for the nominal stellar ages, while the shaded envelopes show the variation of these PDFs for the maximum and minimum stellar ages. The median values and 68% confidence intervals are provided in Table 2.

5.2 Effect of input assumptions

Our results are based on some important assumptions and parameters that need to be evaluated and tested with additional simulations. For these tests, we used as a reference the FGK sample and converted the detection limits into mass using the BEX-COND-hot models. All the results are summarized in Table 2.

Fig. 6 Correlation plots and marginalized probability density functions for fBDB and fPPL∕LN in the parametric model around BA (left), FGK (center), and M (right) stars computed for companions with masses in the range Mp = 1–75 MJup and semimajor axes in the range a = 5–300 au, and using the BEX-COND-hot evolutionary tracks at the optimal stellar ages. Contour lines in the correlation plots correspond to regions containing 68, 95, and 99% of the posterior, respectively. For the FGK subsample, the scale of the axes is different from that of the two other subsamples.

5.2.1 Stellar ages

Some important parameters are the stellar ages and time of planet formation because giant gaseous exoplanets are expected to slowly cool down and therefore eventually decrease in overall luminosity (Burrows et al. 1997; Baraffe et al. 2002; Fortney et al. 2011; Linder et al. 2019). We assumed that the planet age is equal to that of the star. The full age derivation for the sample is presented in Paper I. We study here the variation in PDFs of the frequency of substellar companions when the minimum and maximum ages for all the stars of the sample are compared to the nominal PDF. The corresponding depths of search plots at the differentages are provided in Appendix B. In Fig. 5 the plain lines show the PDFs for the nominal age of the stars, while the shaded regions around each curve show the PDFs computed assuming the minimum and maximum ages for the stars. Generally speaking, the effect of the ages on our results can be considered negligible, with changes of less than 2% in the peak frequency of companions (or upper limit) between thenominal ages and the minimum or maximum ages, mainly because most of the stars in the sample (~80%) are members of young nearby moving groups for which the allowed age range is narrow and well established. For targets that are not part of such moving groups, the age range is generally much larger, but because these targets form only a small fraction of the targets, the overall effect on the PDFs remains low.

5.2.2 Initial entropy

Another major astrophysical assumption is the initial entropy that is used as an input for the evolutionary models. Several studies in the past decade have demonstrated the significant effect of the assumed energy transfer method during the gas accretion phase onto the protoplanetary embryos (Marley et al. 2007; Fortney et al. 2008; Spiegel & Burrows 2012; Marleau & Cumming 2014). An extreme outcome is that the entire energy of the infalling gas is transformed into thermal energy by the accretion shock front without radiative losses, so that the entropy remains high; this leads to a bright planet post-formation. If conversely, the entire energy is radiated away at the shock, the entropy of the postshock gas is much lower, which leads to a faint planet at the end of its formation. These two extreme scenarios are generally known as “hot start” and “cold start”, respectively (Marley et al. 2007), and in reality, a whole range of intermediate initial entropy levels exists that are known as “warm starts” (Spiegel & Burrows 2012; Marleau & Cumming 2014). The most recent advanced models (Marleau et al. 2017, 2019b; Berardo et al. 2017; Berardo & Cumming 2017; Cumming et al. 2018) and global formation calculations (Mordasini 2013; Mordasini et al. 2017) clearly suggest that the classical (very) cold starts first proposed by Marley et al. (2007) are unlikely and would occur for a very small fraction of planets that are formed by CA. The initial luminosity of young Jupiters may strongly correlate with the size of their core (Mordasini 2013), with a realistic core mass associated with high entropy even within nominally cold gas accretion (Mordasini et al. 2017).

In order to test the effect of the post-formation entropy and luminosity, we converted our detection limits into mass using three sets of evolutionary tracks: BEX-COND-warm and BEX-COND-hot, which are described in Sect. 2.5, and COND-2003 (Baraffe et al. 2003). These are the standard tracks that have been used by most studies in the past. They correspond to an even hotter start than BEX-COND-hot. The corresponding depths of search plots for these various models are presented in Appendix B, and a comparison of the PDFs of the frequency of substellar companions is presented in Fig. 7. The effect of varying the model and/or initial entropy for the evolutionary tracks is negligible. This is perfectly in line with the results presented in Mordasini et al. (2017), who showed that the luminosity distributions of planets for the BEX-COND-hot and BEX-COND-warm models are extremely similar, and that the BEX-COND-hot tracks are equivalent to those of the COND-2003 tracks after a few million years (see also Fig. 1 of Marleau et al. 2019a). One should finally note that the COND-2003 models assume initial conditions that are arbitrary and not based on a formation model, in contrast to the BEX models.

The results are presented only for FGK stars, but the same conclusion applies for BA and M stars. We conclude that our results are independent of the choice of initial entropy in the evolutionary models.

5.2.3 Semimajor axis cutoff

Finally, another important parameter is the range of semimajor axes that we used to estimate the companion frequency. Our baseline uses a range extending from 5 to 300 au. As explained in Sect. 2.3, the outer limit of this range is primarily driven by the result of previous direct-imaging surveys, which have already constrained the frequency ofcompanions well down to 300 au. However, the inner limit of the range is more arbitrary and should be driven by the sensitivity of the survey. A lower limit that is too high (e.g., 50 au) will provide reliable constraints because in 50–300 au the sensitivity of SHINE is excellent: more than 100 out of 150 targets are sensitive down to 3 M_Jup (Fig. 1). It will include only about half of the detected companions, however, and will therefore have less statistical significance in a regime that is already dominated by small-number statistics. In contrast, a limit that is too low (e.g., 1 au) will provide looser constraints because the sensitivity at small separations is lower: only five targets provide sensitivity at ~1 au. It is only applicable for masses higher than ~40 M_Jup.

Our final choice of 5 au is set by the detections around β Pic (companion at a = 9 au) and HIP 107412 (companion at a = 6.7 au), but in termsof sensitivity to low masses, fewer than five observations are sensitive to the lowest estimated mass for HD 95086 b (2 M_Jup) at this separation. Starting at ~10 au, about a dozen of our observations have a sensitivity down to 2 M_Jup, and at ~30 au, about a dozen have a sensitivity down to 1 M_Jup. In Fig. 8 we show the variation in PDF of the frequency of substellar companions around FGK stars when the lower limit cutoff in semimajor axis is changed from 5 to 10 au. In the latter case, where one detection is removed from the analysis (HIP 107412 B, around an F5 star), the PDF is slightly modified; the peak frequency for the full parametric model shifts from 5.8 to 5.5%. The two peak frequencies remain fully compatible within their respective 68% confidence intervals: 3.0–10.5% for the 5–300 au analysis and 2.8–9.5% for the 10–300 au analysis. This result demonstrates that our conclusions are reliable in the selected range of semimajor axes. The plots for BA and M stars are also provided in Appendix C.

Although the effect of the semimajor axis cutoff appears to be stronger than that of the stellar ages and the initial entropy, with variations of the peak frequencies up to 5%, it does not change the conclusions we drew in Sect. 5.1. The observed trends remain the same even though detections are removed when a higher cutoff is chosen. This is a strong confirmation of our conclusions.

Fig. 7 Probability density functions of the frequencies of substellar companions around FGK stars based on the parametric model, computed for companions with masses in the range Mp = 1–75 MJup and semimajor axes in the range a = 5–300 au, and using the BEX-COND-hot (plain line), BEX-COND-warm (dashed line) or COND-2003 (dotted line) evolutionary tracks for the mass conversion of the detection limits. The median values and 68% confidence intervals are provided in Table 2.

Fig. 8 Probability density functions of the frequencies of substellar companions around FGK stars based on the parametric model, computed for companions with semimajor axes in the range a = 5–300 au (plain line) or a = 10–300 au (dashed line), and using the BEX-COND-hot evolutionary tracks for the mass conversion of the detection limits. The median values and 68% confidence intervals are provided in Table 2. The same plots for BA and M stars are shown in Fig. C.1.

5.3 Frequency of substellar companions from formation models

In addition to the detection and study of substellar companions, one of the main goals of the SHINE and GPIES campaigns has always been to provide meaningful constraints for planet formation models, or at least to distinguish between different formation scenarios for different categories of objects. In a previous work based on a sample of 200 FGK stars (the NaCo-LP; Vigan et al. 2017), we compared our direct-imaging observations with population synthesis models based on the GI formation scenario. The sensitivity of the NaCo-LP observations was, however, not sufficient to reach a regime of mass and semimajor axis where CA would have been a viable formation scenario (Fig. 2), hence the focus on GI at the time. With the improved SHINE sensitivity at small semimajor axes and low masses, combined with new-generation CA models, it is now possible to compare our observations with outcomes of both GI and CA population synthesis models, as is qualitatively illustrated in Fig. 4.

Here we compare our observations with a combination of CA and GI population synthesis models. From the theoretical point of view, using a superposition of these two formation scenarios is a reasonable assumption. The bottom-up CA formation pathway is very powerful in explaining properties of the exoplanet population within 5–10 au, but it faces great difficulties in explaining the formation of giant planets farther out than 10–20 au because the formation timescales that would be involved are prohibitively long (Alibert et al. 2005; Kennedy & Kenyon 2008). Although the pebble accretion process has been proposed as a way to solve the problem (Lambrechts & Johansen 2012; Levison et al. 2015), simulations show that this mechanism does not really form giant planets (several M_Jup) at large orbitaldistances (Bitsch et al. 2015, their Figs. 4 and 5). Perhaps the only viable scenario to place CA-formed planets on very wide orbits is to invoke scattering between multiple planets in the systems that originate from different initial embryos (Veras et al. 2009; Dawson & Murray-Clay 2013; Marleau et al. 2019a). The populations described in Sect. 3.2.2 include multiple embryos and their subsequent interactions during the early evolution of the protoplanetary disk (20 Myr). Despite the scattering, very few objects are scattered out to distances of several dozen or hundreds of au where some detections are observed (see Fig. 4). The top-down GI formation pathway more readily explains theexistence of gas giants on wide orbits, but the conditions that can lead to disk fragmentation are still not fully understood (Meru & Bate 2011; Paardekooper 2012; Rice et al. 2012, 2014; Young & Clarke 2016). Disk-planet interactions (Kley & Nelson 2012) or planet-planet scattering certainly affect the original semimajor axis distribution of protoplanets and result in exoplanets that cover a wide range of possible masses, sizes, locations, and compositions. These effects are also taken into account in the populations described in Sect. 3.2.1.

Because the population synthesis models described in Sect. 3.2 are currently computed only for solar-mass stars, we based our analysis on the 77 FGK stars from the sample, and the five detections of substellar companions around such stars. Figure 9 shows the PDF of the frequency of substellar companions based on the synthetic population models. The peak f_GI+CA model is located at 5.7%, with a 68% confidence interval of 2.9–9.5%. Interestingly, the corner plot showing the correlation between the two components of the model in Fig. 10 looks different from the plot for the parametric model in Fig. 6. While it is not possible to draw quantitative conclusions here, the CA contribution appears to be greater than the GI part. Both parts of the model are still required to explain the observations, as is visible from the roughly triangular shape of the 2D posterior, but the shape is narrower and more elongated in the direction of CA. The CA part of the model therefore contributes a slightly larger fraction of the full companion population that is required to explain the data. At this stage, it would certainly be necessary toextend our analysis to BA and M stars to confirm if the trends identified in Sect. 5.1 hold when based on physical population models rather than empirical parametric models. This will be explored in future work using the full SHINE sample, with population synthesis models computed for higher and lower mass stars.

As a cross-check with our parametric model we overplot in Fig. 11 the PDFs of the full models. With peak frequencies at 5.8 and 5.7%, and 68% confidence intervals of 3.0–10.5% and 2.9–9.5% for f_{BDB +PPL∕LN} and f_GI+CA, respectively,the results appear to be fully consistent between the two modeling approaches. This is expected because each model includes a combination of planet-like (i.e., bottom-up) and binary star-like (i.e., top-down) formation components.

Fig. 9 Probability density functions of the frequencies of substellar companions around FGK stars based on the population model, computed for companions with masses in the range Mp = 1–75 MJup and semimajor axes in the range a = 5–300 au, and using the BEX-COND-hot evolutionary tracks for the mass conversion of the detection limits. Each plot shows the PDFs for the relative frequencies of the two components of the model (fGI and fCA), and for the total frequency for the full model (fGI+CA). The plain lines show the PDFs for the nominal stellar ages, and the shaded envelopes show the variation of these PDFs for the maximum and minimum stellar ages. The median values and 68% confidence intervals are provided in Table 2.

Fig. 10 Correlation plots and marginalized PDFs for fGI and fCA in the population model around FGK stars, computed for companions with masses in the range Mp = 1–75 MJup and semimajor axes in the range a = 5–300 au, and using the BEX-COND-hot evolutionary tracks at the optimal stellar ages. Contour lines in the correlation plots correspond to regions containing 68, 95, and 99% of the posterior, respectively.

6 Discussion and conclusion

6.1 Comparison to previous works

The SHINE survey is certainly one of the deepest, and it is one of the first to open the low-mass regime at semimajor axes 5–50 au. This enables us to obtain quantitative statistical constraints in that range. It also offers some overlap with the parameter space that has been explored by previous works, which have placed strong statistical constraints on the population of young giant planets on wide orbits. In this section we assess the compatibility of some of these previous works with the SHINE results. The comparison is not completely straightforward because the considered ranges of mass, semimajor axes, and stellar spectral types vary from one study to the next. In addition, numerous studies have so far either made no physical assumptions on the underlying population of planets, often considering flat distributions in mass and semimajor axes, or used power-law parametric models that sometimes were simply extrapolated from RV surveys and truncated to avoid a continuous growth of the number of planets at wide orbits. The latter approach is contradicted by the latest observational results, which show indications for a turnover in the frequency of companions at the snow line (e.g., Fernandes et al. 2019) and a negative power-law distribution in semimajor axis at wide orbital separations (e.g., Nielsen et al. 2019). In most cases, our comparison is therefore more a consistency check than a quantitative comparison.

The comparison of our SHINE results with previous studies is presented in Table 3. For this comparison we have recomputed the frequency of systems based on our parametric model while trying to match the other parameters as closely as possible: mass range, semimajor axis range, or stellar spectral type bins. This is not always possible, and some caveats are inevitable. Nonetheless, the estimates from previous surveys are generally compatible with the new values derived for SHINE for the differentstellar spectral types. To estimate the compatibility, we tested the null hypothesis that the two measurements are equal with a5% risk, as described in Appendix D. In most cases we find that measurements are compatible with each other within this 5% risk. In only two cases are the values not compatible with SHINE for BA stars: the GPIES analysis from Nielsen et al. (2019), which we discuss in more detail in the next paragraph, and the meta-analysis from Bowler (2016). For the latter, it might be argued that their sensitivity in the quoted ranges of mass and semimajor axes around BA stars is marginal at best, which has a strong effect on the frequency that they derive. Because of the sensitivity of SHINE at small semimajor axes, our results can be considered far more robust.

In contrast to Bowler (2016), the sensitivity of the GPIES survey (Nielsen et al. 2019) below 100 au is comparable to that of SHINE. Their frequency estimation for BA stars using a uniform distribution as a prior clearly contradicts our own estimation, but uniform distributions as a prior are not physically realistic. It is more reasonable to compare our derived value with the value they derivedusing their power-law parametric model. Similarly to our parametric model, their model aims at realistically modeling the underlying population of substellar companions with just a few parameters. They modeled the planet population with a distribution of the form $$\frac{\text{d}{N}^2}{\text{d}m\text{d}a}=f{C}_1{m}^{\alpha }{a}^{\beta },$$

where f is the frequency of planetary systems, and m and a the mass and semimajor axis of planets, respectively. Based on the GPIES data acquired for stars in the range 1.5–5 M_⊙, which corresponds to our BA sample, their best fit is obtained for $f={8.9}_{-3.6}^{+5.0}\%$, α = −2.37 and β = −1.99 for planets in the range 3–100 au and 2–13 M_Jup. Within the error bars, their value of f is indeed almost exactly equal to the value of f_{BDB +PPL∕LN} $={8.6}_{-4.5}^{+7.3}\%$ that we derive based on the SHINE data using our own parametric model. We therefore consider that the two surveys agree excellently for early-type stars within the assumptions of our respective parametric models. The agreement can partly be explained by the overlap in terms of targets (67 targets) and of detections in the two surveys (Paper I), but this partial overlap cannot by itself fully explain the agreement. We instead consider the agreement as confirmation that the frequency of substellar systems around young BA stars is indeed around 8%, regardless of the sample that is considered.

Fig. 11 Comparison of the PDF of the frequency of systems with at least one companion for the full parametric and population models, fBDB +PPL∕LN and fGI+CA, respectively.

6.2 Implications for formation theory

Based on our observation of 150 stars that are part of the full SHINE sample, we can already conclude that gas giant companions are more commonly found around higher mass stars, and brown dwarf binary companions are observed more frequently around lower mass stars. The same conclusion is reached by Nielsen et al. (2019) for the GPIES survey. This can be interpreted in the context of our parametric model, which treats the observed companion frequency in terms of mass ratio distributions that are the same for all stellar masses. Within the range of companion masses to which we are sensitive (1–75 M_Jup), the planet part of the companion mass function for higher mass stars samples planets at smaller q, predicting higher frequencies. Conversely, at larger q for lower mass primaries (e.g., >0.1 for 30 M_Jup companions to 0.3 M_⊙ stars), the brown dwarf companion part of the parametric model is expected to dominate, as observed. This suggests a formation framework that is independent of stellar mass, but for which different pathways dominate as a function of q.

We recall that the deuterium burning limit at 13 M_Jup likely plays no specific role in the physical evolution of young very low-mass companions as both GI and CA can produce objects that are above and below that threshold (Mollière & Mordasini 2012; Chabrier et al. 2014). It is likely that binary star-like formation and planet-like formation pathways both contribute to the substellar mass distribution, without strict physical limits between thetwo. In this framework, brown dwarfs constitute the low-end tail of the stellar companion mass ratio distribution extended to the very low-mass regime, as was initially suggested by Metchev & Hillenbrand (2009). This universality of the companion mass distribution has since then received some strong observational support (Reggiani & Meyer 2011, 2013). On the other hand, more planet-like formation models such as GI or CA (with or without pebble accretion) clearly suggest that massive companions at wide orbits can be formed, constituting the high-end tail of their mass distribution (Forgan & Rice 2013; Forgan et al. 2015; Mordasini et al. 2012; Lambrechts & Johansen 2012).

Our parametric analysis predicts only the relative probability that a given mass ratio q may have formedfrom binary star-like or planet-like processes without discontinuities. Ultimately, the characterization of planetary-mass companions might reveal their formation pathway, for instance, if the comparison of volatile abundances shows differences relative to the star, perhaps indicating a disk formation process (Öberg et al. 2011; Piso et al. 2015a; Mordasini et al. 2016). The models also suggest that the overall efficiency of gas giant planet formation and low-q binary formation does not depend strongly on stellar mass, although the estimates of companion frequency strongly depend on the range of q that is considered in the analysis, which in turn depends on the host star mass: while the stellar mass scales the planet mass function in our model (self-similar in q), the normalization constants are roughly consistent throughout the range of stellar masses we studied. In addition, we note that our original sample selection for SHINE (Paper I) has purposely removed all known visual binaries, which means that our observations are far from complete in some ranges of q. This caveat might be circumvented in the future using results from the ESA/Gaia (Gaia Collaboration 2016, 2018) survey to correct the completeness in q.

A comparison of our results to specific realizations of GI population synthesis suggests that this is probably not the dominant channel of planet formation around FGK stars. This conclusion has been reached by Vigan et al. (2017), and our new analysis based on deeper SHINE data only strengthens this result. Perhaps more importantly, the SHINE results show that GI may not be the dominant formation scenario even for the most massive companions at large distances (5–10 M_Jup companions at ~50 au). However, additional input physics and more sophisticated simulations might alter the comparison in the future. Furthermore, some variants of GI such as early fragmentation within infalling disks that are only partially supported by rotation (Stamatellos & Whitworth 2008; Stamatellos et al. 2011; Forgan & Rice 2012) or rapid inward migration might contribute significantly to low-q binary populations. However, overall, the CA population synthesis models appear to be more promising. Alternative initial conditions such as disk lifetimes that depend on host star mass, and other relevant processes such as core formation based on pebble accretion (Alibert 2017; Ndugu et al. 2018) will illuminate the robustness of these predictions.

Like many studies in the past, our results are affected by the choice of evolutionary tracks for the conversion of the detection limits in the luminosity space into mass limits in a physical space. Significant progress is currently made in this field, which provides alternatives to the canonical evolutionary tracks (e.g., Burrows et al. 1997; Baraffe et al. 2003) that usually give good estimates at old ages and high masses, but require further validation at young ages or low masses (e.g., Konopacky et al. 2010; Dupuy & Liu 2017). The current and next-generation population synthesis models are gaining the ability to quantitatively predict the post-formation luminosity of young planets, which in turn drives the early evolution of the planet (Mordasini et al. 2012; Mordasini 2013; Emsenhuber et al. 2020a,b). The predictions of these models provide a robust view of the range of post-formation luminosities that planets can take (Mordasini et al. 2017). Much work still remains to be done to understand and accurately model the physics during the accretion phase, including the accretion geometry onto the planet (e.g., Gressel et al. 2013; Szulágyi 2017; Batygin 2018; Béthune 2019; Schulik et al. 2019) or the radiative properties of the accretion shock (Marleau et al. 2017, 2019b), but the discovery and study of accreting protoplanets such as PDS 70 b and c (Keppler et al. 2018; Haffert et al. 2019; Thanathibodee et al. 2019; Aoyama & Ikoma 2019; Christiaens et al. 2019; Hashimoto et al. 2020) will certainly help. More generally, the continued direct detection of spatially resolved substellar companions that already have dynamical mass estimates either through RV and/or astrometry (e.g., Crepp et al. 2012; Bowler et al. 2018; Peretti et al. 2019; Dupuy et al. 2019; Brandt et al. 2020; ESA/Gaia DR4), will help calibrate the mass-luminosity relationships even more precisely and will give greater confidence in our survey results.

Finally, the combination of results from direct imaging, RV, transit, microlensing, astrometry, and timing variations will provide a more complete picture of exoplanet demographics as a function of stellar mass. Trends in the companion mass ratio distribution as a function of orbital radius might reveal important discontinuities. Future work shouldalso consider more carefully how planet populations depend on stellar multiplicity, and how global planet architectures affect our statistical results, such as the ratio of planet masses and orbital radii in multiplanet systems.

6.3 Summary and perspectives

We have presented the first statistical analysis of the properties of the population of substellar companions at wide orbital separation based on a subset of 150 stars from the SHINE survey. The full details of the sample, the observations, and the data analysis are presented in two companion papers (Papers I; II). Although the size of the current sample is only a fraction of the full SHINE sample, we can already derive some important conclusions. Based on our parametric model presented in Sect. 3.1, we draw the conclusions listed below.

• 1.

  We determine the frequency of systems in which at least one companion has a mass in the range M_p = 1–75 M_Jup and a semimajor axis in the range a = 5–300 au, f_{BDB +PPL∕LN}, to be ${23.0}_{-9.7}^{+13.5}\%$, ${5.8}_{-2.8}^{+4.7}\%$, and ${12.6}_{-7.1}^{+12.9}\%$ for BA, FGK, and M stars, respectively. These values were derived using a conversion of the detection limits into mass using the BEX-COND-hot evolutionary tracks and the nominal age for all the stars in the sample. These values are average estimates over the stated ranges, but the sensitivity at the lowest masses and shorter separations is limited. This means that the uncertainties increase significantly when we focus on low masses and short separations.

• 2.

  The frequency of substellar companions is significantly higher around BA stars than around FGK and M stars, by factors of approximately 4 and 2, respectively. The apparent local minimum in frequency around FGK stars is suggestive and will be examined in detail with our full survey sample in the future.

• 3.

  Our two-component parametric model shows a clear inversion between BA and M stars. While in the case of BA stars the likelihood of the f_PPL∕LN part of the model dominates f_BDB, this former part only becomes an upper limit for M stars. This can be translated physically into a predominance of the planet-like formation pathway for companions detected around early-type stars over the binary star-like formation pathway for the mass ratio range that is sampled as a function of host star type.

• 4.

  The FGK stars are a transition range of spectral types, where observations are better explained by a combination of the two parts of the model, although the contribution of the PPL part is small. While it would be extremely interesting to perform an analysis in smaller bins of spectral type, the current data do not allow this because we have only a few detections.

• 5.

  The input assumptions such as the stellar ages, cutoff in semimajor axis, or evolutionary tracks, have little effect on the frequencies of planetary systems derived from the observations.

• 6.

  We find that the frequency of systems in which at least one companion has a mass in the range M_p = 2–13 M_Jup and a semimajor axis in the range a = 3–100 au around BA stars is ${8.6}_{-4.5}^{+7.3}\%$. This value is fully compatible with the value derived by GPIES (Nielsen et al. 2019) in a survey with similar sensitivity as SHINE, but with a slightly different sample (Paper I). This confirms the reliability of the estimation.

Based on the population model presented in Sect. 3.2, we can also draw the following conclusions for FGK stars:

• 7.

  We determine the frequency of systems in which at least one companion has a mass in the range M_p = 1–75 M_Jup and a semimajor axis in the range a = 5–300 au, f_GI+CA, to be ${5.7}_{-2.8}^{+3.8}\%$ for FGKstars. This value was derived using a conversion of the detection limits into mass using the BEX-COND-hot evolutionary tracks and the nominal age for all the stars in the sample, but again this result is not very sensitive to the input parameters. The same words of caution as in item 1 apply here.

• 8.

  Qualitatively, the contribution of the CA part of the model appears to be larger than the GI part, which means that CA contributes a higher fraction of the full companion population required to explain the data. Simulations extended to BA and M stars are required, however, to determine whether the general trend highlighted in item 2 above holds when weconsider our population model instead of the parametric model.

• 9.

  The values of f_GI+CA and f_{BDB +PPL∕LN} perfectly agree. Although the underlying model is different, the overall frequency values required to explain the observations are almost the same with the two approaches.

The SHINE survey is due to be completed in 2020, but will certainly extend over a few more years to become complete in terms of follow-up for all candidates within at least a 300 au and possibly even farther away. The final sample will include over 600 stars, which will make SHINE the largest high-contrast imaging survey to date, covering from B to M stars in the solar neighborhood. Beyond the reanalysis of the complete SHINE datawith advanced post-processing techniques (Cantalloube et al. 2015; Ruffio et al. 2017; Flasseur et al. 2018), which will hopefully provide improved detection limits at small separations, the full power of the survey will be in the statistical conclusions based on a sample that is almost four times larger than the sample we used here. Some of the prospects for future statistical inference work include an extension of our analysis based on population synthesis models to BA and M stars, the analysis of subsamples such as stars with disks or know infrared excess (e.g., Wahhaj et al. 2013) or stars that belong to nearby young moving groups (e.g., Biller et al. 2013), or an extension of the completeness of the sample in q space using the ESA/Gaia DR4 results.

Acknowledgements

SPHERE is an instrument designed and built by a consortium consisting of IPAG (Grenoble, France), MPIA (Heidelberg, Germany), LAM (Marseille, France), LESIA (Paris, France), Laboratoire Lagrange (Nice, France), INAF – Osservatorio di Padova (Italy), Observatoire de Genève (Switzerland), ETH Zürich (Switzerland), NOVA (Netherlands), ONERA (France) and ASTRON (Netherlands) in collaboration with ESO. SPHERE was funded by ESO, with additional contributions from CNRS (France), MPIA (Germany), INAF (Italy), FINES (Switzerland) and NOVA (Netherlands). SPHERE also received funding from the European Commission Sixth and Seventh Framework Programmes as part of the Optical Infrared Coordination Network for Astronomy (OPTICON) under grant number RII3-Ct-2004-001566 for FP6 (2004–2008), grant number 226604 for FP7 (2009–2012) and grant number 312430 for FP7 (2013–2016). This work has made use of the SPHERE Data Centre, jointly operated by OSUG/IPAG (Grenoble), PYTHEAS/LAM/CeSAM (Marseille), OCA/Lagrange (Nice), Observatoire de Paris/LESIA (Paris), and Observatoire de Lyon/CRAL, and supported by a grant from Labex OSUG@2020 (Investissements d’avenir – ANR10 LABX56). This research has made use of the Direct Imaging Virtual Archive (DIVA), operated at CeSAM/LAM, Marseille, France. A.V. acknowledges funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No. 757561). G.-D.M. acknowledges the support of the DFG priority program SPP 1992 “Exploring the Diversity of Extrasolar Planets” (KU 2849/7-1). C.M., A.E., and G.-D.M. acknowledge the support from the Swiss National Science Foundation under grant BSSGI0_155816 “PlanetsInTime”. Parts of this work have been carried out within the framework of the NCCR PlanetS supported by the Swiss National Science Foundation. A.-M.L. acknowledges funding from Agence Nationale de la Recherche (France) under contract number ANR-14-CE33-0018. T.H. and R.A.T. acknowledge support from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No. 832428). F.Me. acknowledges funding from Agence Nationale de la Recherche (France) under contract number ANR-16-CE31-0013. C.P. acknowledge financial support from Fondecyt (grant 3190691) and financial support from the ICM (Iniciativa Científica Milenio) via the Núcleo Milenio de Formación Planetaria grant, from the Universidad de Valparaíso.

Appendix A Undefined candidates

Figure A.1 provides a cumulative histogram of the number of undefined candidates as a function of projected semimajors axes. Globally, TYC 7879-0980-1 and HIP 63839 are responsible for more than 50% of the total number of undefined candidates: they are both close to the galactic plane and have only one single-epoch observation, which explains why they contribute such a large fraction of the total number of undefined candidates. In total, five targets are responsible for ~80% of undefinedcandidates. With a cutoff at 300 au, a total of 96 undefined candidates remain, which represents 32% of the total number. The way they are handled in the analysis is detailed in Sect. 2.3.

Fig. A.1 Cumulative histogram of the number of undefined candidates as a function of projected semimajor axis. Only the five targets that contribute the highest number of undefined candidates are labeled in the plot for clarity. With a cutoff of 300au for our analysis, only 96 undefined candidates remain (32% of the total number).

Appendix B SHINE depth of search

Figure B.1 shows the depth of search of the SHINE survey for the 150 stars in the sample based on different assumptionsfor the stellar ages and evolutionary models. The depth of search gives the number of stars in the sample around which the survey is sensitive for substellar companions as a function of mass and semimajor axis. We computed it using the nominal, minimum, and maximum ages and for the BEX-COND-warm and BEX-COND-hot (see Eq. (1)), and COND-2003evolutionary models.

Fig. B.1 Depth of search of the SHINE survey for the 150 stars in the sample computed using detection limits converted into mass using three different sets of evolutionary models (left: BEX-COND-warm; center: BEX-COND-hot; right: COND-2003) and for different stellar ages (top: minimum; middle: nominal; bottom: maximum). Each plot gives the numbers of stars around which the survey is sensitive for substellar companions as a function of mass and semimajor axis.

Appendix C Semimajor axis cutoff

Fig. C.1 Probability density functions of the frequencies of substellar companions around BA (left) and M (right) stars based on the parametric model, computed for companions with semi-major axes in the range a = 5–300 au (plain line) or a = 10–300 au (dashed line), and using the BEX-COND-hot evolutionary tracks for the mass conversion of the detection limits. The same plot for FGK stars is shown in Fig. 8.

The effect of the inner limit of the semimajor axis range is not identical for all spectral types due to the different number of detections. In Fig. 8 we show that for FGK stars the effect of changing the lower limit from 5 to 10 au is very weak, even though the detection around HIP 107412 is removed.

Figure C.1 shows the effect of changing the semimajor axis lower limit for BA and M stars. For M stars, the effect is negligible because the only detection around an M star (CD -35 2722) remains untouched since its semimajor axis is constrained in the range 76–216 au. However, the effect is more significant for BA stars, for which the detection around β Pic must be removed because its semimajor axis is tightly constrained within 8.5–9.2 au. The PPL/LN part of the model is the most affected, with a shift of the peak of the PDF from ~15 to ~11% when the cutoff changes from 5 to 10 au. We highlight that (a) β Pic b is clearly in the mass range dominated by the PPL/LN part of the model (in contrast to HIP 107412 B), and (b) the 5–10au range is also where the peak of the PPL/LN distribution is expected. These two elements combined can easily explain that the effect on the PPL/LN part of the model is larger for the BA stars than for the FGK stars. However, despite the larger effect of the semimajor axis cutoff for the BA stars, our conclusions remain unchanged.

Appendix D Comparison of two asymmetric distributions

Let X represent a random variable that has an asymmetric normal distribution. We call X₀ the mode of the distribution (i.e., value for which the distribution is maximum), and σ₋ and σ₊ the standard deviations toward values belowand above X₀, respectively.The PDF ϕ_X of X reads $$\begin{array}{lll}{\varphi }_X(x)\hfill & =\hfill & \frac{\sqrt{2}}{\sqrt{\pi }\left({\sigma }_{-}+{\sigma }_{+}\right)}[\text{mathds}{1}_{]-\infty ,X_0]}\exp \left(-\frac{{(x-X_0)}^2}{2{\sigma }_{-}^2}\right)\hfill \\ \hfill & \hfill & +\text{mathds}{1}_{]X_0,+\infty [}\exp \left(-\frac{{(x-X_0)}^2}{2{\sigma }_{+}^2}\right)],\hfill \end{array}$$(D.1)

where 1_[a,b] is the identity function between a and b, and it is zero elsewhere.

We assume that there are two independent measurements $x=X_0{}_{-{\sigma }_{X,-}}^{+{\sigma }_{X,+}}$ and $y=Y_0{}_{-{\sigma }_{Y,-}}^{+{\sigma }_{Y,+}}$ of the same measure. We then consider that x and y are realizations of the random variables X and Y that follow the asymmetric normal distributions ϕ_X and ϕ_Y defined by Eq. (D.1). We introduce the random variable Z = X − Y for which the probability density function ϕ_Z is $${\varphi }_Z(z)={\displaystyle {\int }_{-\infty }^{+\infty }{\varphi }_X}(t){\varphi }_Y(t+z)\text{d}t.$$(D.2)

From ϕ_Z, we can calculate the probability P(z_threshold) such that |z| < z_threshold $$P(z_{\text{threshold}})={\displaystyle {\int }_{-z_{\text{threshold}}}^{z_{\text{threshold}}}{\varphi }_Z}(t)\text{d}t$$(D.3)

and find the particular value z₉₅ such that $$P(z_{95})=\mathrm{0.95.}$$(D.4)

This value z₉₅ depends on the modes and standard deviations of ϕ_X and ϕ_Y.

As x and y are two independent measurements of the same measurand, we test the null hypothesis x = y. If the mode z_mode of ϕ_Z is such that |z_mode| < z₉₅, then we accept the null hypothesis with a 5% risk.

For each line of Table 3, we compare one SHINE measurement to another. We calculate ϕ_Z from the two asymmetric normal distributions, and then z_mode and z₉₅. If |z_mode| < z₉₅, we conclude that the two measurements are compatible with a 5% risk. In the other case, we conclude that the two measurements are not compatible with a 5% risk.

References

All Tables

All Figures

Fig. 1 Depth of search of the SHINE survey for the 150 stars in the sample. The black and white contour lines give the numbers of stars around which the survey is sensitive to substellar companions as a function of mass and semimajor axis. The mass conversion of the detection limits is based on the nominal stellar ages and on the BEX-COND-hot evolutionary models (Marleau et al. 2019a). The colored circles represent the detected substellar companions in the sample. The color indicates the spectral type of the primary star (BA, FGK, or M). The size of the symbol is proportional to the weight of the detection in the statistical analysis (see Sect. 2.4 and Table 1 for details).

In the text

Fig. 2 Comparison of the sensitivities of the NaCo-LP (Vigan et al. 2017; dashed red lines) and SHINE (this work; solid black lines) surveys (with the current sample), based on the average probability of detecting a companion as a function of its mass and semimajor axis. The analysis is based on detection limits that were converted using the COND-2003 evolutionary tracks for both surveys. The contours for the NaCo-LP are not labeled but are the same as for SHINE, and correspond to equal levels of detection probability. The range of semimajor axes spanning the H2 O and CO snow lines for the stars in the sample are overplotted (see Sect. 2.6 for details).

In the text

Fig. 3 Comparison of the depth of search of the SHINE survey for the 77 FGK stars in the sample with a population of 20 000 draws from our parametric model presented in Sect. 3.1. The contour lines give the numbers of stars around which the survey is sensitive to substellar companions as a function of mass and semimajor axis. The PPL/LN part of the model is represented with shades of red (low density of companions) to yellow (high density of companions), and the BDB part of the model is represented with shades of white (low density of companions) to blue (high density of companions). Only the detections around FGK stars are plotted.

In the text

Fig. 4 Comparison of the depth of search of the SHINE survey for the 77 FGK stars in the sample with the population synthesis models based on the CA and GI formation scenarios presented in Sects. 3.2.1 and 3.2.2, respectively. The contour lines give the numbers of stars around which the survey is sensitive to substellar companions as a function of mass and semimajor axis. The CA companions are represented with shades of red (low density of companions) to yellow (high density of companions), and the GI companions are represented with shades of white (low density of companions) to blue (high density of companions). The apparent lower density of CA objects arises because the vast majority of the CA population is located outside the range of mass and semimajor axis considered in this plot. Only the detections around FGK stars are plotted.

In the text

Fig. 5 Probability density functions of the frequencies of substellar companions around BA (left), FGK (center), and M stars (right) based on the parametric model, computed for companions with masses in the range Mp = 1–75 MJup and semimajor axes in the range a = 5–300 au, and using the BEX-COND-hot evolutionary tracks for the mass conversion of the detection limits. Each plot shows the PDFs for the relative frequencies of the two components of the model (fBDB and fPPL∕LN), and for the total frequency for the full model (fBDB +PPL∕LN). The plain lines show the PDFs for the nominal stellar ages, while the shaded envelopes show the variation of these PDFs for the maximum and minimum stellar ages. The median values and 68% confidence intervals are provided in Table 2.

In the text

Fig. 6 Correlation plots and marginalized probability density functions for fBDB and fPPL∕LN in the parametric model around BA (left), FGK (center), and M (right) stars computed for companions with masses in the range Mp = 1–75 MJup and semimajor axes in the range a = 5–300 au, and using the BEX-COND-hot evolutionary tracks at the optimal stellar ages. Contour lines in the correlation plots correspond to regions containing 68, 95, and 99% of the posterior, respectively. For the FGK subsample, the scale of the axes is different from that of the two other subsamples.

In the text

Fig. 7 Probability density functions of the frequencies of substellar companions around FGK stars based on the parametric model, computed for companions with masses in the range Mp = 1–75 MJup and semimajor axes in the range a = 5–300 au, and using the BEX-COND-hot (plain line), BEX-COND-warm (dashed line) or COND-2003 (dotted line) evolutionary tracks for the mass conversion of the detection limits. The median values and 68% confidence intervals are provided in Table 2.

In the text

Fig. 8 Probability density functions of the frequencies of substellar companions around FGK stars based on the parametric model, computed for companions with semimajor axes in the range a = 5–300 au (plain line) or a = 10–300 au (dashed line), and using the BEX-COND-hot evolutionary tracks for the mass conversion of the detection limits. The median values and 68% confidence intervals are provided in Table 2. The same plots for BA and M stars are shown in Fig. C.1.

In the text

Fig. 9 Probability density functions of the frequencies of substellar companions around FGK stars based on the population model, computed for companions with masses in the range Mp = 1–75 MJup and semimajor axes in the range a = 5–300 au, and using the BEX-COND-hot evolutionary tracks for the mass conversion of the detection limits. Each plot shows the PDFs for the relative frequencies of the two components of the model (fGI and fCA), and for the total frequency for the full model (fGI+CA). The plain lines show the PDFs for the nominal stellar ages, and the shaded envelopes show the variation of these PDFs for the maximum and minimum stellar ages. The median values and 68% confidence intervals are provided in Table 2.

In the text

Fig. 10 Correlation plots and marginalized PDFs for fGI and fCA in the population model around FGK stars, computed for companions with masses in the range Mp = 1–75 MJup and semimajor axes in the range a = 5–300 au, and using the BEX-COND-hot evolutionary tracks at the optimal stellar ages. Contour lines in the correlation plots correspond to regions containing 68, 95, and 99% of the posterior, respectively.

In the text

	Fig. 11 Comparison of the PDF of the frequency of systems with at least one companion for the full parametric and population models, fBDB +PPL∕LN and fGI+CA, respectively.
In the text

	Fig. A.1 Cumulative histogram of the number of undefined candidates as a function of projected semimajor axis. Only the five targets that contribute the highest number of undefined candidates are labeled in the plot for clarity. With a cutoff of 300au for our analysis, only 96 undefined candidates remain (32% of the total number).
In the text

Fig. B.1 Depth of search of the SHINE survey for the 150 stars in the sample computed using detection limits converted into mass using three different sets of evolutionary models (left: BEX-COND-warm; center: BEX-COND-hot; right: COND-2003) and for different stellar ages (top: minimum; middle: nominal; bottom: maximum). Each plot gives the numbers of stars around which the survey is sensitive for substellar companions as a function of mass and semimajor axis.

In the text

Fig. C.1 Probability density functions of the frequencies of substellar companions around BA (left) and M (right) stars based on the parametric model, computed for companions with semi-major axes in the range a = 5–300 au (plain line) or a = 10–300 au (dashed line), and using the BEX-COND-hot evolutionary tracks for the mass conversion of the detection limits. The same plot for FGK stars is shown in Fig. 8.

In the text