The intrinsic multiplicity distribution of exoplanets revealed from the radial velocity method

We follow the general framework of Tremaine & Dong (2012) to derive the intrinsic multiplicity distribution from an exoplanet survey with limited sensitivities. We use $F$ for the intrinsic multiplicity vector, with each component $F_{k}$ representing the fraction of stars with $k$ planets in the target parameter space. ¹¹1Throughout this paper the letter $k$ is used for the number of intrinsic planets, regardless of whether they are detected or not, whereas the letter $j$ is used for the number of planets that are actually detected. The observed multiplicity vector is $N$, with each component $N_{j}$ representing the number of stars with $j$ detected planets. With a statistical sample of $N_{\star}$ stars searched in the survey, a relation can be established

or equivalently

Here $\bar{\mathcal{S}}$ is a matrix that captures the survey detection sensitivity, averaged over all stars in the sample. Each entry of the matrix, $\bar{S}_{jk}$, quantifies the (averaged) probability that the survey identifies $j$ planet detections in a typical system of $k$ planets (Tremaine & Dong, 2012).

Under the general assumption that the detections of individual planets around the same host are independent, the sensitivity matrix $\mathcal{S}$ of the $l$-th star in the sample is given as

where $p_{l}$ is the detection probability of the $l$-th star to a typical planet in the pre-defined parameter space $\Delta\mbox{\boldmath$\theta$}$

Here $\theta$ denotes a set of parameters that define the properties of the planet. Note that the evaluation of $p_{l}$ takes into account not only the detection sensitivity of the star, $s_{l}(\mbox{\boldmath$\theta$})$, but also the intrinsic planet distribution of all stars, $f(\mbox{\boldmath$\theta$})$. The mean sensitivity of the survey is thus

The frequency of planets, or the average number of planets per star, is

With the above equations it is straightforward to show that

where the averaged detection probability

That is, the derivation of the frequency (or occurrence rate) of planets is independent of the multiplicity distribution. This is consistent with the original assumption that the detections of individual planets are independent. Additionally, as the planet frequency only depends on the averaged detection probability $\bar{p}$, it provides a practical way to derive the intrinsic planet distribution $f(\mbox{\boldmath$\theta$})$. This last quantity is central to the calculation of the per-star detection probability $p_{l}$ (Equation 4). We will further demonstrate it in Section 3.

Another quantity of interest to the exoplanet statistics is the frequency of planetary systems, or the fraction of stars with (at least one) planets,

Unlike the derivation of $\bar{n}_{\rm p}$, the derivation of $F_{\rm p}$ depends on planet multiplicity.

The ratio between the frequency of planets and the frequency of planetary systems gives the average number of planets per planetary system, or average multiplicity (Zhu & Dong, 2021)

With a given intrinsic multiplicity vector, $F$, and the sensitivity matrix, $\bar{\mathcal{S}}$, one can compute the expected multiplicity distribution, $\bar{\mbox{\boldmath$N$}}\equiv\{\bar{N}_{j}\}$, based on Equation (1). The detections of individual systems are best described as Poisson processes, and thus the likelihood of having the observed multiplicity distribution, $\mbox{\boldmath$N$}\equiv\{N_{j}\}$, is given by

and its logarithm

The constant term is ignored hereafter. The above likelihood is maximized when $\bar{\mbox{\boldmath$N$}}=\mbox{\boldmath$N$}$, leading to the solution

The above analytic solution does not take into account the physical constraints on the intrinsic multiplicities, namely $0\leq F_{k}\leq 1$. Nevertheless, it provides a simply and useful way to check the numerical results and investigate the impact of input parameters on the derived frequencies.

One can also derive/constrain the intrinsic multiplicities in a numerical way. This allows to impose the physical constraints and derive the uncertainties of the intrinsic multiplicities. The data to be fitted by the model is the observed multiplicity distribution, $(N_{0},N_{1},N_{2},\dots)$, and the free parameters of the model are $(F_{1},F_{2},\dots)$. Note that the first component of $F$, $F_{0}$, is not a free parameter of the model, because it can be derived by $F_{0}=1-\sum_{j=1}F_{j}$. Therefore, the number of observational constraints is always larger than the number of free parameters, and thus the model can be constrained.

