About 30% of Sun-like Stars Have Kepler-like Planetary Systems: A Study of their Intrinsic Architecture
Wei Zhu Cristobal Petrovich Yanqin Wu Subo Dong Jiwei Xie

UTF8gbsn




Canadian Institute for Theoretical Astrophysics, University of Toronto, 60 St. George Street, Toronto, ON M5S 3H8, Canada
Wei Zhu
weizhu@cita.utoronto.ca


Canadian Institute for Theoretical Astrophysics, University of Toronto, 60 St. George Street, Toronto, ON M5S 3H8, Canada
Centre for Planetary Sciences, Department of Physical & Environmental Sciences, University of Toronto at Scarborough, Toronto, Ontario M1C 1A4, Canada


Department of Astronomy and Astrophysics, University of Toronto, 50 St. George Street, Toronto, ON M5S 3H4, Canada


Kavli Institute for Astronomy and Astrophysics, Peking University, Yi He Yuan Road 5, Hai Dian District, Beijing 100871, China


School of Astronomy and Space Science & Key Laboratory of Modern Astronomy and Astrophysics in Ministry of Education, Nanjing University, Nanjing 210093, China


    We constrain the intrinsic architecture of Kepler planetary systems by modeling the observed multiplicities of the transiting planets (tranets) and their transit timing variations (TTVs). We robustly determine that the fraction of Sun-like stars with Kepler-like planets, η_ Kepler, is 30±3%. Here Kepler-like planets are planets that have radii R_ p≳ R_⊕ and orbital periods P<400 days. Our result thus significantly revises previous claims that more than 50% of Sun-like stars have such planets.
    Combining with the average number of Kepler planets per star (∼0.9), we obtain that on average each planetary system has 3.0±0.3 planets within 400 days.
    We also find that the dispersion in orbital inclinations of planets within a given planetary system, σ_i,k, is a steep function of its number of planets, k. This can be parameterized as σ_i,k∝ k^α and we find that -4<α<-2 at 2-σ level.
    Such a distribution well describes the observed multiplicities of both transits and TTVs with no excess of single-tranet systems. Therefore we do not find evidence supporting the so-called “Kepler dichotomy.”
    Together with a previous study on orbital eccentricities, we now have a consistent picture: the fewer planets in a system, the hotter it is dynamically.
    We discuss briefly possible scenarios that lead to such a trend. Despite our Solar system not belonging to the Kepler club, it is interesting to notice that the Solar system also has three planets within 400 days and that the inclination dispersion is similar to Kepler systems of the same multiplicity.
                                     






§ INTRODUCTION
 

    The term “planet occurrence rate” has two different interpretations: the average number of planets per star, and the fraction of stars with planets. These two quantities are different, unless all planetary systems have only one planet. With transiting planets (tranets, ) from surveys such as the Kepler mission <cit.>, one can constrain the first but not the second <cit.>, unless assumptions of the intrinsic architecture (e.g., the orbital inclination distribution and/or the intrinsic multiplicity function) are made. This is because, to determine the average number of planets per star, one needs to compute the probability that individual planet transits, which only involves the orbital period of the planet (or more precisely, the ratio of stellar radius to orbital separation, R_⋆/a). The distribution of orbital periods can be reconstructed from the orbital periods of observed tranets. However, to determine the fraction of stars with planets, one needs to compute the probability that a given star has at least one tranet. This involves the orbital inclinations of all planets around this star, in addition to the orbital periods of these planets. The distribution of planetary inclinations cannot be reconstructed from the tranet sample, because by definition all tranets have ∼90^∘ inclinations.

The commonly accepted result says that more than 50% of Sun-like stars have Kepler-like planets <cit.>. However, they used the transit probability of the innermost planet (as in , or the most easily detected one as in ) as the probability that at least one planet transits. Therefore, their estimates are only valid under the assumption that all detected planets are in multi-planetary systems on coplanar orbits.

Either using the combined constraints of Kepler and Radial Velocity (RV) data <cit.> or transit duration distributions normalized by orbital velocities <cit.>, multi-tranet systems are found to be on average nearly co-planar. However, these methods cannot be applied to the Kepler single-tranet systems, which contribute over half of the detected tranets.

<cit.> modeled the observed multiplicity function of Kepler and found the single-tranet systems are in excess to a single simulated underlying planet population. This is sometimes called the “Kepler dichotomy” <cit.>. However, <cit.> showed that modeling the observed multiplicity function from Kepler data alone cannot arrive at a reliable conclusion on the intrinsic multiplicity function due to its degeneracy with inclination distribution.

The transit timing variation (TTV) technique can help break the degeneracy between intrinsic multiplicity function and inclination distribution <cit.>. Although TTV is the behavior of the transiting planet, it can reveal the existence of the non-transiting companion <cit.>. If there is indeed a large population of intrinsic singles, then transiting planets in the single-tranet systems should have considerably smaller probability to show TTV signals than the transiting planets in multi-tranet systems. However, this is not supported by the large and uniform TTV catalogs. For example, <cit.> found that of the total 260 Kepler planets that showed TTV signals, 121 were in single-tranet systems. The larger TTV catalog of <cit.> that is more sensitive to smaller TTV amplitudes gives a similar result. Both strongly indicate that transiting planets in transit singles and transit multiples have similar probability to show TTV signals, and therefore that there is no large population of intrinsic singles.

Another evidence against the assumption that all Kepler planets are coplanar comes from the study of the distribution of planet eccentricities. Using the distribution of transit durations, <cit.> found that single tranets have on average substantially larger eccentricities than multiple tranets. Because the dispersions of orbital eccentricities and inclinations are generally expected to be correlated <cit.>, this result suggests that systems with fewer number of planets may have larger mutual inclinations.

If a significant fraction of planets are in multi-planet systems with larger mutual inclinations than previously thought, then the probability that one star is seen to have at least one tranet increases. Therefore, the total fraction of Sun-like stars with at least one Kepler-like planet will decrease.

