About 30% of Sun-like Stars Have Kepler-like Planetary Systems: A Study of Their Intrinsic Architecture

To facilitate further discussions, we introduce a few mathematical notations here. Following Tremaine & Dong (2012), we use a matrix G to quantify the detection probability of planetary systems in transit surveys. Each element g_jk denotes the probability that one k-planet system is seen to have j (j ≥ 1) tranets. Thus we have

$$\begin{eqnarray}&&{g}_{{jk}}=0,\quad \mathrm{if}\ j\gt k.\end{eqnarray} \tag{ 1 }$$

If the parameter space that is of interest can fit in at most K planets, then G should be a K × K upper-triangular matrix.⁶

We use symbol ${ \mathcal N }$ for the number of stars in the sample and the vector N = (N₁, ..., N_K) for the observed transit multiplicity function. With the intrinsic multiplicity vector F, the expectation of the transit multiplicity function can be given as

$$\begin{eqnarray}&&{\bar{N}}_{j}=\displaystyle \sum _{k=j}^{K}{g}_{{jk}}{ \mathcal N }{f}_{k}=\displaystyle \sum _{k=1}^{K}{g}_{{jk}}{ \mathcal N }{f}_{k},\quad \mathrm{or}\,\bar{{\boldsymbol{N}}}={\boldsymbol{G}}\cdot ({ \mathcal N }{\boldsymbol{F}}).\end{eqnarray} \tag{ 2 }$$

The second equality in the summation form has used the fact that g_jk = 0 if k < j. The fraction of stars with planets, F_p, and the average number of planets per star, ${\bar{n}}_{{\rm{p}}}$, are given by

$$\begin{eqnarray}&&{F}_{{\rm{p}}}=\displaystyle \sum _{k=1}^{K}{f}_{k},\quad {\bar{n}}_{{\rm{p}}}=\displaystyle \sum _{k=1}^{K}{{kf}}_{k},\end{eqnarray} \tag{ 3 }$$

respectively. The ratio of these two quantities, ${\bar{n}}_{{\rm{p}}}/{F}_{{\rm{p}}}$, gives the average number of planets per planetary system, which we call the average multiplicity.

We also introduce the matrix T to quantify the detection probability of TTVs. Each element t_jk represents the probability that one k-planet system has j tranets and at least one of them shows detectable TTV signals. The TTV multiplicity function is given by M = (M₁, ..., M_K) and the expectation of this is given by

$$\begin{eqnarray}&&{\bar{M}}_{j}=\displaystyle \sum _{k=j}^{K}{t}_{{jk}}{ \mathcal N }{f}_{k},\mathrm{or}\,\bar{{\boldsymbol{M}}}={\boldsymbol{T}}\cdot ({ \mathcal N }{\boldsymbol{F}}).\end{eqnarray} \tag{ 4 }$$

In our model, whether a planet transits or not is determined by the transit parameter  ≡ R_⋆/a and its orbital inclination I_p. We do not take into account the minor impact of the planet size. Below we describe the distributions of and I_p. We also describe the criteria we use in generating multi-planet systems and detecting TTV signals.

3.2.1. Distribution of Transit Parameters

Following Tremaine & Dong (2012), we model the distribution of transit parameter as

$$\begin{eqnarray}&&\displaystyle \frac{{dN}}{d\mathrm{ln}\epsilon }\propto \displaystyle \frac{{(\epsilon /{\epsilon }_{0})}^{a}}{1+{(\epsilon /{\epsilon }_{0})}^{b}},\end{eqnarray} \tag{ 5 }$$

which is essentially a broken power law but with smooth transition at ₀. Instead of using the values for from transit modelings that are not well constrained for some tranets, we re-compute them based on the orbital periods and LAMOST stellar parameters (R_⋆ and M_⋆). In this way this parameter is better constrained and its lower and upper boundaries are compatible with our sample selection criteria (Section 2). This reconstructed distribution using all tranets in our sample is shown as black dots in Figure 4, in which we also show the distribution from only transit singles for a reference.

Zoom In Zoom Out Reset image size

We then model this distribution with the smoothed broken power-law form (Equation (5)). After correcting for the geometric transit probability (∝), we determine the underlying distribution to be

$$\begin{eqnarray}&&\displaystyle \frac{{dN}}{d\mathrm{ln}\epsilon }=0.36\displaystyle \frac{{(\epsilon /{\epsilon }_{0})}^{0.04}}{1+{(\epsilon /{\epsilon }_{0})}^{3.18}}\quad ({\epsilon }_{0}=0.074)\end{eqnarray} \tag{ 6 }$$

for 0.004 <  < 0.6, and zero elsewhere. This yields a logarithmically flat distribution below ₀, a result in agreement with previous studies (e.g., Dong & Zhu 2013; Petigura et al. 2013).

3.2.2. Distribution of Planetary Inclinations

For multi-planet systems, the planetary inclination relative to the observer, I_p, depends on the inclination of the system invariable plane, I, the planet inclination with respect to this invariable plane, i, and a nuisance parameter ϕ (i.e., the phase angle)

$$\begin{eqnarray}&&\cos {I}_{{\rm{p}}}=\cos I\cos i-\sin I\sin i\cos \phi .\end{eqnarray} \tag{ 7 }$$