The observed multiplicity distribution can in principle be extended to infinity, which cannot be handled in practice, so where to truncate it? Theoretically, mass conservation of the planetary material and long-term dynamical stability arguments do not allow a system to host infinite number of planets. However, such arguments do not give deterministic criterion on the maximum number of planets per system, so we rely on the statistical argument and regard it as an issue of model selection. The likelihood of the model is optimized under Equation (15), and the inclusion of higher multiplicities with zero detections no longer enhances it. Therefore, both intrinsic and observed multiplicity distributions are truncated at the highest multiplicities with non-zero detections. ²²2A few criteria, such as the Bayesian Information Criterion (BIC), have been proposed to quantify the superiority of one model over another. However, the commonly used form of BIC assumes that the number of data points substantially exceeds the number of model parameters. This is not true in the present case. We will discuss the robustness of our results against the choice of the maximum multiplicity in the modeling, $k_{\rm max}$, in Section 4.2.

We derive the uncertainties on the intrinsic multiplicities in the following way. We generate a fine grid around the maximum-likelihood solution in Equation (15) and calculate the likelihood value at each grid point. The $n$-$\sigma$ confidence interval is then defined by grid points whose likelihoods are no more than $\Delta\ln{\mathcal{L}}=n^{2}/2$ from the maximum value.

The CLS sample contains 178 RV planet detections around 719 stars (Rosenthal et al., 2021a; Fulton et al., 2021). The CLS stars cover a broad range in stellar properties, including FGKM-type dwarf and giant stars. For this study, we will focus on the Sun-like stars, which are defined as stars with effective temperature, $T_{\rm eff}$, in the range of 4700–6500 K and surface gravity $\log{g}>4.0$. Applying these cuts to the CLS sample, we are left with 474 Sun-like stars. These include 92 known planet hosts, and the total number of CLS planet detections is 142. Figure 1 shows the selected Sun-like stars and planet hosts in comparison to the overall CLS sample.

We reconstruct the CLS completeness curves for the selected stars based on the per-star injection–recovery simulations of Rosenthal et al. (2021a). The averaged completeness curves are shown in Figure 2 together with the CLS planet detections around Sun-like hosts. Our sample has $(382,57,26,5,2,2)$ systems with $(0,1,2,3,4,5)$ planet detections.

In principle, one can use the overall observed multiplicity distribution and the estimated completeness curves to constrain the intrinsic multiplicity vector of planets across the whole mass–semi-major axis parameter space. However, given the limited sensitivity of CLS, such an attempt is determined to be very sensitive to assumptions about the unknown distributions of low-mass, long-period planets. We therefore limit ourselves to certain regions of the parameter space. Specifically, we identify the following types of planets (see Figure 2 for an illustration):

• •

  Close-in giant planets: planets with semi-major axis $a<1\,$AU and minimum mass $m\sin{i}$ in the range 0.3–13$\,M_{\rm J}$;

• •

  Close-in intermediate-mass planets: planets with semi-major axis $a<1\,$AU and minimum mass $m\sin{i}$ in the range 10$\,M_{\oplus}$–0.3$\,M_{\rm J}$;

• •

  Cold giant planets: planets with semi-major axis $a$ in the range 1–10 AU and minimum mass $m\sin{i}$ in the range 0.3–13$\,M_{\rm J}$;

We also study the combinations of the above planet types, namely close-in giant and intermediate-mass planets and all giant planets (including close-in and cold giant planets). Additionally, for the purpose of future RV observations, we also study the planet distribution in the RV semi-amplitude $K$ vs. semi-major axis plane and derive the multiplicity distributions for planets with $K>3\,{\rm m\,s^{-1}}$ and $a<10\,$au as well as with $K>1\,{\rm m\,s^{-1}}$ and $a<1\,$au. The observed multiplicity distributions for different types of planets are summarized in Table 1.