For similar reasons, the statistical studies based on the radial velocity (RV) samples also overestimated the fraction of stars with planets, by using the probability to detect the most detectable planet for the probability to detect at least one planet <cit.>. Although this overestimation is less severe for giant planets because of their low multiplicity rate, it can significantly reduce the fraction of stars with lower-mass planets (such as super-Earths).

In this study, we combine the information of transiting planets and their non-transiting companions as inferred by TTVs, to constrain the intrinsic architecture of planetary systems. We focus on Sun-like stars in this paper. In Section <ref> we construct the transit and TTV multiplicity functions based on a homogeneous sample. We then forward model these functions to constrain the intrinsic architecture in Section <ref>. Our results are presented in Section <ref>, and discussed in more details in Section <ref>.




§ KEPLER-LAMOST SAMPLE
 

To select Sun-like (FGK-type dwarf) stars for our study, we rely on the spectroscopic data from the Large Sky Area Multi-Object Fiber Spectroscopic Telescope (LAMOST, also known as Goushoujing Telescope, ), which had surveyed over 30% of all Kepler targets by 2017 (DR4), with no bias toward planet hosts <cit.>. The derived stellar parameters are accurate at least for main-sequence stars, as previous studies <cit.> have shown. The sample selection is similar to <cit.> and <cit.>. In short, we find 30,759 stars with effective temperature T_ eff in the range 4700-6500 K and stellar surface gravity logg>4.0 (in cgs unit), based on the stellar parameters derived by the LAMOST official pipeline <cit.>. We then cross-match this stellar catalog with the planet candidate table from Kepler data release 23 <cit.>, and find 1635 KOIs. Then we remove KOIs that meet the following criteria:

    
  * Identified by <cit.>, <cit.>, and <cit.> as false positives; 484 are removed.
    
  * KOIs with transit S/N<7.1 according to <cit.>; 108 are removed.
    
  * KOIs with P>400days; 34 are removed.
    
  * KOIs with R_p>20 R_⊕; here R_p is computed using LAMOST stellar parameters. 108 are removed.
    
  * KOIs that are in single-tranet systems and have large False Positive Probabilities (FPP >68%, ); 74 are removed.

The last criterion is not applied to the multi-tranet systems, which overall have very low false positive rates <cit.>. As <cit.> pointed out, their FPPs for multi-tranet systems could have been inflated by the effect of unidentified TTVs. Indeed, Kepler-23b <cit.> and Kepler-50b <cit.> are both confirmed planets, but have FPP>0.68 according to <cit.>. This criterion is applied to the single-tranet systems, because of their overall high false positive rates and low TTV fractions. Of the 74 single-tranets removed, only two have TTV signals according to <cit.>, the inclusion or exclusion of which does not affect our results.

In the end, we have 827 planets (or planet candidates) around 589 stars. The transit multiplicity function, denoting the number of systems as a function of number of tranets in each system, is (N_1, N_2, N_3, N_4 ,N_5, N_6)=(432, 99, 42, 11, 3, 2), and no system with more than six tranets. We show in Figure <ref> the radii and orbital periods for tranets in our sample. The observed transit multiplicity function is illustrated in Figure <ref>.







lcccccc
The two-sample KS test p values for different combinations of subsamples and different (planetary and stellar) parameters.


Input subsamples     R_ p     P     R_⋆     M_⋆     [Fe/H]     T_ eff
Transit singles & Transit multis     0.055     0.003^a     0.80     0.52     0.66     0.55 

Transit singles & TTV singles             0.09     0.15     0.96     0.15 

Transit multis & TTV multis             0.10     0.08     0.99     0.63 

TTV singles & TTV multis     0.08     0.09     0.014^b     0.005^b     0.95     0.074 

^a This small p value is likely due to the geometric effect. See Appendix <ref> for more discussions. 

^b Given that TTV singles are statistically similar to transit singles and that TTV multis are statistically similar to transit multis, these two small p values are likely due to a random sampling effect.


To find out how many of these tranets show TTV signals, we cross-match with the TTV catalog of <cit.>, which was produced based on a uniform search among over 2600 KOIs with relatively high S/N. We find that there are (23, 12, 7, 4, 0, 1) systems in our sample that have (1, 2, 3, 4, 5, 6) tranets and at least one of the tranet shows TTV signals. We dub this the TTV multiplicity function, and also show it in Figure <ref>.

In our analysis that follows, we will make use of all the transiting planets in our sample, which includes the subsample exhibiting TTV signals. In what follows, we argue that the subsample exhibiting TTVs and those that do not are drawn from a common population, validating our usage of both subsets within a common statistical framework. Based on whether one system has single or multiple tranets and whether any of the tranets show TTV signals, we can divide the whole sample into four subsamples: transit singles, transit multis, TTV singles, and TTV multis. In the transit singles (transit multis) subsample, we do not exclude planets/stars in the TTV singles (TTV multis) subsample. Because the transiting planets outnumber the TTV planets by an order of magnitude, the inclusion or exclusion of the ones with TTVs does not make any noticeable difference. Figure <ref> shows the cumulative distributions of planetary (radius and orbital period) and stellar (radius, mass, [Fe/H], and T_ eff) properties of planets/stars in these subsamples, and Table <ref> provides the two-sample Kolmogorov-Smirnov (KS) test p values between selected subsamples. There are a number of notable features. First, transit singles and transit multis are statistically similar in almost every index, in particular of the stellar parameters, suggesting that most of the transit singles are likely drawn from the same underlying population and gone through similar formation processes as the transit multis. This is an indication that a substantial fraction of the transit singles are in fact intrinsic multiples, a conclusion that we come to endorse later in the paper.