The distribution of I is isotropic (∝sin I for 0° ≤ I ≤ 180°) and the distribution of ϕ is random between 0° and 360°. The inclination i quantifies the flatness of the multi-planet system and we assume that it is related to the number of planets in the system k. We parameterize this dependence as a power law between the dispersion of i (or more accurately, sin i) and k, and choose the normalization at k = 5:

$$\begin{eqnarray}&&{\sigma }_{i,k}\equiv \sqrt{\langle {\sin }^{2}i\rangle }={\sigma }_{i,5}{\left(\displaystyle \frac{k}{5}\right)}^{\alpha }.\end{eqnarray} \tag{ 8 }$$

It is written in this form, so that the normalization factor, ${\sigma }_{i,5}$, can be determined separately from the distribution of transit duration ratios of planet pairs in five-planet systems (Section 3.3). For such high-multiple planetary systems, the observed multiplicity very likely reflects the intrinsic multiplicity. We decide to use k = 5 for the normalization term for two reasons. First, there are only a few k ≥ 6 planetary systems found by Kepler, the number being so small that the mutual inclination dispersion cannot be well constrained. Second, although there are more four-planet systems than five-planet systems, the fraction of contamination from intrinsically higher multiplicities is also larger for four-planet systems than for five-planet systems. The power-law index α quantifies the steepness of this inclination dispersion function and can be constrained from the transit and TTV multiplicity functions.

With this inclination dispersion σ_i,k, the planetary inclination with respect to the invariable plane is then modeled as a Fisher (1953) distribution (Fabrycky & Winn 2009; Tremaine & Dong 2012),

$$\begin{eqnarray}&&P(i| {\kappa }_{k})=\displaystyle \frac{{\kappa }_{k}\sin i}{2\sinh {\kappa }_{k}}{e}^{{\kappa }_{k}\cos i}.\end{eqnarray} \tag{ 9 }$$

The parameter κ_k is related to the dispersion parameter σ_i,k via

$$\begin{eqnarray}&&{\sigma }_{i,k}^{2}=\langle {\sin }^{2}i\rangle =\displaystyle \frac{2}{{\kappa }_{k}}\left(\coth {\kappa }_{k}-\displaystyle \frac{1}{{\kappa }_{k}}\right).\end{eqnarray} \tag{ 10 }$$

This Fisher distribution provides a smooth transition from an isotropic distribution (κ_n ≪ 1) to a Rayleigh distribution (κ_n ≫ 1). The latter one is commonly used for compact multi-planet systems (e.g., Fabrycky et al. 2014).

The inclination dispersion σ_i,k, by its definition given by Equation (8), has a maximum value of $\sqrt{2/3}$, which can be achieved only when the distribution of i becomes isotropic. For any given ${\sigma }_{i,5}$, the upper bound on the inclination dispersion therefore sets a limit on α (and vice versa).

3.2.3. Stability Criterion

We describe the stability criterion used in generating multi-planet systems. For intrinsic multiples (k ≥ 2), we inject planets one by one. The transit parameter of the first planet is randomly drawn from the distribution specified by Equation (6), and then the orbital period is derived via

$$\begin{eqnarray}&&P={\left(\displaystyle \frac{{R}_{\odot }}{\mathrm{au}}\right)}^{3/2}{\epsilon }^{-3/2}{\left(\displaystyle \frac{{\rho }_{\star }}{{\rho }_{\odot }}\right)}^{-1/2}\mathrm{year}.\end{eqnarray} \tag{ 11 }$$

Throughout our simulations, we fix ρ_⋆ = ρ_⊙. This is because the distribution of has absorbed the variance of ρ_⋆ and therefore there is no need to assume a separate distribution for ρ_⋆. For any additional planet, the transit parameter and the orbital period are randomly assigned in a similar way but with the restriction that the new planet must be far away from any previously injected planets such that the system remains dynamically stable. For the latter, we use the Deck et al. (2013) stability criterion, which requires that for any given planet pair,

$$\begin{eqnarray}&&\displaystyle \frac{{P}_{\mathrm{out}}}{{P}_{\mathrm{in}}}\gt 1+2.2{q}^{2/7}=1.16.\end{eqnarray} \tag{ 12 }$$

The critical value is derived by assuming a characteristic planet-to-star mass ratio q = 10⁻⁴. The actual choice of this characteristic value has a very marginal impact on the modeling output, primarily because of the weak dependence on q. Furthermore, even though the chosen q value is larger than the typical planet-to-star mass ratio (∼10⁻⁵) of Kepler planets, it is extremely rare to have a period ratio below 1.3 in actual observations (Lissauer et al. 2011; Fabrycky et al. 2014).

3.2.4. TTV Detection Criteria

We acquire external constraints on the normalization factor ${\sigma }_{i,5}$ of the inclination dispersion relation (Equation (8)). This is necessary because otherwise we would end up with a strong correlation between ${\sigma }_{i,5}$ and α.

Because we are constraining ${\sigma }_{i,5}$ separately, we are not limited to the planetary systems in our current sample. In fact, our sample only contains 3 five-tranet systems and these provide $3\times {C}_{5}^{2}=30$ planet pairs. Instead, we find 15 five-tranet systems whose hosts are Sun-like stars from the California Kepler Survey (Petigura et al. 2017), which give us 150 planet pairs to constrain the inclination dispersion ${\sigma }_{i,5}$.

We adopt a similar approach as Fabrycky et al. (2014) to constrain the inclination dispersion. Each observed planet pair gives a quantity (Steffen et al. 2010)