The method outlined in Section 2 to derive the intrinsic planet multiplicity has assumed that the properties (e.g., orbital periods, mass, etc) and detections of individual planets around the same host are mutually independent. Given data with high enough quality, the detections of multiple planets in the same system are indeed independent. In practice, however, the detectability of a new planet is dependent on the properties of the detected planets. One well-known example is the degeneracy between one eccentric planet and two near-circular planets in 2:1 resonant orbits (e.g., Anglada-Escudé et al., 2010; Wittenmyer et al., 2019). The recent analysis by Boisvert et al. (2018) found that about 25% of published RV systems might suffer from this degeneracy. This fraction may be taken as an upper limit on the impact of the degeneracy in the CLS sample, given that CLS on average has more data points and covers longer durations. For systems that may actually suffer from this degeneracy, the hidden planets are usually not massive enough, ³³3Specifically, the hidden planet is lighter than the observed one by a factor of $\sim e/2^{1/3}$, where $e$ is the orbital eccentricity of the reported planet. and thus our main conclusion about the giant planet multiplicity will not be affected. Future studies and observations will be needed to understand the impact of the degeneracy on the derived intrinsic multiplicity of low-mass planets.

One indirect way to check the validity of the assumption about independence is to compare the parameters of planets from systems with only one planet detection and at least two planet detections. Such a comparison is shown in Figure 3. The two-sample Kolmogorov-Smirnov (KS) tests yield $p$-values on the order of unity, suggesting that the two distributions under comparison can be considered to be drawn from the same underlying distribution. In other words, the planetary parameters are statistically similar regardless of the number of planet detections. Therefore, it is reasonable to assume that the planet properties and detections are largely independent in multi-planet systems. See also Tremaine & Dong (2012) for further discussion on this separability issue.

Detected planets in the Kepler multi-planet systems appear to have similar properties (e.g., Lissauer et al., 2011; Ciardi et al., 2013). Whether or not this is of astrophysical origin is still under debate (e.g., Millholland et al. 2017; Weiss et al. 2018; He et al. 2019; Zhu 2020; Murchikova & Tremaine 2020), and we refer to Zhu & Dong (2021) for further discussions about this issue. Regardless of the origin, it has marginal impact on the method and results of the current work for two reasons. First, even if the observed intra-system similarity in Kepler systems is astrophysical, such a correlation will be much weaker in typical RV systems, which have on average more massive and fewer planets. Second, as explained earlier in the paper, we study the multiplicity distribution of specific planet types whose physical properties are already restricted to a relatively narrow range (e.g., all giant planets have very similar sizes). In other words, the hypothetical planets in the hypothetical multi-planet systems implied in our method already have similar properties, even though they are generated in an independent way. The primary parameter that is affected by assumptions like the intra-system similarity is the detection probability $p$. We will discuss in Section 4.1 how our results will be changed with different average detection probabilities.

We first derive the planet frequency, $\bar{n}_{\rm p}$, over the chosen parameter space, in order to determine the mean detection probability $\bar{p}$ as well as the detection probabilities $\{p_{l}\}$ for individual stars.

We follow the general approach of Zhu & Dong (2021) to derive the planet frequencies. The parameter space is first divided into a combination of grids, with each grid covering an area that is small compared to the overall parameter space but large compared to the positional uncertainties of enclosed planets. Within each grid, the planet distribution function can be treated as a constant, and the average detection probability is given by the mean efficiency of the survey detection. The posterior distribution of the planet frequency within each grid is given by a Gamma distribution. For any grid with at least two planet detections, a log-flat prior on the planet frequency is adopted to determine the median and the 68% confidence interval. Otherwise, the 95% upper limit is derived under a linear prior. The results are illustrated in Figures 4 and 5 for planets in the mass vs. semi-major axis plane and the RV semi-amplitude vs. semi-major axis plane, respectively.