Another notable feature is that the TTV planets prefer to have slightly larger planetary radii, and are at slightly longer periods, than the rest. This is expected. To enable TTV detections, the planet transits tend to be deeper and have longer periods (; see also Equation (<ref>)). However, the stellar properties of systems with and without TTVs (i.e., TTV singles vs. transit singles, and TTV multis vs. transit multis) are statistically similar, since their two-sample KS test p values given in Table <ref> are all above the standard threshold (0.05). Therefore, there is no reason to suspect that planets in the TTV sample are drawn from a different population than the tranets are. In this work, we explicitly assume that the TTV planets, despite their relative proximity to mean-motion resonances (MMRs), are not special and that their abundances can be used to constrain the overall planet population. Besides the similarity in stellar parameters, another supporting evidence is that, as new techniques are invented to detect lower-amplitude TTVs, more systems appear to show TTV signals <cit.>. The transition from showing and not showing TTVs is smooth rather than sharp.

Figure <ref> shows that, although both the transit and TTV multiplicity functions are monotonically decreasing (subject to statistical noises) as the transit multiplicity increases, the TTV multiplicity function has a weaker dependence on the transit multiplicity. This is because TTV is relatively insensitive to the inclination variations. Although such a feature prevents from using TTV as a characterization technique to precisely constrain the mutual inclination values (e.g., ), it indeed helps to use TTV as a detection technique to probe planet with a broader range of inclination values than the transit technique. The different slopes of the two multiplicity functions are the key to uncover the inclination distribution in multi-planet systems.





§ FORWARD MODELING THE OBSERVED MULTIPLICITY FUNCTIONS
 

The transit and TTV multiplicity functions, as illustrated in Figure <ref>, are both monotonically decreasing (subject to Poisson noises) as the intrinsic multiplicity k increases, but they behave quantitatively differently, with the TTV multiplicity function less dependent on k. Below we show that the transit and TTV multiplicity functions can be simultaneously well described when the inclination dispersion of the k-planet system is a power-law function of the intrinsic multiplicity k. Using this relation, we can constrain the intrinsic multiplicity vector 𝐅≡ (f_1, f_2, ⋯, f_k), where f_k is the fraction of Sun-like stars with k Kepler-like planets. This section describes the model we use to fit the observed transit and TTV multiplicity functions.




 §.§ Notations
 

To facilitate further discussions, we introduce a few mathematical notations here. Following <cit.>, we use a matrix 𝐆 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

    g_jk = 0 ,   if j>k .

If the parameter space that is of interest can fit in at most K planets, then 𝐆 should be a K× K upper-triangular matrix.
[Note that the <cit.> extended their notations to j=k=0 and thus their 𝐆 matrix had (K+1)× (K+1) dimensions.]

We use symbol 𝒩 for the number of stars in the sample, and the vector 𝐍=(N_1, ⋯, N_K) for the observed transit multiplicity function. With the intrinsic multiplicity vector 𝐅, the expectation of the transit multiplicity function can be given as

    N̅_j = ∑_k=j^K g_jk𝒩f_k = ∑_k=1^K g_jk𝒩 f_k,   or 𝐍̅ = 𝐆· (𝒩𝐅) .

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, n̅_ p, are given by

    F_ p=∑_k=1^K f_k,  n̅_ p=∑_k=1^K kf_k ,

respectively. The ratio of these two quantities, n̅_ p/F_ p, gives the average number of planets per planetary system, which we call the average multiplicity.

We also introduce the matrix 𝐓 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_1,⋯, M_K), and the expectation of this is given by

    M̅_j = ∑_k=j^K t_jk𝒩f_k,  or 𝐌̅ = 𝐓· (𝒩𝐅) ;





 §.§ Model Ingredients
 



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.



  §.§.§ Distribution of Transit Parameters



Following <cit.>, we model the distribution of transit parameter ϵ as

    N/lnϵ∝(ϵ/ϵ_0)^a/1+(ϵ/ϵ_0)^b ,

which is essentially a broken power law but with smooth transition at ϵ_0. 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 <ref>). This reconstructed ϵ distribution using all tranets in our sample is shown as black dots in Figure <ref>, in which we also show the distribution from only transit singles for a reference.

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

    N/lnϵ = 0.36 (ϵ/ϵ_0)^0.04/1+(ϵ/ϵ_0)^3.18   (ϵ_0=0.074)

for 0.004<ϵ<0.6, and zero elsewhere. This yields a logarithmically flat distribution below ϵ_0, a result in agreement with previous studies <cit.>.



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

    cosI_ p = cosIcosi - sinIsinicosϕ .

The distribution of I is isotropic (∝sinI 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, sini) and k, and choose the normalization at k=5

    σ_i,k≡√(⟨sin^2 i ⟩) = (k/5)^α .

It is written in this form, so that the normalization factor, , can be determined separately from the distribution of transit duration ratios of planet pairs in five-planet systems (Section <ref>). 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 <cit.> distribution <cit.>

    P(i|κ_k) = κ_k sini/2sinhκ_k e^κ_kcosi .

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

    σ^2_i,k = ⟨sin^2i⟩ = 2/κ_k(κ_k-1/κ_k) .

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

The inclination dispersion σ_i,k, by its definition given by Equation (<ref>), has a maximum value of √(2/3), which can be achieved only when the distribution of i becomes isotropic. For any given , the upper bound on the inclination dispersion therefore sets a limit on α (and vice versa).



  §.§.§ 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 (<ref>), and then the orbital period is derived via

    P = (R_⊙/ au)^3/2ϵ^-3/2(ρ_⋆/ρ_⊙)^-1/2 .

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 <cit.> stability criterion, which requires that for any given planet pair,

    P_ out/P_ in > 1+2.2 q^2/7 = 1.16 .

The critical value is derived by assuming a characteristic planet-to-star mass ratio q=10^-4. The actual choice of this characteristic value has 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^-5) of Kepler planets, it is extremely rare to have a period ratio below 1.3 in actual observations <cit.>.



  §.§.§ TTV Detection Criteria
 