$$\begin{eqnarray}&&\xi \equiv \displaystyle \frac{{T}_{\mathrm{dur},\mathrm{in}}/{P}_{\mathrm{in}}^{1/3}}{{T}_{\mathrm{dur},\mathrm{out}}/{P}_{\mathrm{out}}^{1/3}},\end{eqnarray} \tag{ 16 }$$

where T_dur is the transit duration (from first to fourth contact) and P is the orbital period. The subscripts "in" and "out" denote the values of the inner and outer planets, respectively. With Kepler's third law, one can easily see that this parameter ξ is essentially the ratio of orbital-velocity normalized transit durations, which is most sensitive to the inclination i and is marginally dependent on other parameters such as the orbit eccentricity e (Fabrycky et al. 2014). We can also write ξ in terms of the planet-to-star radius ratio r and the transit impact parameter b

$$\begin{eqnarray}&&\xi =\sqrt{\displaystyle \frac{{(1+{r}_{\mathrm{in}})}^{2}-{b}_{\mathrm{in}}^{2}}{{(1+{r}_{\mathrm{out}})}^{2}-{b}_{\mathrm{out}}^{2}}}.\end{eqnarray} \tag{ 17 }$$

Because both T_dur and P are much better measured than r or b in observations, the expression given by Equation (16) is therefore used in constructing the distribution of ξ from the data. For the 15 five-planet systems, we use the values of T_dur and P from the the most recent Kepler data release (DR25; Thompson et al. 2018). We then compute the ξ values for individual planet pairs and show the distribution as the black histogram on the left panel of Figure 5.

Zoom In Zoom Out Reset image size

We then model this ξ distribution and attempt to constrain the inclination dispersion ${\sigma }_{i,5}$. We ignore the dependence of ξ on the eccentricity e and simply assume circular orbits for all planets in such five-planet systems. While simplifying the modeling, this still remains a very reasonable assumption. First, planets in such high-multiple and stable systems are not expected to have large eccentricities. Second, Fabrycky et al. (2014) showed that this ξ parameter alone could not constrain e very well, and that the correlation between inclination dispersion and eccentricity is weak.

For a given value of ${\sigma }_{i,5}$, we produce the ξ distribution following the method of Fabrycky et al. (2014). First, we randomly draw an impact parameter for the outer planet, b_out,sim, uniformly from the range [0, b_out,max], where b_out,max is the impact parameter the planet would require in order that the total S/N ($\sqrt{{b}_{\mathrm{out},\max }/{b}_{\mathrm{out}}}$ times the actual S/N) would be equal to the S/N threshold (7.1). Then we randomly draw b_in,sim from a normal distribution with mean ${b}_{\mathrm{out},\mathrm{sim}}{({P}_{\mathrm{in}}/{P}_{\mathrm{out}})}^{2/3}$ and dispersion ${\sigma }_{i,5}$/(d_in + r_in), where d ≡ R_⋆/a is the stellar radius scaled to the planet–star separation and is taken as the value from DR25. Such a normal distribution in impact parameters reproduces a Rayleigh distribution with dispersion ${\sigma }_{i,5}$ in inclinations. If $| {b}_{\mathrm{in},\mathrm{sim}}| \leqslant {b}_{\mathrm{in},\max }$, we consider this simulated planet to be detectable. Otherwise we repeat the previous step until the above condition is met. Once a simulated planet pair is generated, we compute ξ using Equation (17). We repeat the whole process and generate 250 simulated pairs for each planet pair in the sample. The ξ distribution for the given ${\sigma }_{i,5}$ is then generated from the ensemble of all simulated pairs.

We then compute the likelihood for the actual ξ distribution to be drawn from the simulated one. We note that this is more accurate than the approach of Fabrycky et al. (2014), which used the p value of Kolmogorov–Smirnov test. We compute the likelihood as

$$\begin{eqnarray}&&{ \mathcal L }=\displaystyle \prod _{j=1}^{N}\int {P}_{\mathrm{sim}}(\mathrm{ln}\xi )\exp \left[-\displaystyle \frac{{(\mathrm{ln}{\xi }_{j}-\mathrm{ln}\xi )}^{2}}{2{\sigma }_{\mathrm{ln}\xi ,j}^{2}}\right]d\,\mathrm{ln}\,\xi ,\end{eqnarray} \tag{ 18 }$$

where P_sim(ln ξ) is the probability distribution of simulated ξ in logarithmic space, ξ_j is the observed value of ξ of the jth planet pair, and σ_lnξ,j is the fractional uncertainty of ξ_i. We compute σ_lnξ,j based on the measured uncertainties on T_dur,in and T_dur,out.

We repeat the Monte Carlo simulation and compute the likelihood ${ \mathcal L }$ for 100 values of ${\sigma }_{i,5}$ equally spaced from 01 to 40, and show the results on the right panel of Figure 5. As the scatter plot shows, the likelihood reaches its maximum around ${\sigma }_{i,5}$ = 08. The simulated ξ distribution for this best-fit ${\sigma }_{i,5}$ is also shown on the left panel of Figure 5. To find the different confidence levels of ${\sigma }_{i,5}$, we use a spline function to smooth the $\mathrm{ln}\,({ \mathcal L }/{{ \mathcal L }}_{\max })$ versus ${\sigma }_{i,5}$ scatter plot, and find the 1σ, 2σ, and 3σ confidence intervals, defined as $\mathrm{ln}({ \mathcal L }/{{ \mathcal L }}_{\max })\geqslant -{n}^{2}/2$ to be 065–096, 053–116, and 042–144 for n = 1, 2, and 3, respectively.