For our definition of close-in giant, close-in intermediate-mass, and cold giant planets, the planet frequency, $\bar{n}_{\rm p}$, is determined to be $0.076\pm 0.013$, $0.21\pm 0.05$, and $0.22\pm 0.03$, respectively. In the $K$ vs. $a$ plane, the planet frequency is $0.42\pm 0.04$ for planets with $K>3\,{\rm m\,s^{-1}}$ and $a<10\,$au and $0.58\pm 0.09$ for planets with $K>1\,{\rm m\,s^{-1}}$ and $a<1\,$au. The average detection probabilities for these planet classes can be derived according to Equation (7). Similarly, the average detection probability for planets defined in the joint parameter space (i.e., close-in CLS planets and all giant planets) as well as those in the RV semi-amplitude vs. semi-major axis plane can also be derived. The median values of $p$ for all planet classes are included in Table 1.

In order to derive of the intrinsic multiplicity distribution, we need to calculate the per-star detection probability, $p_{l}$, which is the detection probability of a typical planet around the $l$-th star in the sample. Making use of the injection–recovery simulations, this is given by

Here $f(\mbox{\boldmath$\theta$})$ is the rate map in the specified parameter space (i.e., Figure 4 for $m\sin{i}$ vs. $a$ and Figure 5 for $K$ vs. $a$), and $s_{l}(\mbox{\boldmath$\theta$}_{j^{\prime}})$ evaluates whether or not the $j^{\prime}$-th injected planet around the $l$-th star is detected in the simulation. The summations are performed over all simulated planets that are within the pre-defined parameter space. We show in Figure 6 the distributions of the per-star detection probabilities for the five planet types considered here. As it indicates, this detection probability parameter varies substantially among all stars in the sample, primarily due to the different numbers and quality levels of RV observations for individual stars.

These per-star detection probabilities enable the calculation of the mean sensitivity matrix, $\bar{\mathcal{S}}$, a key ingredient in the derivation of the intrinsic multiplicity distribution.

We model the observed multiplicity distribution in the CLS Sun-like star sample following the method detailed in Section 2. The results are summarized in Table 2 and illustrated in Figure 7. Below we discuss the individual planet classes in details.

The per-star detection probabilities, $\{p_{l}\}$, are crucial in the derivation of the intrinsic multiplicity distribution as well as the overall frequencies. As defined by Equation (4), this per-star detection probability takes into account both the intrinsic planet distribution as well as the survey completeness. While the latter is (or can be) known fairly precisely, the former is usually not. As a result, the detection probability may not be so well determined.

Our derivation of the intrinsic multiplicity distribution makes use of the per-star detection probabilities, and thus it is less susceptible to the variation in the intrinsic planet distribution. Nevertheless, we investigate the impact of the mean detection probability parameter, $\bar{p}$, on the derived frequencies, as the intrinsic multiplicity distribution can be approximated as

The latter step makes use of the property of the sensitivity matrix (Tremaine & Dong, 2012). Figure 8 shows the derived frequencies—frequency of planetary systems, frequency of planets, and average multiplicity—as functions of the average detection probability $\bar{p}$ for the “all giant planet” case. We may convert the uncertainty in the frequency of planets, $\bar{n}_{\rm p}$, into the uncertainty of the average detection probability, given the direct connection between these two quantities (Equation 7). This gives $\bar{p}_{\rm all~{}giants}=0.65\pm 0.11$. The variation of $\bar{p}$ at such a level will not lead to substantial changes in the other two frequencies, both of which are rather insensitive to the variation of $\bar{p}$.

In our maximum-likelihood analysis, the multiplicity distribution is truncated at the highest multiplicity with non-zero detections. As explained in Section 2.2, this is because the inclusion of higher multiplicities in the data and the modeling introduces additional degrees of freedom but no benefit in the goodness of fit. Here we use the case of “all giant planet” to further demonstrate this point.