For any system with at least one tranet, we determine whether there is detectable TTV signal. Previous studies <cit.> have shown that in systems with detected TTVs, it is almost always the case that the TTV signal comes from first-order MMRs between two neighboring planets. We use this empirical result and only consider the closest one (if the tranet is innermost or outermost) or two planets in the TTV detection. We only consider the TTV signals from first-order (J:J-1=2:1, 3:2, 4:3, and 5:4) MMRs. Higher-order MMRs are too weak, and the additional first-order MMRs are either not allowed by or too close to the stability limit. We consider the TTV signal to be detectable, if all the following conditions are met for at least one of the chosen first-order MMRs:

    
  * The orbital period of the tranet is less than 200 days.
        [Although <cit.> included P<300days in their initial selections, they also required at least six transits observed in the Kepler window.]
    
  * The super period of the planet pair
        
    P_ sup≡P_ inP_ out/|JP_ in-(J-1)P_ out|

        is in the range 100-3000 days. Here P_ in and P_ out are the orbital periods of the inner and outer planets, respectively.
    
  * The TTV amplitude indicator
        
    P/Δ > 1.3× 10^3  days

        where P is the orbital period of the tranet, and Δ is the fractional separation to period commesurability <cit.>.
        
    Δ≡| P_ out/P_ inJ-1/J - 1 |  .


The first two criteria are used to mimic the conditions in actual TTV detections <cit.>. The threshold used in the third criterion, 1.3×10^3 days, is approximately the median value of P/Δ in identified TTV pairs of <cit.>. For the characteristic mass ratio (q=10^-4), this corresponds to a 30-min TTV amplitude, which is also the median amplitude of all TTV planets in <cit.>. Our results are insensitive to the numerical threshold adopted here, as we will show in Section <ref>. Specifically, reducing (or increasing) the threshold value on the right-hand side of Equation (<ref>) will allow more (or less) pairs to become TTV eligible. However, since this raises the fraction of TTV sample uniformly across all systems, this normalization change does not affect our results.




 §.§ Constraining σ_i,5 from Transit Duration Ratios
 



We acquire external constraints on the normalization factor  of the inclination dispersion relation (Equation <ref>). This is necessary because otherwise we would end up with a strong correlation between  and α.

Because we are constraining  separately, we are not limited to the planetary systems in our current sample. In fact, our sample only contains three five-tranet systems, and these provide 3× C_5^2=30 planet pairs. Instead, we find 15 five-tranet systems whose hosts are Sun-like stars from the California Kepler Survey <cit.>, which give us 150 planet pairs to constrain the inclination dispersion .

We adopt a similar approach as <cit.> to constrain the inclination dispersion. Each observed planet pair gives a quantity <cit.>

    ξ≡T_ dur,in/P_ in^1/3/T_ dur,out/P_ out^1/3 ,

where T_ dur is the transit duration (from first to fourth contact) and P is the orbital period. The subscripts “in” and “out” denotes 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 <cit.>. We can also write ξ in terms of the planet-to-star radius ratio r and the transit impact parameter b

    ξ = √((1+r_ in)^2-b_ in^2/(1+r_ out)^2-b_ out^2)

Because both T_ dur and P are much better measured than r or b in observations, the expression given by Equation (<ref>) 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; ). We then compute the ξ values for individual planet pairs, and show the distribution as the black histogram on the left panel of Figure <ref>. 