For the subsequent modeling of the power-law index α and intrinsic multiplicity vector F, we will only consider values of ${\sigma }_{i,5}$ from the 3σ confidence interval. The smoothed $\mathrm{ln}\,({ \mathcal L }/{{ \mathcal L }}_{\max })$ is also used as our prior on ${\sigma }_{i,5}$ unless specified otherwise.

We use Monte Carlo simulations to compute the G and T matrices for a grid of ${\sigma }_{i,5}$ and α. We note that Tremaine & Dong (2012) provided an analytical formalism to compute the G matrix. However, the time-limiting factor is always the computation of the T matrix, which may not be done in the analytical way. We compute the T matrix within a Monte Carlo, and this automatically produces the G matrix.

For the intrinsic singles (k = 1), there is no TTV signal (t₁₁ = 0), and only one parameter (g₁₁) needs to be calculated. This can be done analytically (Tremaine & Dong 2012): ${g}_{11}=\langle \epsilon \rangle =0.03$.

For intrinsic multiples (k ≥ 2), we first randomly draw the inclination of the invariable plane, I, from an isotropic distribution and then inject planets one by one with the stability criterion given in Section 3.2.3 imposed. For each planet, we assign randomly a phase parameter, ϕ, from the uniform distribution and the planet inclination with respect to the invariable plane, i, from the Fisher distribution (Equation (9)), whose parameter κ_k is determined by Equations (8) and (10). The inclination of the planet with respect to the line of sight is then given by Equation (7). Finally, whether a planet is a tranet or not is determined by the transit criterion

$$\begin{eqnarray}&&\epsilon \gt \cos {I}_{{\rm{p}}}.\end{eqnarray} \tag{ 19 }$$

For any system with at least one tranet, we invoke the TTV criteria in Section 3.2.4 to determine whether there is detectable TTV signals.

With the above procedures, we are able to generate a k-planet system, compute the number of tranets, and determine whether any of the tranets shows detectable TTV signals. Given that we only have up to six tranets in one system, our simulations run up to six-planet systems. The impact of higher multiples to our results will be discussed in Section 5.1. We repeat the whole process and generate a large number of planetary systems, until each element in the T matrix is determined to <2% precision.

For a given intrinsic multiplicity vector F and matrices G and T, the probability to see the transit and TTV multiplicity functions as shown in Figure 2 is given by

$$\begin{eqnarray}&&{ \mathcal L }=\displaystyle \prod _{k=0}^{K}\displaystyle \frac{{\bar{N}}_{k}^{{N}_{k}}\exp (-{\bar{N}}_{k})}{{N}_{k}!}\times \displaystyle \prod _{k=1}^{K}\displaystyle \frac{{\bar{M}}_{k}^{{M}_{k}}\exp (-{\bar{M}}_{k})}{{M}_{k}!},\end{eqnarray} \tag{ 20 }$$

where ${\bar{N}}_{k}$ and ${\bar{M}}_{k}$ are given by Equations (2) and (4), respectively. Note that here we have extended the transit multiplicity function down to k = 0, with ${f}_{0}=1-{\sum }_{k\geqslant 1}\,{f}_{k}$. In practice, we ignore the constants in the log of the above likelihood and try to maximize the following quantity to find the best model parameters (${\sigma }_{i,5}$, α, and F):

$$\begin{eqnarray}&&\mathrm{ln}{ \mathcal L }=\displaystyle \sum _{k=0}^{K}({N}_{k}\mathrm{ln}{\bar{N}}_{k}-{\bar{N}}_{k})+\displaystyle \sum _{k=1}^{K}({M}_{k}\mathrm{ln}{\bar{M}}_{k}-{\bar{M}}_{k}).\end{eqnarray} \tag{ 21 }$$

Given the large number of dimensions, we choose the Markov Chain Monte Carlo (MCMC) algorithm that is implemented in emcee (Foreman-Mackey et al. 2013) as the optimization method. This way we also obtain the posterior distributions of F_k (k = 1, ..., 6) for given sets of ${\sigma }_{i,5}$ and α.

We find the maximum likelihood values for each grid point in the (${\sigma }_{i,5}$, α) plane following the procedure detailed in the previous section. After imposing the prior probability on ${\sigma }_{i,5}$ from the transit duration ratios (Figure 5), we can determine the nσ contours, where n (=1, 2, 3) refers to the number of σ and the contour is defined by $\mathrm{ln}({ \mathcal L }/{{ \mathcal L }}_{\max })=-{n}^{2}/2$. Consequently, we can construct the posterior distribution of the power-law index α in a similar way. These results are shown in Figure 6.

Zoom In Zoom Out Reset image size

Figure 6 indicates that the power-law index α is well constrained to be close to its lower bound −4, which is given by the normalization factor ${\sigma }_{i,5}$ ≈ 08 and ${\sigma }_{i,2}\leqslant \sqrt{2/3}$ (Section 3.2). Specifically, α < −3 at the 1σ level, α < −2 at the 2σ level, and α ≥ 0 can be securely excluded. Therefore, the more planets a system has, the smaller the planetary inclination dispersion is, and this dispersion is a steep function of the intrinsic multiplicity.