The CLS Sun-like sample has zero systems with more than four giant planets inside $10\,$au, and thus we have adopted the maximum multiplicity $k_{\rm max}=4$ in our maximum-likelihood analysis. The results have been reported in the previous section, and the likelihood contours are shown in Figure 9. For comparisons, we also perform the analysis with $k_{\rm max}=5$. The resulting likelihood contours are also illustrated in Figure 9. The inclusion of the higher multiplicity with zero planet detections does not change the maximum-likelihood solution or the uncertainties on all model parameters except for $F_{4}$, of which the lower limits are slightly modified. This is because the introduction of $F_{5}$ permits $F_{4}=0$, which is otherwise not allowed with $k_{\rm max}=4$. The integrated fraction of planetary systems and the average multiplicity remain largely unaffected, too.

Regarding the frequencies of giant planets around Sun-like stars, we find

Here “HJ”, “WJ”, and “CJ” denote hot Jupiters (giant planets within 0.1 au), warm Jupiters (giant planets between 0.1–1 au), and cold Jupiters (giant planets between 1–10 au), respectively. The HJ rate we derive is higher than previous RV studies, which collectively points to a value of $\sim 1\%$ (e.g., Cumming et al., 2008; Mayor et al., 2011; Wittenmyer et al., 2020). However, given the large statistical uncertainties and the possibly different stellar distributions (e.g., metallicity, binarity) in these studies, we do not investigate further the cause of this discrepancy.

Regarding the WJ and CJ rates, a very direct comparison is with Fulton et al. (2021), who used the whole CLS sample and derived the frequency of planets around FGKM stars. ⁵⁵5Note the different definition of a “giant planet” used in Fulton et al. (2021) and here. Fulton et al. (2021) set the lower limit of a “giant planet” to be $30\,M_{\oplus}$. Adopting their best-fit broken power-law model, we find $\sim 0.05$ and $\sim 0.14$ for the frequencies of WJs and CJs, respectively. While the WJ rate is consistent with ours, their CJ rate is substantially lower. This is probably because Fulton et al. (2021) included M-type stars in their analysis. In the CLS sample, there are only 13 cold giant planets found around 245 M-type stars, whereas the numbers are 57 and 474 for Sun-like stars. Our rates of WJs and CJs are also in general agreement with studies using different samples (e.g., Mayor et al., 2011; Fernandes et al., 2019; Wittenmyer et al., 2020). For example, the re-analysis of the Mayor et al. (2011) sample by Fernandes et al. (2019) found $\sim 0.04$ and $\sim 0.17$ for frequencies of WJs and CJs, respectively. These values are estimated based on the asymmetric broken power-law model used in Fernandes et al. (2019).

The relative frequencies between HJ and WJ can provide useful insights into the formation of these close-in giant planets (see Dawson & Johnson 2018 and references therein). In the CLS Sun-like star sample, the overall rate of WJ is higher than the overall rate of HJ, and the rate of WJ per log interval is within a factor of two of that of HJ. Both features pose challenges to the theoretical models of HJ formation.

One approach to derive the frequency of planetary systems in RV surveys, which has been frequently used in the literature, is to include only one out of the multiple planet detections in each system and regard the resulting frequency of planets as the frequency of planetary systems (e.g., Cumming et al., 2008; Mayor et al., 2011). In the mathematical framework of the present work, this is

Here $\tilde{p}$ is the average detection probability derived from the planet sample that includes at most one planet detection for each star. It will not be too different from the value of $\bar{p}$, given the dominating number of single-detections in typical RV samples.

For our approach, which takes the planet multiplicity fully into account, a reasonable approximation of the frequency of planetary systems is

Here the correction factor, $p_{m}$, quantifies the probability of detecting at least one planet in a system with $\bar{m}_{\rm p}$ planets