We then model this ξ distribution and attempt to constrain the inclination dispersion . 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-multiply and stable systems are not expected to have large eccentricities. Second, <cit.> has shown 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 , we produce the ξ distribution following the method of <cit.>. 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 (√(b_ out,max/b_ 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_ out,sim(P_ in/P_ out)^2/3 and dispersion /(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  in inclinations. If |b_ in,sim|≤ b_ 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 (<ref>). We repeat the whole process and generate 250 simulated pairs for each planet pair in the sample. The ξ distribution for the given  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 <cit.>, which used the p-value of Kolmogorov-Smirnov test. We compute the likelihood as

    ℒ = ∏_j=1^N∫ P_ sim(lnξ) exp[-(lnξ_j-lnξ)^2/2σ_lnξ,j^2]dlnξ ,

where P_ sim(lnξ) is the probability distribution of simulated ξ in logarithmic space, ξ_j is the observed value of ξ of the j-th 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 ℒ for 100 values of  equally spaced from 0.1^∘ to 4.0^∘, and show the results on the right panel of Figure <ref>. As the scatter plot shows, the likelihood reaches its maximum around =0.8^∘. The simulated ξ distribution for this best-fit  is also shown on the left panel of Figure <ref>. To find the different confidence levels of , we use a spline function to smooth the ln(ℒ/ℒ_ max) vs.  scatter plot, and find that the 1-σ 2-σ, and 3-σ confidence intervals, defined as ln(ℒ/ℒ_ max) ≥ -n^2/2 to be 0.65^∘-0.96^∘, 0.53^∘-1.16^∘, and 0.42^∘-1.44^∘ for n=1, 2, and 3, respectively.

For the subsequent modeling of the power-law index α and intrinsic multiplicity vector 𝐅, we will only consider values of  from the 3-σ confidence interval. The smoothed ln(ℒ/ℒ_ max) is also used as our prior on  unless specified otherwise.




 §.§ Monte Carlo Simulations


We use Monte Carlo simulations to compute the 𝐆 and 𝐓 matrices for a grid of  and α. We note that <cit.> provided an analytical formalism to compute the 𝐆 matrix. However, the time-limiting factor is always the computation of the 𝐓 matrix, which may not be done in the analytical way. We compute the 𝐓 matrix within a Monte Carlo, and this automatically produces the 𝐆 matrix.

For the intrinsic singles (k=1), there is no TTV signal (t_11=0), and only one parameter (g_11) needs to be calculated. This can be done analytically <cit.>: g_11=⟨ϵ⟩ = 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 <ref> 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 (<ref>)), whose parameter κ_k is determined by Equations (<ref>) and (<ref>). The inclination of the planet with respect to the line of sight is then given by Equation (<ref>). Finally, whether a planet is a tranet or not is determined by the transit criterion

    ϵ>cosI_ p .

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

With the above procedures, we are able to generate an 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 <ref>. We repeat the whole process and generate a large number of planetary systems, until each element in the 𝐓 matrix is determined to <2% precision.




 §.§ Modeling the Observed Multiplicity Functions


For a given intrinsic multiplicity vector 𝐅 and matrices 𝐆 and 𝐓, the probability to see the transit and TTV multiplicity functions as shown in Figure <ref> is given by

    ℒ = ∏_k=0^K N̅_k^N_kexp(-N̅_k)/N_k!×∏_k=1^K M̅_k^M_kexp(-M̅_k)/M_k! .

where N̅_k and M̅_k are given by Equations (<ref>) and (<ref>), respectively. Note that here we have extended the transit multiplicity function down to k=0, with f_0=1-∑_k≥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 (, α, and 𝐅)

    lnℒ = ∑_k=0^K (N_k lnN̅_k - N̅_k) + ∑_k=1^K (M_k lnM̅_k - M̅_k) .

Given the large number of dimensions, we choose the Markov Chain Monte Carlo (MCMC) algorithm that is implemented in  <cit.> as the optimization method. This way we also obtain the posterior distributions of F_k (k=1, ⋯, 6) for given sets of  and α.





§ RESULTS
 












 §.§ Planetary Inclination Dispersion


We find the maximum likelihood values for each grid point in the (, α) plane following the procedure detailed in the previous section. After imposing the prior probability on  from the transit duration ratios (Figure <ref>), we can determine the n-σ contours, where n (=1, 2, 3) refers to the number of σ and the contour is defined by ln(ℒ/ℒ_ 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 <ref>.

Figure <ref> indicates that the power-law index α is well constrained to be close to its lower bound -4, which is given by the normalization factor ≈0.8^∘ and σ_i,2≤√(2/3) (Section <ref>). Specifically, α<-3 at 1-σ level, α<-2 at 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 <ref> 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 <cit.> or a different TTV catalog <cit.> is used. Therefore, our constraints on the power-law index parameter is robust.

Figure <ref> 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.




 §.§ Intrinsic Multiplicity Vector


We stack the Markov chains from all MCMC runs, and discard entries that are more than 3-σ away (Δχ^2>9) from the best one. With this combined Markov chain, we can then investigate the constraints on individual components of the intrinsic multiplicity vector 𝐅.
[We confirm that this approach produces very similar but smoother posteriors than directly using the shape of the χ^2 (i.e., lnℒ) curve.]

The left panel of Figure <ref> 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_1, 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 <ref>.

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




 §.§ Overall Planet Occurrence Rates
 

Although the individual components of the intrinsic multiplicity vector 𝐅 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 n̅_ p, are found to be well constrained. This is not surprising for n̅_ p, because this quantity only depends on the distribution of the transit parameter ϵ <cit.>. In fact, as we prove in Appendix <ref>, 

    n̅_ p=∑_j=1^K jN_j/𝒩⟨ϵ⟩ .

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 n̅_ p=0.90±0.03, which agrees with our model outputs. Our constraint on n̅_ p also agrees with previous studies <cit.>.

As the left panel of Figure <ref> 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.
[This notation follows the well-accepted term η_⊕.]
The posterior distribution gives

    η_ Kepler = 30±3% .

This is a factor of ∼two lower than previous estimates <cit.>. With the determinations of n̅_ 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 <ref>. Our work is the first to determine this average multiplicity.

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 <ref>, but the conclusion is that, this results from some property of the matrix 𝐆 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.





§ DISCUSSION
 




 §.§ The Determination of η_ Kepler
 



Here we answer the question from Section <ref>: why η_ Kepler could be constrained so well?

Following Section <ref>, the total number of transiting systems reads

    ∑_j=1^K N̅_j = 𝒩∑_k=1^K f_k ∑_j=1^k g_jk .

Unfortunately, the quantity ∑_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 to have at least one tranet out of k (k≥ 1) planets is no less than the probability to see one tranet in the single-planet systems. Mathematically, this implies ∑_j=1^k g_jk≥ g_11. Therefore, the above equation gives an upper limit on F_ p

    F_ p≡∑_k=1^K f_k ≤1/𝒩 g_11∑_k=1^K N̅_k .

For our sample, this gives η_ Kepler≤ 62%.

A second relevant quantity is the number of transit singles

    N̅_1 = 𝒩∑_k=1^K f_k g_1k .

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+∑_j=1^k g_jk, remains roughly constant for broad ranges of , α, and k. 
This can be seen in Figure <ref>, which illustrates the dependence of this quantity on various model parameters. Here because only the 𝐆 matrix is involved, we use the deterministic approach of <cit.> to compute 𝐆. This approach requires to truncate the Legendre series at a certain threshold l_ max. Unlike <cit.> who used l_ max=50, we choose a much higher threshold, l_ max=3000, which is necessary in order to have the elements of 𝐆 for flat and multi-planet systems (i.e., upper left corner of the left panel of Figure <ref>) converge.

The conservation of the quantity g_1k+∑_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 ∑_j=1^k g_jk, compensate each other if the number of planets increases while the inclination dispersion decreases. Figure <ref> also suggests that the weighted mean

    ⟨ g_1k+∑_j=1^k g_jk⟩≡1/F_ p∑_k=1^K f_k ( g_1k+∑_j=1^k g_jk) ≈ 0.11 ,

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

    F_ p = 1/𝒩⟨ g_1k+∑_j=1^k g_jk⟩( N_1 + ∑_j=1^K N_j ) ,

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 significant contribution to the weighted mean. That is, there is no large population of intrinsic singles. This has been discussed qualitatively in Section <ref>, and here we provide a simple quantitative argument. Of all 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, 23 are in transit singles and 26 (excluding the double counting from tranet pairs both showing TTVs) are 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. One therefore finds that at most 19% (or 82 systems) of the transit singles are intrinsic singles. With the mean transit probability (⟨ϵ⟩=0.03) and the total number of surveyed stars (30759), the fraction of intrinsic singles (f_1) is at most 9%. Although this upper limit is derived based on a subset of the Kepler catalog and a specific TTV table <cit.>, the conclusion remains the same even if one uses the multiplicity fraction of the overall Kepler catalog or a different TTV table <cit.>. Even with this upper limit for f_1, the weighted mean of g_1k+∑_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 <cit.>, 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 <cit.>.

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+∑_j=1^K f_jk≲ 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, n̅_ p, is determined independently from the inclination distribution (Equation (<ref>)), the average multiplicity is also robustly measured.

Equation (<ref>) also points toward a robust and straightforward way to determine the total fraction of stars with planets, which has practical applications. As <cit.> have pointed out, 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.




 §.§ Only 30% of Sun-like Stars Host Kepler Planets


In Section <ref> we have explained qualitatively why previous studies <cit.> overestimated η_ Kepler. Now we explain it in a more quantitative way. We focus on the transit studies here, but the conclusion should apply to RV studies <cit.> 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 <cit.> as the example of previous studies. <cit.> 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 by 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 <cit.>. 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 <cit.>.

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 <cit.>, 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 <ref>),
[Because transit singles outnumber the transit multiples significantly. The ϵ distribution will be essentially the same regardless of whether the innermost tranets <cit.> or the most detectable tranets <cit.> of transit multiples are included.]
we find that the probability that a typical planet transits is g_11=0.025. According to Equation (<ref>) 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 <ref> under the very general condition.






 §.§ Inclination Dispersions and Multiplicities: Comparisons with Theories


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 <ref> shows, a system with k≥5 planets within 400 days is pretty flat, with 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 <cit.>, 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.

The formation of super-Earths remain unresolved. In  the standard core accretion theories <cit.>, 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 <cit.>. While multiple scenarios have been proposed to break the planets out of MMRs <cit.>, 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 ).

An alternative scenario, first proposed by <cit.>, favors in-situ formation, in which these planets are locally assembled. 
[Here, migration by the gas disk is artificially suppressed. One possibility for this is late assembly <cit.>.]
This is similar to the conventional story for the formation of terrestrial planets (e.g., ). Using N-body simulations, <cit.> 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, <cit.> 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 ⟨ kf_k σ_i,k⟩∼ 4^∘ (⟨ f_k σ_i,k⟩∼ 6^∘) with 2≤ k≤7. Second, <cit.> 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_1=f_2=0 and f_3≪ f_4,f_5 in <cit.>, 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 <cit.>. 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 <cit.> the critical number of mutual Hill radii that a system can have is K_ crit=2+k+27/5×σ_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. 
[The extrapolation might break down at k=2 as K_ crit is somewhat insensitive to changes in mutual inclinations for ≲40^∘ <cit.>.]
Further assuming that the systems extend one decade (∼ 0.1-1 AU) in semi-major axis, then we 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-0.8^∘(1+k/2) (with q=10^-5). This boundary follows a decreasing trend with k, which is similar to but less steep than our result in Figure <ref>. 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, <cit.> did not explicitly present their eccentricity/inclination dispersions as functions of intrinsic multiplicities, which would have made a pivotal comparison against our Figure <ref>. An analytic version of this formation process was put forward by <cit.>, 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 <cit.>. 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 <cit.>. However, as the lower middle panel of Figure <ref> 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 the planet-planet mutual inclinations <cit.>. 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.




 §.§ Solar System vs. Kepler Systems


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 planet detectable by Kepler. As was first pointed out by <cit.>, the Kepler systems contain more solid mass in their inner regions than is expected from a minimum-mass Solar nebula (MMSN). The typical mass for Kepler planets is ∼3 M_⊕ <cit.>. 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.
[The Solar system has giant planets beyond 1 au, which is atypical for Sun-like stars with Solar metallicity. If the outer planetary system is also taken into account, our Solar system no longer belongs to the majorities <cit.>.]

Even though none of the Solar system planets is detectable by Kepler, it is interesting to notice that our own system shares at least two similarities with the Kepler systems, as has been shown in Figure <ref>. 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, 3.5^∘, 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 <cit.>. 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.




We thank the anonymous referee for comments. We would also like to thank Andy Gould for comments on an earlier version of the manuscript.

This paper includes data collected by the Kepler mission. Funding for the Kepler mission is provided by the NASA Science Mission directorate.

This paper also uses data from the LAMOST survey. Guoshoujing Telescope (the Large Sky Area Multi-Object Fiber Spectroscopic Telescope, LAMOST) is a National Major Scientific Project built by the Chinese Academy of Sciences. Funding for the project has been provided by the National Development and Reform Commission. LAMOST is operated and managed by the National Astronomical Observatories, Chinese Academy of Sciences.
CP acknowledges support from the Jeffrey L. Bishop Fellowship and from the Gruber Foundation Fellowship. 
SD acknowledges Project 11573003 supported by National Natural Science Foundation of China (NSFC).



§ SIMPLE ESTIMATION OF AVERAGE NUMBER OF PLANETS PER STAR
 



As <cit.> pointed out, the average number of planets per star n̅_ p is directly given by the the total number of tranets ∑_j j N_j, the total number of stars 𝒩, and the average transit parameter ⟨ϵ⟩. See also our Equation (<ref>). With our mathematical notations, this is equivalent to

    ∑_j=1^k jg_jk = k⟨ϵ⟩ .

In other words, the number of tranets one expects to see in a k-planet system is proportional to k, regardless of the details of g_jk. This can be proved mathematically using the expression of 𝐆 in <cit.>. Below we provide another simple and robust proof, and discuss the associated assumption.

The transit probability is essentially the fractional area on a unit sphere where transit happens (ϵ>cosI_p). This can be shown graphically by projecting the 3D sphere onto a 2D plane, as done in Figure <ref> for 1-planet and 2-planet cases. See also <cit.> for the 3-planet case. For k=1, Equation (<ref>) reduces to the definition of ⟨ϵ⟩. For k=2, with the notations for different areas in Figure <ref> the left-hand-side of Equation (<ref>) is (A_1+A_2+A_3+A_4+A_5+A_6)+2(A_7+A_8)=(A_1+A_7+A_2+A_8+A_3)+(A_4+A_7+A_5+A_8+A_6)=g_11(p1) + g_11(p2), where p1 and p2 denote the planet 1 and 2, respectively. Therefore, the left-hand-side equals to the right-hand-side as long as the transit parameters (or approximately the separations) of the two planets are statistically no different. This is the only assumption that goes into Equation (<ref>). The k≥3 cases can be easily proved similarly.

Is the separation distribution of planets in singles different from the separation distribution of planets in multiples? If using tranets, one indeed sees a difference between these two distributions. For example, see the upper middle panel of Figure <ref> as well as Figure <ref>. However, this is not super surprising, as the detections of tranets around the same star are correlated due to the geometric effect. Even so, the difference in these two distributions is not very prominent. Perhaps a better sample to use is the sample of planets from RV observations. As <cit.> have discussed, the distributions of the semi-major axes of RV planets are statistically the same in single- and multiple-planet systems. See in particular the lower left panel of their Figure 2.




[Agol et al.(2005)]Agol:2005 Agol, E., Steffen, J., Sari, R., & Clarkson, W. 2005, , 359, 567 
[Ballard & Johnson(2016)]Ballard:2016 Ballard, S., & Johnson, J. A. 2016, , 816, 66 
[Borucki et al.(2010)]Borucki:2010 Borucki, W. J., Koch, D., Basri, G., et al. 2010, Science, 327, 977 
[Brakensiek & Ragozzine(2016)]Brakensiek:2016 Brakensiek, J., & Ragozzine, D. 2016, , 821, 47 
[Burke et al.(2015)]Burke:2015 Burke, C. J., Christiansen, J. L., Mullally, F., et al. 2015, , 809, 8 
[Chatterjee & Ford(2015)]CF2015 Chatterjee, S. & Ford, E. B. 2015, , 803, 33
[Chiang & Laughlin(2013)]Chiang:2013 Chiang, E., & Laughlin, G. 2013, , 431, 3444 
[Coughlin et al.(2016)]Coughlin:2016 Coughlin, J. L., Mullally, F., Thompson, S. E., et al. 2016, , 224, 12 
[Cui et al.(2012)]Cui:2012 Cui, X.-Q., Zhao, Y.-H., Chu, Y.-Q., et al. 2012, Research in Astronomy and Astrophysics, 12, 1197 
[Cumming et al.(2008)]Cumming:2008 Cumming, A., Butler, R. P., Marcy, G. W., et al. 2008, , 120, 531 
[De Cat et al.(2015)]Decat:2015 De Cat, P., Fu, J. N., Ren, A. B., et al. 2015, , 220, 19 
[Deck et al.(2013)]Deck:2013 Deck, K. M., Payne, M., & Holman, M. J. 2013, , 774, 129 
[Delisle et al.(2014)]Delisle2014 Delisle, J.-B., Laskar, J., & Correia, A. C. M. 2014, , 566, A137
[Dong et al.(2017)]Dong:2017 Dong, S., Xie, J.-W., Zhou, J.-L., Zheng, Z., & Luo, A. 2017, arXiv:1706.07807 
[Dong et al.(2014)]Dong:2014 Dong, S., Zheng, Z., Zhu, Z., et al. 2014, , 789, L3 
[Dong & Zhu(2013)]DongZhu:2013 Dong, S., & Zhu, Z. 2013, , 778, 53 
[Fabrycky & Winn(2009)]Fabrycky:2009 Fabrycky, D. C., & Winn, J. N. 2009, , 696, 1230 
[Fabrycky et al.(2014)]Fabrycky:2014 Fabrycky, D. C., Lissauer, J. J., Ragozzine, D., et al. 2014, , 790, 146 
[Fang & Margot(2012)]Fang:2012 Fang, J., & Margot, J.-L. 2012, , 761, 92 
[Figueira et al.(2012)]Figueira:2012 Figueira, P., Marmier, M., Boué, G., et al. 2012, , 541, A139 
[Fischer & Valenti(2005)]Fischer:2005 Fischer, D. A., & Valenti, J. 2005, , 622, 1102 
[Fisher(1953)]Fisher:1953 Fisher, R. 1953, Proceedings of the Royal Society of London Series A, 217, 295 
[Ford et al.(2012)]Ford:2012 Ford, E. B., Fabrycky, D. C., Steffen, J. H., et al. 2012, , 750, 113 
[Foreman-Mackey et al.(2013)]emcee:2013 Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, , 125, 306 
[Fressin et al.(2013)]Fressin:2013 Fressin, F., Torres, G., Charbonneau, D., et al. 2013, , 766, 81 
[Fulton et al.(2017)]Fulton:2017 Fulton, B. J., Petigura, E. A., Howard, A. W., et al. 2017, , 154, 109 
[Goldreich & Schlichting(2015)]GS2015 Goldreich & Schlichting, H. E. 2014, AJ, 147, 32
[Hadden & Lithwick(2017)]Hadden:2017 Hadden, S., & Lithwick, Y. 2017, , 154, 5 
[Hansen & Murray(2013)]Hansen2013 Hansen, B. M. S., & Murray, N. 2013, ApJ, 775, 53
[Hansen & Murray(2015)]Hansen:2015 Hansen, B. M. S., & Murray, N. 2015, , 448, 1044 
[Holczer et al.(2016)]Holczer:2016 Holczer, T., Mazeh, T., Nachmani, G., et al. 2016, , 225, 9 
[Holman & Murray(2005)]Holman:2005 Holman, M. J., & Murray, N. W. 2005, Science, 307, 1288 
[Huang et al.(2017)]Huang:2017 Huang, C. X., Petrovich, C., & Deibert, E. 2017, , 153, 210 
[Ida et al.(1993)]Ida:1993 Ida, S., Kokubo, E., & Makino, J. 1993, , 263, 875
[Ida & Lin(2004)]IdaLin:2004 Ida, S., & Lin, D. N. C. 2004, , 616, 567 
[Izidoro et al.(2017)]Izidoro:2017 Izidoro, A., Ogihara, M., Raymond, S. N., et al. 2017, , 470, 1750 
[Johansen et al.(2012)]Johansen:2012 Johansen, A., Davies, M. B., Church, R. P., & Holmelin, V. 2012, , 758, 39 
[Kokubo & Ida(1998)]kokubo1998 Kokubo, E., & Ida, S. 1998, Icar, 131, 171
[Kokubo & Ida(2000)]kokubo2000 Kokubo, E., & Ida, S. 2000, Icar, 143, 15
[Lai & Pu(2017)]Lai:2017 Lai, D., & Pu, B. 2017, , 153, 42 
[Lee et al.(2014)]lee2014 Lee E. J., Chiang E., & Ormel C. W., 2014, 797, 95

[Lissauer et al.(2011)]Lissauer:2011 Lissauer, J. J., Ragozzine, D., Fabrycky, D. C., et al. 2011, , 197, 8 
[Lissauer et al.(2012)]Lissauer:2012 Lissauer, J. J., Marcy, G. W., Rowe, J. F., et al. 2012, , 750, 112 
[Liu et al.(2017)]Liu2017 Liu, B., Ormel, C., & Lin, D. N. C. 2017, , 601, 15L
[Lithwick et al.(2012)]Lithwick:2012 Lithwick, Y., Xie, J., & Wu, Y. 2012, , 761, 122 
[Luo et al.(2015)]Luo:2015 Luo, A.-L., Zhao, Y.-H., Zhao, G., et al. 2015, Research in Astronomy and Astrophysics, 15, 1095 
[Marcy et al.(2014)]Marcy:2014 Marcy, G. W., Isaacson, H., Howard, A. W., et al. 2014, , 210, 20 
[Mayor et al.(2011)]Mayor:2011 Mayor, M., Marmier, M., Lovis, C., et al. 2011, arXiv:1109.2497 
[Morbidelli et al.(2012)]Morbidelli:2012 Morbidelli, A., Lunine, J. I., O'Brien, D. P., Raymond, S. N., & Walsh, K. J. 2012, Annual Review of Earth and Planetary Sciences, 40, 251 
[Mordasini et al.(2009)]Mordasini:2009 Mordasini, C., Alibert, Y., & Benz, W. 2009, , 501, 1139 
[Morton et al.(2016)]Morton:2016 Morton, T. D., Bryson, S. T., Coughlin, J. L., et al. 2016, , 822, 86 
[Mullally et al.(2015)]Mullally:2015 Mullally, F., Coughlin, J. L., Thompson, S. E., et al. 2015, , 217, 31 
[Mustill et al.(2017)]Mustill:2017 Mustill, A. J., Davies, M. B., & Johansen, A. 2017, , 468, 3000 
[Ofir et al.(2018)]Ofir:2018 Ofir, A., Xie, J.-W., Jiang, C.-F., Sari, R., & Aharonson, O. 2018, , 234, 9 
[Owen & Wu(2017)]Owen:2017 Owen, J. E., & Wu, Y. 2017, , 847, 29 
[Petigura et al.(2017)]Petigura:2017 Petigura, E. A., Howard, A. W., Marcy, G. W., et al. 2017, , 154, 107 
[Petigura et al.(2013)]Petigura:2013 Petigura, E. A., Howard, A. W., & Marcy, G. W. 2013, Proceedings of the National Academy of Science, 110, 19273 
[Petrovich(2015)]petrovich2015 Petrovich C., 2015, ApJ, 808, 120
[Pu & Wu(2015)]Pu:2015 Pu B., & Wu Y., 2015, ApJ, 807, 44
[Pu & Lai(2018)]Pu:2018 Pu, B., & Lai, D. 2018, arXiv:1801.06220 
[Raymond et al.(2004)]raymond2004 Raymond, S. N., Quinn, T., & Lunine, J. I. 2004, Icar, 168, 1
[Ren et al.(2016)]Ren:2016 Ren, A., Fu, J., De Cat, P., et al. 2016, , 225, 28 
[Rowe et al.(2014)]Rowe:2014 Rowe, J. F., Bryson, S. T., Marcy, G. W., et al. 2014, , 784, 45 
[Spalding & Batygin(2016)]Spalding:2016 Spalding, C., & Batygin, K. 2016, , 830, 5 
[Spalding et al.(2018)]Spalding:2018 Spalding, C., Marx, N. W., & Batygin, K. 2018, , 155, 167 
[Steffen et al.(2010)]Steffen:2010 Steffen, J. H., Batalha, N. M., Borucki, W. J., et al. 2010, , 725, 1226 
[Steffen et al.(2013)]Steffen:2013 Steffen, J. H., Fabrycky, D. C., Agol, E., et al. 2013, , 428, 1077 
[Thompson et al.(2017)]Thompson:2017 Thompson, S. E., Coughlin, J. L., Hoffman, K., et al. 2017, arXiv:1710.06758 
[Tremaine(2015)]tremaine2015 Tremaine, S. 2015, , 807, 157
[Tremaine & Dong(2012)]Tremaine:2012 Tremaine, S., & Dong, S. 2012, , 143, 94 
[Winn & Fabrycky(2015)]Winn:2015 Winn, J. N., & Fabrycky, D. C. 2015, , 53, 409 
[Xie et al.(2014)]Xie:2014 Xie, J.-W., Wu, Y., & Lithwick, Y. 2014, , 789, 165 
[Xie et al.(2016)]Xie:2016 Xie, J.-W., Dong, S., Zhu, Z., et al. 2016, Proceedings of the National Academy of Science, 113, 11431 
[Youdin(2011)]Youdin:2011 Youdin, A. N. 2011, , 742, 38 
[Zhao et al.(2012)]Zhao:2012 Zhao, G., Zhao, Y.-H., Chu, Y.-Q., Jing, Y.-P., & Deng, L.-C. 2012, Research in Astronomy and Astrophysics, 12, 723 
[Zhu et al.(2016)]Zhu:2016 Zhu, W., Wang, J., & Huang, C. 2016, , 832, 196 
[Zhu & Wu(2018)]Zhu:2018 Zhu, W., & Wu, Y. 2018, arXiv:1805.02660 



