The SPHERE infrared survey for exoplanets (SHINE). III. The demographics of young giant exoplanets below 300 au with SPHERE

The SHINE project is a 500-star survey performed with SPHERE on the VLT for the purpose of directly detecting new substellar companions and understanding their formation and early evolution. Here we present an initial statistical analysis for a subsample of 150 stars that are representative of the full SHINE sample. Our goal is to constrain the frequency of substellar companions with masses between 1 and 75 MJup and semimajor axes between 5 and 300 au. We adopt detection limits as a function of angular separation from the survey data for all stars converted into mass and projected orbital separation using the BEX-COND-hot evolutionary tracks and known distance to each system. Based on the results obtained for each star and on the 13 detections in the sample, we use a MCMC tool to compare our observations to two different types of models. The first is a parametric model based on observational constraints, and the second type are numerical models that combine advanced core accretion and gravitational instability planet population synthesis. Using the parametric model, we show that the frequencies of systems with at least one substellar companion are $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. We also demonstrate that a planet-like formation pathway probably dominates the mass range from 1-75 MJup for companions around BA stars, while for M dwarfs, brown dwarf binaries dominate detections. In contrast, a combination of binary star-like and planet-like formation is required to best fit the observations for FGK stars. Using our population model and restricting our sample to FGK stars, we derive a frequency of $5.7_{-2.8}^{+3.8}\%$, consistent with predictions from the parametric model. More generally, the frequency values that we derive are in excellent agreement with values obtained in previous studies.

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 −5 -10 −6 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 that are 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.

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

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 secondepoch 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.

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 1 (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 noise in 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.

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 −1 , with a minimum and maximum of 13 and 454 mas yr −1 , 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 −1 , 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 followup 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 A&A 651, A72 (2021) Notes. (a) With a semimajor axis of ∼620 au and no additional published constraints, HIP 78530 is not taken into account into most of our simulations, which use a cutoff at 300 au. 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 directimaging 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. 2012Vigan et al. , 2017Rameau 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.

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 a statistical 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 account the 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, 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 2 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.

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 Cumming et al. 2018;Marleau et al. 2017Marleau et al. , 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 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.

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 4 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  Table 1 for details). 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 2 O and CO snow lines for the stars in the SHINE sample. The snow lines are estimated based on a parametric disk temperature profile 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. The average evaporation temperatures for H 2 O and CO have been reported in  (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 H 2 O and CO snow lines for the stars in the sample are overplotted (see Sect. 2.6 for details). Ö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.

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.

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 log a = 0.45 and σ = 0.52 for M stars (Meyer et al. 2018;Fernandes et al. 2019), log a = 0.58 and σ = 0.69 for FGK stars, and log a = 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 assumed a 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, log a = 1.70 and σ = 1.68 for FGK stars, and log a = 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 .

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.

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 starforming 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 terrestrialmass 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.

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) andEmsenhuber 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 −2 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.

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.

Statistical tools
We used a statistical tool based on the MCMC sampling method described in Fontanive et al. (2018Fontanive et al. ( , 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 1 and f 2 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 1 and f 2 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 3 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 1 = N 2 = 10 4 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 1 and N 2 companions, respectively. The expected total number of detections λ around the observed targets is then given by 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 1 , and the second term the fraction from population 2 with companion frequency f 2 . 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 1 and f 2 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 L 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, A72, page 9 of 22 A&A 651, A72 (2021) weighting each part of the model by taking into account the relative companion fractions of each population, f 1 and f 2 . 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 (L nk ).
The final likelihood L is then computed as where L nk is the integral computed above for the nth observed giant planet or brown dwarf companion at the kth iteration, and L 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 1 and f 2 . 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 1 and f 2 . In each simulation presented below, the code was run with 10 3 walkers taking 10 4 steps each. The initial 500 steps were discarded to remove the so-called burnin phase, as a mean acceptance fraction was reached after some hundred steps.

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

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 massluminosity 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.  Table 2. 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.

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.

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 different ages 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 the nominal 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.

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. , 2019bCumming 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 . 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, andCOND-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.

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  Table 2. the result of previous direct-imaging surveys, which have already constrained the frequency of companions 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 terms of 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.

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 orbital distances (Bitsch et al. 2015, their Figs. 4 and 5). Perhaps the only viable scenario to place CA-formed planets on very wide orbits A72, page 13 of 22 Probability density FGK f GI f CA f GI + CA 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 M p = 1-75 M Jup and semimajor axes in the range a = 5-300 au, and using the BEX-CONDhot 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 ( f GI and f CA ), and for the total frequency for the full model ( f GI+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. 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 the existence 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. 2012Rice et al. , 2014Young & 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 to extend 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.

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 Notes. The "Mass" and "S.m.a." columns give the ranges of companion masses and semimajor axes, respectively. (a) Compatibility between the results from SHINE and from the previous work. We assumed one asymmetric normal distribution for each measurement, and we tested the null hypothesis that the two measurements are equal with a 5 % risk, as described in Appendix D. A check mark indicates that the null hypothesis is accepted, and a cross mark that it is not. (b) The SHINE analysis is always truncated at 300 au. (c) In contrast to confidence intervals that are expressed at 68% confidence level, all upper limits are expressed at 95% confidence level. 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 different stellar spectral types. To estimate the compatibility, we tested the null hypothesis that the two measurements are equal with a 5% 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 derived using 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 A72, page 15 of 22 A&A 651, A72 (2021) a distribution of the form 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 +5.0 −3.6 %, α = −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 +7.3 −4.5 % 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.

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 lowmass 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 the two. 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. 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 formed from binary starlike 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 lowq 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 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; 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 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 . 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. , 2019b, but the discovery and study of accreting protoplanets such as PDS 70 b and c 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 should also 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.

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 +13.5 −9.7 %, 5.8 +4.7 −2.8 %, and 12.6 +12.9 −7.1 % 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 +7.3 −4.5 %. 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 +3.8 −2.8 % for FGK stars. 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 we consider 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 data with 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. Figure B.1 shows the depth of search of the SHINE survey for the 150 stars in the sample based on different assumptions for 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-2003 evolutionary models. 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-10 au 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.