To understand what information is driving the constraint on α, we show in Figure 7 the best-fit models with fixed values of α. As this figure shows, a steep inclination dispersion function is required primarily because of the large number of TTV singles, which contribute nearly half of systems with TTVs. Although our sample only contains 23 TTV singles, the ratio between numbers of TTV singles and TTV multiples remains essentially the same even if a much larger TTV catalog (Holczer et al. 2016) or a different TTV catalog (Ofir et al. 2018) is used. Therefore, our constraints on the power-law index parameter is robust.

Zoom In Zoom Out Reset image size

Figure 7 seems to suggest that there may be more TTV singles than even the steepest model curve (α = −4) can account for. However, with our current sample size it is not statistically significant that would require a more complicated model than our current one.

We stack the Markov chains from all MCMC runs and discard entries that are more than 3σ away (Δχ² > 9) from the best one. With this combined Markov chain, we can then investigate the constraints on individual components of the intrinsic multiplicity vector F.⁸

The left panel of Figure 8 shows the full posterior distributions of individual components f_k, which is the fraction of Sun-like stars with k (1 ≤ k ≤ 6) Kepler-like planets. As we can see, the majority of the components are not constrained very well because of the strong degeneracies between neighboring components. However, it is notable that the first component, f₁, is consistent with zero. That is, there can be effectively zero intrinsic singles and nearly all the transit singles are in fact intrinsic multiples with the additional planets non-transiting. This is a result we have discussed qualitatively in Section 1.

Zoom In Zoom Out Reset image size

For future practical usage (such as to predict the yield of future transit missions), we nevertheless report measurements and uncertainties of individual components f_k. This is done by taking the median, 16%–84%, and 5%–95% of the posterior distributions. The result is shown in the right panel of Figure 8. Note that this plot seems to suggest a sharp drop in occurrence rate from low multiples (k ≤ 4) to high multiples (k ≥ 5), but this feature is artificial and comes purely from the way these values are derived.

Although the individual components of the intrinsic multiplicity vector F are not well constrained, the overall occurrence rates, meaning the total fraction of stars with planets F_p and the average number of planets per star ${\bar{n}}_{{\rm{p}}}$, are found to be well constrained. This is not surprising for ${\bar{n}}_{{\rm{p}}}$, because this quantity only depends on the distribution of the transit parameter (Youdin 2011; Tremaine & Dong 2012). In fact, as we prove in the Appendix,

$$\begin{eqnarray}&&{\bar{n}}_{{\rm{p}}}=\displaystyle \frac{{\displaystyle \sum }_{j=1}^{K}\,{{jN}}_{j}}{{ \mathcal N }\langle \epsilon \rangle }.\end{eqnarray} \tag{ 22 }$$

That is, the average number of planets per star is given by the total number of tranets, the total number of stars, and the average probability that one planet transits. For our sample, the above equation gives ${\bar{n}}_{{\rm{p}}}=0.90\pm 0.03$, which agrees with our model outputs. Our constraint on ${\bar{n}}_{{\rm{p}}}$ also agrees with previous studies (e.g., Fressin et al. 2013; Petigura et al. 2013).

As the left panel of Figure 9 indicates, the total fraction of Sun-like stars with Kepler-like planets is also well constrained. To avoid the confusion with the general fraction F_p (for arbitrary planet sizes and orbital distances), we introduce η_Kepler for this specific fraction.⁹ The posterior distribution gives

$$\begin{eqnarray}&&{\eta }_{\mathrm{Kepler}}=30\pm 3 \% .\end{eqnarray} \tag{ 23 }$$

This is a factor of ∼2 lower than previous estimates (Fressin et al. 2013; Petigura et al. 2013; Winn & Fabrycky 2015). With the determinations of ${\bar{n}}_{{\rm{p}}}$ and η_Kepler, we also find that on average each Kepler-like planetary system has 3.0 ± 0.3 planets, as shown in the right panel of Figure 9. Our work is the first to determine this average multiplicity.

Zoom In Zoom Out Reset image size

How could the total fraction be constrained so well even though the individual components f_k were not? We will provide the explanation in Section 5.1, but the conclusion is that this results from some property of the matrix G that relates the nature of the transit probabilities. Therefore, our result that only 30% of Sun-like stars host Kepler-like planets is robust and, in particular, it is not sensitive to the details of TTV multiplicity function or our TTV modelings.

Here, we answer the question from Section 4.3: why could η_Kepler be constrained so well?

Following Section 3.1, the total number of transiting systems reads

$$\begin{eqnarray}&&\displaystyle \sum _{j=1}^{K}{\bar{N}}_{j}={ \mathcal N }\displaystyle \sum _{k=1}^{K}{f}_{k}\displaystyle \sum _{j=1}^{k}{g}_{{jk}}.\end{eqnarray} \tag{ 24 }$$

Unfortunately, the quantity ${\sum }_{j=1}^{k}{g}_{{jk}}$ (i.e., the probability to see at least one tranet) is not a constant. Otherwise the determination of F_p would be straightforward. However, unless the period distribution of planets in multiple systems is dramatically different from the period distribution of planets in single systems, the probability of having at least one tranet out of k (k ≥ 1) planets is no less than the probability of seeing one tranet in the single-planet systems. Mathematically, this implies ${\sum }_{j=1}^{k}{g}_{{jk}}\geqslant {g}_{11}$. Therefore, the above equation gives an upper limit on F_p,

$$\begin{eqnarray}&&{F}_{{\rm{p}}}\equiv \displaystyle \sum _{k=1}^{K}{f}_{k}\leqslant \displaystyle \frac{1}{{ \mathcal N }{g}_{11}}\displaystyle \sum _{k=1}^{K}{\bar{N}}_{k}.\end{eqnarray} \tag{ 25 }$$