The estimator given by Equation (17) yields unbiased rates only when the average multiplicity $\bar{m}_{\rm p}$ approaches unity. We provide in the last two rows of Table 2 the two frequency estimators for our chosen planet types. In all cases, our new estimator, $F_{\rm p}^{\rm multi}$, yields better agreement with the true frequency.

The estimator given by Equation (18) has broader implications to other statistical analyses that involve planet multiplicity issues. For Sun-like stars, the sensitivity region of the Kepler transit mission is almost fully enclosed in the box-shaped parameter space defined by a lower radius limit of $\sim 1\,R_{\oplus}$ and an outer period limit of $\sim 1\,$yr. According to the estimated planet frequencies in this box-shaped region (e.g., Zhu & Dong, 2021), the probability of having any one planet inside the Kepler sensitivity region is $\bar{p}_{\rm kep}\approx 0.8$. For Kepler-like planetary systems with an average multiplicity of three (Zhu et al., 2018), the chance to see at least one planet out of the box-shaped parameter space to be inside the Kepler sensitivity region is therefore $1-(1-\bar{p}_{\rm kep})^{3}=0.99$. In other words, the fraction of Sun-like stars with planets of radii above $R_{\oplus}$ and periods within 1 yr is almost the same as the fraction of Sun-like stars with at least one planet inside the Kepler sensitivity range. For comparison, the sensitivity correction on the planet-by-planet basis would lead to an increase in the frequency of planetary systems by nearly a factor of two (e.g., Yang et al., 2020).

In addition to providing an unbiased estimate of the total frequency of planetary systems, the approach undertaken in the present work also yields the intrinsic multiplicity distribution $F$, which is more informative to theoretical models as well as other statistical studies. For example, for giant planets within $3\,$au, based on which Cumming et al. (2008) derived a total frequency of $10.5\%$, our approach yields the following total rates based on the CLS Sun-like sample

as well as the intrinsic multiplicity distribution

The overall rate ($F_{\rm p}^{\rm C08}$) is higher than Cumming et al. (2008) probably because of the difference in the stellar samples.

By including in the analysis the intrinsic multiplicity distribution, our method does not require the presences of multiple planets in the same system to be independent (although the properties of them should be, largely, independent). Therefore, it allows us to investigate whether or not the existence of a certain type of planet depends on the existence of some other type of planet in the same system. Considering the relatively large uncertainties of the derived frequencies, here we focus on confirming previously found (anti-)correlations rather than reporting tentative evidence for new ones.

Under the assumption that the presences of close-in and cold giant planets are not correlated, the frequency of planetary systems with any giant planets within 10 au should be $1-(1-7.2\%)\times(1-17\%)=23\%$. This is higher than what is given by the maximum likelihood analysis over the joint parameter space ($19.2\%$, see Table 2). Therefore, the assumption is rejected/disfavored and the close-in and cold giant planets are positively correlated. This is consistent with previous findings that close-in giant planets frequently have distant companions (e.g., Knutson et al., 2014; Bryan et al., 2016). Similarly, our results also suggest positive correlations between close-in planets (close-in giant and intermediate-mass planets) and cold giant planets, which also confirms previous findings (e.g., Zhu & Wu, 2018; Bryan et al., 2019; Rosenthal et al., 2021b).

It is noticed that the fraction of stars with close-in giant and intermediate-mass planets equals the summation of frequencies of “close-in giant” and “close-in intermediate-mass” planets. Indeed, under the assumption that the presences of these latter two classes of planets are uncorrelated, one would expect that the frequency of close-in CLS planetary system should be $1-(1-0.072)\times(1-0.15)=21\%$, which is consistent with the fraction of planets in the joint parameter space ($21\pm 4\%$). In other words, the two planet populations—close-in giant planets and close-in intermediate-mass planets—are uncorrelated. This may be explained by the combination of two different features of close-in planets. On the one hand, the majority of HJs and some WJs are lonely (e.g., Steffen et al., 2012; Hord et al., 2021). On the other hand, at least half of WJs do have nearby, small companions (Huang et al., 2016).