For our sample, this gives η_Kepler ≤ 62%.

A second relevant quantity is the number of transit singles,

$$\begin{eqnarray}&&{\bar{N}}_{1}={ \mathcal N }\displaystyle \sum _{k=1}^{K}{f}_{k}{g}_{1k}.\end{eqnarray} \tag{ 26 }$$

Again, the quantity g_1k, the probability to see one tranet in a k-planet system, is not conserved: the more planets a system has and the hotter (i.e., larger inclination dispersion) the system is, the larger this quantity will be.

However, we find that the combination of the two quantities, ${g}_{1k}+{\sum }_{j=1}^{k}{g}_{{jk}}$, remains roughly constant for broad ranges of ${\sigma }_{i,5}$, α, and k. This can be seen in Figure 10, which illustrates the dependence of this quantity on various model parameters. Here, because only the G matrix is involved, we use the deterministic approach of Tremaine & Dong (2012) to compute G. This approach requires truncating the Legendre series at a certain threshold l_max. Unlike Tremaine & Dong (2012), who used l_max = 50, we choose a much higher threshold, l_max = 3000, which is necessary in order to have the elements of G for flat and multi-planet systems (i.e., upper-left corner of the left panel of Figure 10) converge.

Zoom In Zoom Out Reset image size

The conservation of the quantity ${g}_{1k}+{\sum }_{j=1}^{k}{g}_{{jk}}$ means that the more planets a system has and the hotter the system is, the larger this joint probability is, and the two separate probabilities, g_1k and ${\sum }_{j=1}^{k}{g}_{{jk}}$, compensate each other if the number of planets increases while the inclination dispersion decreases. Figure 10 also suggests that the weighted mean

$$\begin{eqnarray}&&\left\langle {g}_{1k}+\displaystyle \sum _{j=1}^{k}{g}_{{jk}}\right\rangle \equiv \displaystyle \frac{1}{{F}_{{\rm{p}}}}\displaystyle \sum _{k=1}^{K}{f}_{k}\left({g}_{1k}+\displaystyle \sum _{j=1}^{k}{g}_{{jk}}\right)\approx 0.11,\end{eqnarray} \tag{ 27 }$$

so that the total fraction of stars with planets can be given by

$$\begin{eqnarray}&&{F}_{{\rm{p}}}=\displaystyle \frac{1}{{ \mathcal N }\langle {g}_{1k}+{\displaystyle \sum }_{j=1}^{k}{g}_{{jk}}\rangle }\left({N}_{1}+\displaystyle \sum _{j=1}^{K}{N}_{j}\right).\end{eqnarray} \tag{ 28 }$$

With numbers from our sample, this estimator directly gives η_Kepler = 30%, as long as the following three conditions are met.

First, the k = 1 term does not provide a significant contribution to the weighted mean. That is, there is no large population of intrinsic singles. This has been discussed qualitatively in Section 1 and here we provide a simple quantitative argument. Of all the tranets in our sample, 432 and 395 are in transit singles and transit multis, respectively. Of the subset of those that show detected TTV signals, there are 23 in transit singles and 26 (excluding the double counting from tranet pairs both showing TTVs) in transit multis. For the transiting planets in transit multis, the probability of showing TTV is 26/395 = 6.6%. For those in transit singles, the same probability is 23/432 = 5.3% < 6.6%, which we interpret as the blending of intrinsic singles in the transit singles. Therefore, one finds that at most 19% (or 82 systems) of the transit singles are intrinsic singles. With the mean transit probability ($\langle \epsilon \rangle =0.03$) and the total number of surveyed stars (30759), the fraction of intrinsic singles (f₁) is at most 9%. Although this upper limit is derived based on a subset of the Kepler catalog and a specific TTV table (Holczer et al. 2016), the conclusion remains the same even if one uses the multiplicity fraction of the overall Kepler catalog or a different TTV table (Ofir et al. 2018). Even with this upper limit for f₁, the weighted mean of ${g}_{1k}+{\sum }_{j=1}^{K}{g}_{{jk}}$ only varies by 10%.

Second, the more planets there are in the system, the smaller the planet inclination dispersion is. Together with previous investigations of the multi-tranet systems (Lissauer et al. 2012; Tremaine & Dong 2012; Fabrycky et al. 2014), the result that there is no large population of intrinsic singles suggests that a significant fraction of planetary systems must have large planet–planet mutual inclinations. Although it is in principle possible that these are high multiples, it is much more likely that such systems with high mutual inclinations are low multiples, and indirect evidence from the planet eccentricity study also supports this (Xie et al. 2016).

Finally, there is no large population of very high (k ≥ 7) multiples. Given the previous result, these k ≥ 7 multiples should be very flat. With our inner (∼1 day) and outer (400 days) period boundaries, the weighted quantity ${f}_{1k}\,+{\sum }_{j=1}^{K}{f}_{{jk}}\lesssim 0.25$, and the probability to see seven tranets is at least R_⊙/(1.06 au) = 0.0044. Therefore, the fact that we do not have any 7-tranet system in a sample of 30,759 stars sets an upper limit on the fraction of systems with k ≥ 7, f_≥7 < 2.2% (95% confidence level). Including these very high multiples in the calculation of the weighted mean changes it by 10% in the opposite direction as the intrinsic singles does.

The above arguments explain why the total fraction can be constrained very well even though the individual components cannot, and more importantly, confirm that our determination of the total fraction of Sun-like stars with Kepler-like planets, η_Kepler, is fairly robust, and does not depend on the details of our modeling or the TTV multiplicity function. Furthermore, because the average number of planets per star, ${\bar{n}}_{{\rm{p}}}$, is determined independently from the inclination distribution (Equation (22)), the average multiplicity is also robustly measured.

Equation (28) also points toward a robust and straightforward way to determine the total fraction of stars with planets, which has practical applications. As Zhu et al. (2016) noted, the fraction of stars with planets is better than the average number of planets per star for the purpose of quantifying the correlation between planet formation efficiencies and stellar properties (such as metallicity). It is nevertheless the latter that has been broadly used.

In Section 1 we explained qualitatively why previous studies (Fressin et al. 2013; Petigura et al. 2013) overestimated η_Kepler. Now we explain it in a more quantitative way. We focus here on the transit studies, but the conclusion should apply to RV studies (e.g., Mayor et al. 2011) as well given that the same statistical approach was used.

There are two primary differences in our study and previous studies: the parameter space under investigation and the statistical method. Here, we discuss the impact of the former. We use Fressin et al. (2013) as the example of previous studies. Fressin et al. (2013) took into account both the geometric transit probability and the pipeline detection efficiency, and concluded that 52% of Sun-like stars should have at least one planet with R_p > 0.8 R_⊕ and P < 85 days. Our result that only 30% of Sun-like stars host Kepler planets comes from studying the Kepler planets as a whole and only accounting for the geometric transit probability. However, the parameter space we study is inclusive of the parameter space investigated in Fressin et al. (2013). Specifically, the pipeline detection efficiency would only increase the number of systems with planets of R_p > 0.8 R_⊕ and P < 85 days by 10%. Such a small change is far from what is needed to explain the discrepancy between our work and Fressin et al. (2013).

All previous studies on the fraction of stars with planets used the probability to detect the most detectable planet for the probability to detect at least one planet. In those transit studies (Fressin et al. 2013; Petigura et al. 2013), by reducing the number of tranets in all systems to one, they assumed that the resulting average number of planets per star should be the total fraction of stars with planets. Using the distribution of the transit parameter determined by transiting planets in transit singles (the red curve in Figure 4),¹⁰ we find that the probability that a typical planet transits is g₁₁ = 0.025. According to Equation (22), the resulting average number of planets per star is 0.77. Taking this value for η_Kepler would mean an overestimation by a factor of 2.6. Note that this value (77%) also exceeds the upper limit (62%) we derived in Section 5.1 under the very general condition.

The bulk of the Kepler planets are the so-called super Earths (planets with radii between Earth and Neptune) and the majority of such planets discovered by Kepler reside within ∼100 days. Such super Earths are absent in our own solar system. We find that the average number of such planets a system hosts is ∼3, and that the dispersion of planetary inclinations is a steep function of the intrinsic multiplicity. As Figure 11 shows, a system with k ≥ 5 planets within 400 days is pretty flat, with an inclination dispersion within 1°. But as the number of planets reduces, the system puffs up and the planets can be mutually inclined up to ∼10^◦. An isotropic distribution for planets in two-planet systems is not completely ruled out by the data. Together with the statistical result of eccentricities from transit durations (Xie et al. 2016), we now have a consistent picture that systems with fewer planets are dynamically hotter. In the following, we briefly cast these results in light of formation theories that have been proposed.

Zoom In Zoom Out Reset image size

The formation of super Earths remain unresolved. In the standard core accretion theories (Ida & Lin 2004; Mordasini et al. 2009), in which the cores of these planets are formed at large distances when the gas disk is still fully present, these planets are migrated rapidly to near the inner edge of the proto-planetary disks. It is thought that MMRs naturally get set up between planets that are migrating with different speeds, preventing them from being engulfed by the host stars. However, such a story predicts an abundance of MMRs, in contrast with the observed period ratios in which MMRs feature minimally (Lissauer et al. 2011; Fabrycky et al. 2014). While multiple scenarios have been proposed to break the planets out of MMRs (Delisle et al. 2014; Goldreich & Schlichting 2014; Chatterjee & Ford 2015; Liu et al. 2017), it is unclear that the overall period distribution can be explained. It is also unclear in such a framework how to account for the diverse multiple systems and their current dynamical states (but see Izidoro et al. 2017).

An alternative scenario, first proposed by Hansen & Murray (2013), favors in situ formation, in which these planets are locally assembled.¹¹ This is similar to the conventional story for the formation of terrestrial planets (e.g., Kokubo & Ida 1998, 2000; Raymond et al. 2004). Using N-body simulations, Hansen & Murray (2013) studied planet assembly, starting from a large number of proto-planets with a total mass of 20 M_⊕ within 400 days. They further assumed that all impacts are fully accretional and that there is no external source of dissipation (by gas or by planetesimals). They presented a number of statistical properties for the resulting systems, which we proceed to compare against here.

First, Hansen & Murray (2013) found that the median number of planets in a system is 4, an average eccentricity dispersion of σ_e = 0.11, and an average inclination dispersion σ_i ∼ a few degrees. These are somewhat similar to our findings of 3 planets per system, and $\langle {{kf}}_{k}{\sigma }_{i,k}\rangle \sim 4^\circ$ ($\langle {f}_{k}{\sigma }_{i,k}\rangle \sim 6^\circ$) with 2 ≤ k ≤ 7. Second, Hansen & Murray (2013) noted that their simulations predicted too few single-tranet systems. This remains a problem even after this work. The mismatch seems to be mainly due to the differences in their predicted multiplicities: f₁ = f₂ = 0 and f₃ ≪ f₄, f₅ in Hansen & Murray (2013), compared to f_k ∼ 5%–10% for k = 1, 2, 3, 4 and much smaller for k > 4 in our work (Figure 8).

We show that lower k systems are dynamically hotter (also see Xie et al. 2014, for the eccentricity aspect). What is the reason behind this?

One possibility is that this is a natural outcome of the stability requirement. One naively expects that a more packed system (higher k) has to have a lower dispersion to avoid dynamical instability. From Pu & Wu (2015), the critical number of mutual Hill radii that a system can have is ${K}_{\mathrm{crit}}=2+k+27/5\times {\sigma }_{i,k}{[3/(2q)]}^{1/3}$, with q = M_p/M_⋆. Here, we have assumed that σ_e ∼ 2σ_i and that the dependence of K_crit on σ_i,k obtained for k = 5 can be applied to smaller k values.¹² Further assuming that the systems extend one decade (∼0.1–1 au) in semimajor axis, we then find that the total number of Hill spheres is ln (10)[3/(2q)]^1/3, and that the average separation (in unit of Hill radii) between planets in a k-planet system is ln (10)[3/(2q)]^1/3/k. By equating this average separation and the critical number K_crit, we get that the critical inclination dispersion is σ_i,k ∼ 24°/k − 08(1 + k/2) (with q = 10⁻⁵). This boundary follows a decreasing trend with k, which is similar to but less steep than our result in Figure 11. However, given the approximate nature of the arguments above, we are unable to quantitatively compare this prediction with our result, specifically at k = 2, where the extrapolation from larger k likely breaks down. Further theoretical work will be required to address whether stability is responsible for the observed inclination dispersion trend.

It is also possible that the correlation between σ_i,k and k reflects the formation process via giant impacts. Unfortunately, Hansen & Murray (2013) did not explicitly present their eccentricity/inclination dispersions as functions of intrinsic multiplicities, which would have made a pivotal comparison against our Figure 11. An analytic version of this formation process was put forward by Tremaine (2015), which used the ergodic approximation and dynamical stability. This model predicts that the dispersion in eccentricities ∝1/k, which is similar to the stability argument given above if we assume σ_i ∝ σ_e. In practice, this model does not provide any prediction for the inclinations as it has assumed coplanar orbits.

Some theoretical works suggest that distant giant planets could drive up the mutual inclinations of the inner planetary systems and/or decrease their multiplicity by dynamical instabilities (Johansen et al. 2012; Huang et al. 2017; Lai & Pu 2017; Mustill et al. 2017; Pu & Lai 2018). If this is the primary channel for the observed features, one would expect that the stellar hosts of single tranets (low j and/or large σ_i) should reflect the well-known giant planet-metallicity correlation (e.g., Fischer & Valenti 2005). However, as the lower middle panel of Figure 3 shows, the metallicity distribution of the stars with single tranets is similar to that of stars with multiple tranets. A similar result is found by looking at the sample of single-tranet systems with TTVs, which should be indicative of multi-planet systems with large inclination dispersions. Therefore, these observations indicate that the perturbation by the unseen distant giant planet is likely not the primary cause of the observed inclination dispersions.

Alternatively, host stars with a large quadrupole moment can also excite planet–planet mutual inclinations (Spalding & Batygin 2016; Spalding et al. 2018). Unfortunately, our current sample does not have enough hot (T_eff > 6200 K) stars, for which this mechanism is expected to be most efficient, to test this hypothesis.

We have determined the intrinsic architecture of Kepler planetary systems and now we discuss briefly how our solar system fits in this revised picture.

First of all, with η_Kepler = 30%, our solar system belongs to the majority of Sun-like stars that do not have any planets detectable by Kepler. As was first pointed out by Chiang & Laughlin (2013), the Kepler systems contain more solid mass in their inner regions than is expected from a minimum-mass solar nebula. The typical mass for Kepler planets is ∼3 M_⊕ (Marcy et al. 2014; Fulton et al. 2017; Hadden & Lithwick 2017; Owen & Wu 2017). So these systems contain from 3 to 18 M_⊕ of solid masses, in contrast to our own 2 M_⊕ (within 400 days). As far as the inner (≲1 au) planetary system is concerned, it seems plausible that the difference arises because our system is a slightly low-metallicity version of the typical Kepler system.¹³

Even though none of the solar system planets are detectable by Kepler, it is interesting to note that our own system shares at least two similarities with the Kepler systems, as has been shown in Figure 11. First, similar to the typical Kepler system, we also have three planets within 400 days. Second, if only considering these three planets (Mercury, Venus, and Earth), the average inclination relative to the invariable plane, 35, is consistent with the inclination dispersions of 3-planet Kepler systems. It is important to note, however, that the higher mass Kepler systems are typically thought to form in a substantially shorter amount of time (prior to the dispersal of the natal disk) as compared to the solar system terrestrial planets (Morbidelli et al. 2012; Chiang & Laughlin 2013). It is therefore unclear whether the similarities between the solar system and Kepler systems are purely coincidental or are representative of a more fundamental behavior of general planetary systems.