A Flat Inner Disk Model as an Alternative to the Kepler Dichotomy in the Q1 to Q16 Planet Population

TL;DR

This study demonstrates that assuming flat inner planetary disks can explain the observed distribution of transiting planets around Kepler stars without requiring a dichotomous planet formation ability, challenging previous models.

Contribution

It introduces a flat inner disk model that reproduces Kepler planet multiplicity data without the Kepler dichotomy assumption.

Findings

Flat disk model matches observed planet multiplicity.

Average of about 2 planets per star in the specified period and radius range.

Contrasts with Solar System planet distribution.

Abstract

We use simulated planetary systems to model the planet multiplicity of Kepler stars. Previous studies have underproduced single planet systems and invoked the so called Kepler dichotomy, where the planet forming ability of a Kepler star is dichotomous, producing either few or many transiting planets. In this paper we show that the Kepler dichotomy is only required when the inner part of planetary disks are just assumed to be flared. When the inner part of planetary disks are flat, we reproduce the observed planet multiplicity of Kepler stars without the need to invoke a dichotomy. We find that independent of the disk model assumed, the mean number of planets per star is approximately 2 for orbital periods between 3 and 200 days, and for planetary radii between 1 and 5 Earth radii. This contrasts with the Solar System where no planets occupy the same parameter space.

Equations

$$\widehat{{\mathbf{N}}_{\mathbf{k}}}=[\widehat{N_1},\widehat{N_2},\widehat{N_3},..],$$

$$4100\text{K}<T_{\text{eff}}$$

$${\sigma }_{{\text{CDPP}}_{45}}$$

$$R_{\ast }$$

$$T_{\text{baseline}}$$

$$f_{\text{duty}}$$

$$3\text{days}<$$

$$1{\text{R}}_{\oplus }<$$

$$M$$

$$\widehat{{\mathbf{N}}_{\mathbf{k}}}=[631,127,45,9,3,1].$$

$$MES=\frac{0.84\delta (c+s\sqrt{\delta })}{{\sigma }_{\text{cdpp}}}\sqrt{n_{\text{tr}}}$$

$${\eta }_{\text{detect}}(MES)=\frac{1}{c^b\Gamma (b)}{\int }_0^{MES-\beta }{x}^{b-1}{e}^{-x/c}dx$$

$${\eta }_{\text{transit}}=\frac{R_{\ast }}{a},$$

$$\eta (P,R_p)={\eta }_{\text{transit}}\times {\eta }_{\text{detect}}$$

$$\frac{df}{dPd{R}_p}=C{F}_{0}g(P,R_p)$$

$$=$$

$$b=\frac{a}{R_{\ast }}\cos i\le 1$$

$$\xi =\frac{T_{\text{dur,in}}/P_{\text{in}}^{1/3}}{T_{\text{out,in}}/P_{\text{out}}^{1/3}},$$

$${\chi }^2=\sum _{k=1}^n\frac{(\widehat{{\mathbf{N}}_{\mathbf{k}}}-{\hat{\mathbf{N}}}_{\mathbf{k},\mathbf{ij}}{)}^2}{{\hat{\mathbf{N}}}_{\mathbf{k},\mathbf{ij}}}.$$

$${\chi }_{\text{combined}}^2\approx -2\sum _{m=1}^M\ln p_m$$

$$\frac{df}{dPd{R}_p}=C{F}_{0}g(P,R_p)$$

$$=$$

$${\int }_{R_{\text{min}}}^{R_{\text{max}}}{\int }_{P_{\text{min}}}^{P_{\text{max}}}Cg(P,R_p)dPd{R}_p=1$$

$$\ln (L)\propto \left[\sum _{i=1}^{N_{\text{pl}}}\ln \left(C{F}_{0}g(P,R_P)\right)\right]-N_{\text{exp}}$$

$$N_{\text{exp}}=C{F}_0{\int }_{R_{\text{min}}}^{R_{\text{max}}}{\int }_{P_{\text{min}}}^{P_{\text{max}}}\left[\sum _{j=1}^{N_{\ast }}{\eta }_j(P,R_p)\right]g(P,R_P)dPd{R}_p$$

$$R_{H,ij}={\left(\frac{m_i+m_j}{3M_{\ast }}\right)}^{1/3}\frac{a_i+a_j}{2}$$

$${\Delta }_{ij}=\frac{a_j-a_i}{R_{H,ij}}$$

$${\Delta }_{ij}+{\Delta }_{jk}<18$$

$$i_{\text{flare}}=⟨\theta ⟩+\varphi \cos (\Omega ).$$

$$i_{\text{flat}}=⟨\theta ⟩+\arcsin (Z_0/a)\cos (\Omega )$$

$$2R_{\ast }\approx v_{\text{orb}}{T}_{\text{dur}}$$

$$2\sqrt{R_{\ast }^2-b^2}=v_{\text{orb}}{T}_{\text{dur}}$$

$$2\sqrt{R_{\ast }^2-b_{\text{in}}^2}=T_{\text{dur,in}}{v}_{\text{orb,in}}$$

$$2\sqrt{R_{\ast }^2-b_{\text{out}}^2}=T_{\text{dur,out}}{v}_{\text{orb,out}}$$

$$v_{\text{orb,in}}\propto \frac{a_{\text{in}}}{P_{\text{in}}}\propto P_{\text{in}}^{-1/3}$$

Full text

††pagerange: A Flat Inner Disk Model as an Alternative to the Kepler Dichotomy in the Q1$-$Q16 Planet Population–E††pubyear: 2016

1 Introduction

The Kepler Q1$-$Q16 catalog (Mullally et al., 2015) uses 47 months of Kepler data collected from $\sim$190,000 stars. This has resulted in the detection of over 4,000 planet candidates orbiting $\sim$3,200 stars. While the majority of the $\sim$3,200 stars contain a single detected planet, transit signals from multiple planets have been detected around 656 of these stars. Comparisons to the architecture of the Solar System are limited, due to the relatively smaller periods and planetary radii that Kepler can efficiently sample.

1.1 The Kepler Dichotomy

The mutual inclination distribution between planets around Kepler stars has been well studied (Lissauer et al. (2011); Fang & Margot (2012); Figueira et al. (2012); Tremaine & Dong (2012); Johansen et al. (2012); Weissbein & Steinberg (2012); Ballard & Johnson (2016); Fabrycky et al. (2014), see Appendix E). The majority of of these studies show good agreement between simulated planets and the Kepler sample when the orbital planes of simulated planets are closely aligned. Specifically, when the mutual inclinations between planets are drawn from a Rayleigh distribution (a ‘flared’ planetary disk) with a mode of the flare angle between $\sim 1^{\circ}-5^{\circ}$.

In contrast to the agreement for mutual inclinations, some studies report a significant underproduction of simulated systems with a single detected planet (Lissauer et al., 2011; Johansen et al., 2012; Hansen & Murray, 2013; Ballard & Johnson, 2016). These studies underproduce the number of simulated stars with a single detected transiting planet by a factor of $\sim 3$.

The underproduction of simulated systems with a single detected transiting planet has led to the proposal of dichotomous planetary systems in the Kepler field, the so-called Kepler Dichotomy. One population of planetary systems is required to either suppress planet formation, or be “dynamically hot” (Hansen & Murray, 2013), where mutual inclinations between planets are increased, or where planets are more likely to be ejected from the system. For the host stars in these planetary systems, the probability of detecting multiple transiting planets is reduced, leading to a higher proportion of stars with a single detected transiting planet in this population.

Potential explanations for the dynamically hot planetary system population include dynamical instability caused by high mass planets (Johansen et al., 2012; Lai & Pu, 2016), instability or suppressed planet formation caused by stellar binaries (Ballard & Johnson, 2016), varying surface density profiles and disk masses (Moriarty & Ballard, 2015), varying strengths of gas depletion or spin-orbit misalignment between the star and planet (Spalding & Batygin, 2016). Ballard & Johnson (2016) show that to account for the excess of detected single-planet transiting systems around M dwarfs, these stars with a reduced probability of multiple transiting planets need to account for $\sim 55$ per cent of M dwarfs in the Kepler field.

1.2 Detected Transiting Planets

We define the true planetary system multiplicity vector $\bm{N_{k}}$ as the number of stars which are host to $k$ planets. For the Kepler mission (and any transit survey), the observed vector $\bm{\hat{N_{k}}}$ will be significantly lower than the true $\bm{N_{k}}$, due to the low probability of planets transiting their host star, and since Kepler can only efficiently detect planets across a small fraction of parameter space. The observed planetary system multiplicity vector, hereafter simply referred to as the multiplicity vector, is given by

[TABLE]

where $\hat{N_{1}}$, $\hat{N_{2}}$ and $\hat{N_{3}}$ are the the number of stars with 1, 2 and 3 detected transiting planets respectively, and so on. For the Kepler Q1$-$Q16 catalog (see Section 2.1), $1\leq k\leq 6$ and $\bm{\hat{N_{k}}}=[$2608, 413, 141, 52, 18, 3].

1.3 Mutual Inclinations

For two or more planets in the same planetary system, the mutual inclination between those planets is defined as the angle between their orbital planes. The probability of multiple planets transiting the same star is non-negligible for small mutual inclinations only, generally on the order of a few degrees. Planets in the system with larger mutual inclinations, relative to the transiting planets, require alternative detection methods.

In general, the true inclination of the orbital plane of a transiting planet cannot be determined from a transit lightcurve alone. The transit method is only sensitive to the line-of-sight component of the inclination $i$ (Figure 10). In Figure 1, the distribution of the detectable inclination component for a set of simulated planets is shown in the bottom panels. The orthogonal component of inclination, typically not detectable by the transit method, represents the y-axis of the middle panels. The true mutual inclination between a pair of planets is given by $\sqrt{\Delta\theta^{2}+\Delta\theta_{\text{(y-z)}}^{2}}$.

1.4 An Alternative Disk Model

In the studies mentioned in Section 1.1, the true mutual inclinations between simulated planets are drawn from a Rayleigh distribution with mode $\sigma_{\Delta\phi}$. The Rayleigh distribution is composed of two Gaussian distributed components, with standard deviations equal to $\sigma_{\Delta\phi}$. We can visualize the inclination distribution by considering one of these Gaussian components, i.e. viewing systems edge-on at an arbitrary plane perpendicular to the invariable plane, as in the top panels of Figure 1. Rayleigh distributed mutual inclinations represent a ‘flared disk’ model, where a planet’s height above the invariable plane111 The mode of the Rayleigh distribution of inclinations relative to an invariable plane $\sigma_{\phi}$, is related to the Rayleigh distribution of mutual inclinations $\sigma_{\Delta\phi}$, by $\sigma_{\phi}\approx\sigma_{\Delta\phi}/\sqrt{2}$.

tends to increase with increasing semi-major axis. Planet inclinations relative to the invariable plane do not depend on semi-major axis.

In this paper, we use a ‘flat disk’ model, where a planet’s height above the invariable plane does not depend on semi-major axis, and planet inclinations relative to the invariable plane tend to decrease with increasing semi-major axis (as seen in the right-side panels of Figure 1). Hansen & Murray (2013) tested the in situ assembly of close in planets, and found that planets with small semi-major axes tended to have larger inclinations, particularly $<0.1$ AU.

We apply this flat disk model to the typical semi-major axis space probed by Kepler, i.e. the interior part of planetary disks, as shown in Figure 2. In general, this represents planets with semi-major axes much less than the semi-major axes of the inner Solar System planets. We improve on previous modeling efforts by removing the flared disk assumption. We show that for a flat inner planetary disk there is no need to invoke a dichotomous planetary system population, where one population of host stars have a decreased probability of hosting multiple transiting planets.

In Section 2 we define our stellar and planetary samples, based on minimizing false positives and false negatives. In Section 3 we estimate the transit and detection completeness across our parameter space. In Section 4 we estimate the underlying orbital period and planet radius distributions. We outline the process of producing model planetary populations in Section 5. In Section 6, we compare the simulated detections in our model systems to the Kepler Q1$-$Q16 candidates, for both flat and flared disk models. In Section 7 we discuss the results from our model planet populations, including estimates for the mean number of planets per star within our parameter space.

2 Sample selection

When we generate a model planetary system, the stellar properties for that system are assigned from a random Kepler star in our stellar sample. We produce our input stellar sample in the following way.

We begin with the 198,917 stars from the Kepler Q16 stellar catalog222 http://exoplanetarchive.ipac.caltech.edu/cgi-bin/TblSearch/nph-tblSearchInit?app=ExoTbls&config=keplerstellar. We limit our sample to low-noise Solar type stars, similar to the majority of previous studies mentioned in Section 1.1. We apply the following cuts to the input catalog:

[TABLE]

where $T_{\text{eff}}$ and $R_{*}$ are the stellar effective temperature and radius respectively. The 4.5 hour CDPP (combined differential photometric precision, Christiansen et al. (2012)) of the star, a measure of the combined instrumental and stellar noise, is given by $\sigma_{\text{CDPP}_{45}}$. $T_{\text{baseline}}$ is the timespan of observations for each star and $f_{\text{duty}}$ is the fraction of valid observations over $T_{\text{baseline}}$. Note that the combination of $T_{\text{baseline}}$ > 1000 days and $f_{\text{duty}}$ > 0.6 generally ensures at least 3 transits for orbital periods up to 200 days.

The above stellar cuts, in addition to removing stars without a reported mass, results in our input stellar sample of 63,128 stars. In later sections, the stellar properties of each simulated planetary system are assigned from a randomly drawn star in this sample.

To minimize the detection incompleteness and false-positives in our observed planet sample, to which we will compare our simulations, we select only those planets with a high pipeline detection efficiency. We begin with the 1,338 Kepler Objects of Interest (KOIs) labeled ‘candidate’ by the Q1$-$Q16 pipeline333 http://exoplanetarchive.ipac.caltech.edu/cgi-bin/TblView/nph-tblView?app=ExoTbls&config=koi. An additional 4,401 KOIs are labeled ‘not dispositoned’. We update these dispositions using the Kepler Q17 catalog for reference. This results in 2,812 KOIs changing from ‘not dispositioned’ to ‘candidate’.

These planets form our initial sample of 4,150 planet candidates from the Q1$-$Q16 catalog (Mullally et al., 2015). To conform with our input stellar sample, we remove planets around host stars outside of our stellar parameter space defined by Equation 2. This reduces our sample of observed planets from 4,150 to 1,790.

We set an upper orbital period limit of 200 days to avoid the increase in false-positives towards the Kepler orbital period of $\sim\,$372 days (Mullally et al., 2015), and to remain consistent with the Kepler pipeline completeness calculations (Christiansen et al., 2015).

The Kepler pipeline is known to have an increasing false-negative rate with decreasing orbital period for orbital periods $\lesssim 3$ days. This is largely due to the pipeline harmonic filter, which can remove transit signals which are on the same timeframe as the expected stellar noise (Christiansen et al., 2015). In addition, a small fraction of the fitted planetary radii for planets with orbital periods $\lesssim$ 10 days can be significantly lower than the true planet radius, diluting the transit signals for some of these planets. We choose a orbital period lower limit of 3 days, in order to retain a sample of Kepler stars with $\geq 4$ transiting planets.

The Kepler pipeline reports a summary statistic for the strength of a transit detection, the Multiple Event Statistic (MES). A lower limit planet radius of $1\;R_{\oplus}$ and a lower limit MES of 10 are chosen since false-positives are dominated by low MES ($\lesssim 8$) detections (Mullally et al., 2015). An upper planet radius limit of $5\;R_{\oplus}$ is chosen to avoid increasing false-positives with planet size, and since the mass-radius relation becomes degenerate for larger planetary radii.

To summarize, we only retain the Kepler Q1$-$Q16 candidates which meet the following criteria:

[TABLE]

This results in our observed sample of 1,077 candidates in 816 planetary systems, within the parameter space outlined in Equations 2 and 3. The observed planetary system multiplicity vector $\bm{\hat{N_{k}}}$ (Equation 1) for our parameter space is given by

[TABLE]

3 Transit and detection efficiency

When attempting to estimate the underlying multiplicity vector $\bm{N_{k}}$ given the observed $\bm{\hat{N_{k}}}$ (Equation 1), there exists a degeneracy between the underlying multiplicity distribution and the underlying mutual inclination distribution. For example, an observed $\bm{\hat{N_{k}}}$ could be reproduced by systems which contain many planets with a large dispersion in mutual inclinations, or by systems containing fewer planets but with a small mutual inclination dispersion. These two underlying distributions must be modeled simultaneously.

We estimate the underlying inclination and multiplicity distributions of Kepler systems in the Q1$-$Q16 catalog, within the period and radius parameter space where Kepler can more reliably detect transiting planets (given by Equations 2 and 3). We produce sets of $10^{6}$ simulated planetary systems across a grid of inclination and multiplicity distributions. For each set of model assumptions, we estimate the probability of Kepler detecting each simulated planet. By comparing the $\bm{\hat{N_{k}}}$ for a set of simulated systems to the $\bm{\hat{N_{k}}}$ for the Q1$-$Q16 Kepler catalog, we can estimate the underlying architecture between planetary orbital planes and the distribution of the number of planets per star.

3.1 Pipeline detection efficiency

Detection incompleteness and false-positives are important issues when comparing the detected planets around simulated and observed stars. Previous studies did not have the benefit of the Kepler pipeline detection completeness provided by transit injection and recovery experiments (Christiansen et al., 2015), shown in Figure 3.

We follow the approximation of the pipeline MES by (Burke et al., 2015), which includes a limb-darkening approximation and accounts for nonzero impact parameters.

[TABLE]

where $c=1.0874$ and $s=1.0187$ for G dwarfs, $\delta=(R_{p}/R_{*})^{2}$, and $n_{\text{tr}}$ is the number of transits. Values of $\sigma_{\text{cdpp}}$ are reported for 14 different transit durations from 1.5 hours to 15 hours for each Kepler star (Burke et al., 2015). The $\sigma_{\text{CDPP}}$ value chosen for Equation 5 is interpolated from the 14 reported CDPP values, to match the transit duration of the planet.

The number of transits for a planet is estimated by $n_{\text{tr}}=(T_{\text{baseline}}\times f_{\text{duty}})/P$, where $P$ is the planet period, $T_{\text{baseline}}$ is the total observing time for the Q1$-$Q16 catalog ($\sim\,$1426 days) and $f_{\text{duty}}$ is the duty cycle for the observed star; the fraction of valid observations over the observing baseline. Note that $T_{\text{baseline}}$ and $f_{\text{duty}}$ are reported for each star, accounting for systematics such as the differences in CCD detectors and pixels.

We define $\eta_{\text{detect}}$ as the Kepler pipeline completeness, shown in Figure 3. The pipeline completeness as a function of the multiple event statistic is approximately represented by the $\Gamma$ cumulative distribution function

[TABLE]

where $\Gamma$ is the Gamma function. For our sample of FGK dwarfs, $b=4.35$, $c=1.05$ and $\beta=4.093$ (Christiansen et al., 2015).

To calculate the total probability of transit detection $\eta(P,R_{p})$, we must also take into account the geometric transit probability $\eta_{\text{transit}}$ of a planet,

[TABLE]

where $a$ is the semi-major axis of the planet. The product of these two equations gives the total transit detection probability $\eta(P,R_{p})$, the probability of the planet transiting $\eta_{\text{transit}}$, and the probability of the transiting planet being detected by the Kepler Q1$-$Q16 pipeline $\eta_{\text{detect}}$,

[TABLE]

Given an input star from our stellar sample (Equation 2) and using Equations 5 to 8, we can estimate the total transit and detection completeness of a simulated planet with planet parameters of period, radius and transit duration. Across our planetary parameter space used in this paper (Equation 3), we estimate the mean total transit and detection completeness $\langle\eta(P,R_{p})\rangle$, by taking the mean value of $\eta(P,R_{p})$ at each grid point over all stars in our stellar sample. This is shown in Figure 4, where $\langle\eta(P,R_{p})\rangle$ ranges from $\sim 0$ to a maximum of $\sim 0.1$. Transit and detection probabilities $>5$ per cent only exist for planets with orbital periods $\lesssim 8$ days. It can be seen that the pipeline detection probability becomes important for planetary radii less than $2.5$ Earth radii.

4 Underlying planet distributions

Our simulated results of the number of stars with $k$ detected transiting planets is reliant on an input planet radius and orbital period distribution, which has been one of the primary goals of the Kepler mission. The planet radius distribution is often modeled as a broken power law (Youdin, 2011; Howard et al., 2012; Burke et al., 2015), with a logarithmic plateau at $\lesssim 2.5R_{\oplus}$. This logarithmic plateau is also seen when the pipeline efficiency is probed using transit injection and recovery experiments (Petigura et al., 2013).

For orbital periods between 50 and 300 days, a single power law is sufficient to describe the orbital period distribution (Burke et al., 2015). Our parameter space includes planets with orbital periods less than 50 days, where the transit and detection completeness is more dynamic, particularly for periods $\lesssim 15$ days (Figure 4). For this parameter space, a single power law is not sufficient, and we model the orbital period distribution as a broken power law.

The planet radius and orbital period distributions are combined into a planet distribution function (PLDF), in this case composed of a broken power law for the distribution of orbital periods, and a broken power law for the distribution of planetary radii.

Our PLDF has 7 free parameters, $F_{0}$, $\beta_{1}$, $\beta_{2}$, $P_{\text{brk}}$, $\alpha_{1}$, $\alpha_{2}$, $R_{\text{brk}}$, where $F_{0}$ is the number of planets per star within our parameter space, and $R_{\text{brk}}$ and $P_{\text{brk}}$ are the transition points between the two power laws for the planet radius and orbital period respectively,

[TABLE]

where the power law exponents $\alpha_{1}$, $\alpha_{2}$ and $\beta_{1}$, $\beta_{2}$ represent the exponents for the orbital period and the planet radius distribution respectively, either side of the power law breaks.

For each set of model parameters in an underlying planet distribution function, an expected number of planet detections is computed by convolving the planet distribution function with $\langle\eta(P,R_{p})\rangle$ (Figure 4). The number of expected detections for an underlying planet distribution function is then compared to the number of Kepler Q1$-$Q16 detections by maximizing the Poisson likelihood of the PLDF. The maximum likelihood derivation for our PLDF (Equation 9) is shown in Appendix B, with maximum likelihood parameters of

Maximum likelihood PLDF parameters

$F_{0}$ $\beta_{1}$ $\beta_{2}$ $P_{\text{brk}}$ $\alpha_{1}$ $\alpha_{2}$ $R_{\text{brk}}$

0.852 1.007 -0.932 15.332 -1.168 -4.906 2.740

indicating breaks in the power law distributions at $\sim$15 days and $\sim$2.7 $R_{\oplus}$ for orbital periods and planetary radii respectively. The break in the orbital period distribution corresponds to a peak in the distribution, whereas the break in the planet radius distribution corresponds to the logarithmic planet plateau for $R_{p}\lesssim 2.7$ $R_{\oplus}$. These results appear to be consistent with Foreman-Mackey et al. (2014), where breaks in the logarithmic orbital period and planet radius rates are indicated at $\sim 20$ days and $\sim 2-3$ Earth radii respectively.

We can marginalize our maximum likelihood PLDF in terms of orbital period and planet radius. This is shown in Figures 5 and 6 respectively, where the thick red lines represent our marginalized maximum likelihood PLDF, which is the estimated underlying planet distribution. The dashed red lines indicate the corresponding transit detected distribution, after applying the mean total transit and detection probability $\langle\eta(P,R_{p})\rangle$ for the stars in our sample.

When producing model planetary systems in Section 5, we assign orbital periods and planetary radii by drawing randomly from the maximum likelihood underlying distributions, shown in Figures 5 and 6.

4.1 Parameterizing Planet multiplicity

For a set of model planetary systems, we need to assume a distribution for the inherent number of planets per star within our parameter space, $N_{pl}$. In this paper we trial two different $N_{pl}$ distributions.

The first trial distribution is a modified Poisson distribution $N_{pl,\text{Poi}}$ (Fang & Margot, 2012). Each star is assigned a random number of planets, drawn from a Poisson distribution with mean $\mu$. Stars drawn with zero planets are redrawn from the same Poisson distribution, until all $10^{6}$ model planetary systems are populated with planets, resulting in a final mean $\geq\mu$.

The second trial distribution is a modified Exponential distribution $N_{pl,\text{exp}}$ (Gaidos et al., 2016) and is produced in the same way as $N_{pl,\text{Poi}}$, except that stars are assigned a number of planets drawn from an exponential distribution with mean $\mu$. The mode of the exponential distribution is always 0, resulting in a natural tendency for more planetary systems to contain a single transiting planet rather than multiple transiting planets. It has been shown that when an exponential distribution is used to model the inherent number of planets per star, no Kepler dichotomy is required (Gaidos et al., 2016). We include this trial distribution for comparative purposes.

5 Producing model planetary systems

In our simulations, we assume two different planetary disk models, and two different distributions for the number of planets per star $N_{pl}$, resulting in simulations with four unique combinations of model assumptions. For a given set of model assumptions, we generate populations of planetary systems across a grid. The mean of the number of planets per star $\mu$ ranges from 0.5 to 3.5 in steps of 0.1. For the flared disk model, the mode of the Rayleigh distributed mutual inclinations ranges from 0 to 5 degrees in steps of 0.1 degrees. Similarly for the flat disk model, the standard deviation of the height above the invariable plane ranges from 0 to 5 $R_{*}$ in steps of 0.1 $R_{*}$.

This results in a total of 1,500 grid points for each set of model assumptions, with $10^{6}$ model planetary systems generated at each grid point. Each model planetary system is produced as follows:

{easylist}\ListProperties

(Space*=0.1cm) & A random star is chosen from our sample of 63,128 Kepler stars outlined in Section 2, and its mass and radius are assigned to the star in the model system.

The angle to the system’s invariable plane relative to the observer $\langle\theta\rangle$ (Figure 10), is chosen from a random point on a sphere: $0\leq\cos\left[\langle\theta\rangle\right]\leq 1$.

The number of planets in the system $n_{i}$ is drawn randomly, according to the assumed $N_{pl}$ distribution from Section 4.1, with a mean value based on the current grid point.

The radii of the $n_{i}$ planets are drawn from the underlying distribution in Section 4 ($1\;\text{R}_{\oplus}<R_{p}<5\;\text{R}_{\oplus}$), and converted to their corresponding masses444$M_{p}\approx\left(R_{p}/a\right)^{b}$, where a$$\sim1.11 and b$$\sim2.41..

The periods of the $n_{i}$ planets are drawn randomly from the underlying distribution in Section 4 ($3\;\text{days}<P<200\;\text{days}$), and converted to their corresponding semi-major axes, using the stellar properties of the assigned Kepler star.

The dynamical stability of the system is estimated by testing the stability of sets of 3 sequential planets, or pairwise if $n_{i}=2$. If any set of planets in the system is deemed unstable, the system is labeled unstable and new planet periods for all $n_{i}$ planets are redrawn as in step 5. See Section C.1 for a complete description of estimating the stability of a system, including termination criteria.

The inclinations of the planetary orbital planes relative to the observer, are determined from $\langle\theta\rangle$ (step 2) and the assumed planetary disk model (Section 1.4), and the parameter value at the current grid point. See Appendix C.2 for a complete description of how planet inclinations are assigned for flat and flared disk models.

The above steps generate the $10^{6}$ model systems according to the assumed disk model, the assumed planet multiplicity distribution, the current grid point parameters, and the underlying planet period and radius distributions. The final step is to estimate which simulated planets would be detected by the Kepler Q1$-$Q16 pipeline, and compare this detected sample to the observed Q1$-$Q16 detections.

5.1 Determining transiting and detected planets

Once planetary inclinations are assigned, the model system is complete and we test for transiting planets. We define a transiting planet by its impact parameter $b$, where a planet is defined to transit if:

[TABLE]

where $R_{*}$, $a$ and $i$ were determined from steps 1, 5 and 7 respectively.

For each simulated transiting planet, we estimate the Multiple Event Statistic ($MES$, Equation 5). The $MES$ is dependent on stellar properties, along with the planet’s orbital period, radius and transit duration. Circular orbits are assumed when estimating transit durations. The planet’s $MES$ is then used to estimate the pipeline detection efficiency $\eta_{\text{detect}}$ (Equation 6). For each simulated transiting planet, a uniform random number $Y_{m}$ is drawn between [math] and $1$. A simulated planet is labeled as detected if it transits, and if its pipeline detection efficiency $\eta_{\text{detect}}>Y_{m}$. All simulated planets which meet this criteria are added to the detected planet sample for the grid point, ${\bm{X_{\bm{ij}}}}$, where $i$ and $j$ represent the current grid point.

6 comparing Simulated and observed planet detections

The simulations outlined in Section 5 were performed across a grid for the 4 sets of model assumptions. For each grid point, the simulated planet detections ${\bm{X_{\bm{ij}}}}$ are used to generate two distributions, the system multiplicity vector $\bm{\hat{N}_{\bm{k,ij}}}$ (Equation 1), and the distribution of orbit normalized transit duration ratios $\bm{\xi_{\bm{ij}}}$ (Steffen et al., 2010). Unlike the $\bm{\hat{N}_{\bm{k,ij}}}$ distribution, the $\bm{\xi_{\bm{ij}}}$ distribution only consists of model systems with two or more detected transiting planets. For a pair of planets orbiting the same star,

[TABLE]

where $T_{\text{dur}}$ and $P$ are the transit durations and the periods for the inner and outer planets, given by the subscripts in and out respectively. For each unique planet pair in a system, $\xi$ is calculated, giving $n_{i}(n_{i}+1)/2$ values of $\xi$ for a star with $n_{i}$ planets. For each grid point we generate the ensemble $\bm{\xi}_{ij}$ distribution by calculating the $\xi$ value for each unique pair of simulated transit detections orbiting the same star, across all $10^{6}$ model systems. For a deeper discussion of $\xi$, see Appendix D.

The $\bm{\hat{N}_{\bm{k,ij}}}$ and $\bm{\xi_{\bm{ij}}}$ distributions are compared to the $\bm{\hat{N_{k}}}$ and $\bm{\xi}$ distributions of the Kepler Q1$-$Q16 candidates, and are used to assess the goodness of fit at each grid point. We perform a $\chi^{2}$ goodness of fit test (Equation 12) comparing the simulated $\bm{\hat{N}_{\bm{k,ij}}}$ to the observed $\bm{\hat{N_{k}}}$ for our parameter space (Equation 4). We scale $\bm{\hat{N}_{\bm{k,ij}}}$ such that $\sum\bm{\hat{N}_{\bm{k,ij}}}=\sum\bm{\hat{N_{k}}}$. To compensate for the poor quality of the $\chi^{2}$ test with low cell counts, values less than 5 are merged into their adjacent cells.

[TABLE]

Similarly, we perform a two-sample Kolmogorov$-$Smirnov (KS) test between the simulated $\bm{\xi}_{ij}$ distribution at each grid point, and the $\bm{\xi}$ distribution of the observed Q1$-$Q16 Kepler candidates within our parameter space.

6.1 Flared Disk and Poisson distributed planets per star

In the top panel of Figure 7, $\bm{\hat{N}_{\bm{k,ij}}}$ is compared to $\bm{\hat{N_{k}}}$ at each grid point, under the assumption of a flared planetary disk and a Poisson distributed number of planets per star. The $\chi^{2}$ values are represented by the 1$\sigma$, 2$\sigma$ and 3$\sigma$ values relative to the best fit. As expected, no good fit is found to the $\bm{\hat{N_{k}}}$ distribution, as is the case in the majority of previous studies (Lissauer et al., 2011; Johansen et al., 2012; Ballard & Johnson, 2016; Gaidos et al., 2016).

The bottom panel displays the resulting p-values from the KS test between the $\bm{\xi}_{ij}$ and $\bm{\xi}$ distributions. The orbital normalized transit duration ratios favor mutual inclinations with a mode between $1.5-4$ degrees, consistent with all previous studies shown in Appendix E. The mean number of planet per star cannot be determined from the $\bm{\xi_{\bm{ij}}}$ distribution alone.

It is clear from Figure 7 that the best-fit regions (dark red) of the two tests do not appear consistent. Comparing multiplicity vectors favors near coplanar mutual inclinations, with a mode $\lesssim 1^{\circ}$ (top panel). However, modes $\lesssim 1.5^{\circ}$ are ruled out by comparing orbit normalized transit duration ratios (bottom panel).

6.2 Flat Disk and Poisson distributed planets per star

For the set of simulations with a Poisson distributed $N_{pl}$ and a flat disk model, the two tests appear more consistent (Figure 8). Unlike for the flared disk model, the best fit $\bm{\hat{N}_{\bm{k,ij}}}$ is a good match to $\bm{\hat{N_{k}}}$, giving a $\chi^{2}/\text{dof}$ of 0.9, indicating that no Kepler dichotomy is required.

While Gaussian disk thicknesses $>2\;R_{*}$ are supported by comparing multiplicity vectors, comparing orbit normalized transit durations refines the disk thickness to $1\;R_{*}\lesssim Z_{0}\lesssim 2\;R_{*}$. There is significant overlap between the two tests within this region.

6.3 Combining independent tests

The results from comparing the $\bm{\hat{N_{k}}}$ and $\bm{\xi}$ distributions can be combined in order to estimate the overall best-fit parameters, for a given set of model assumptions. The $p$ values from each test are combined using Fisher’s method into a single test statistic,

[TABLE]

where $M$ is the number of tests combined and $p_{m}$ is the p-value of the $m$th test. The degrees of freedom is given by $2M$, where in this case $M=2$. We use this combined statistic to produce a probability grid, $P_{\text{combined}}$, across the parameter space for each set of model assumptions. $P_{\text{combined}}$ is derived from the likelihood $P_{ij}\propto\exp(-\chi_{\text{combined},ij}^{2}/2)$ and the requirement $\sum P_{ij}=1$.

Figure 9 displays the probability grids for each set of model assumptions, along with the best-fit point and the $1\sigma$ and $2\sigma$ probability contours. Panel $a)$ of Figure 9 combines the tests of Figure 7, and Panel $b)$ combines the tests of Figure 8. A similar process is involved for panels $c)$ and $d)$, where the number of planets per star $N_{pl}$ is drawn from an exponential distribution. The intermediate figures for these two panels are not displayed for succinctness.

7 Results and Discussion

7.1 Flared Disk Model

Our result for a flared disk with a Poisson distributed number of planets per star, appears to be compatible with the majority of previous analyses. We find a mean number of planets per star $\mu=2.0^{+0.3}_{-0.2}$, over our parameter space, [3 days $<P<200$ days] and [1 $R_{\oplus}<R_{p}<5\;R_{\oplus}$]. For a similar orbital period and planet radius parameter space, Gaidos et al. (2016) find $\mu=2.2\pm 0.3$ for Kepler M dwarfs. Fang & Margot (2012) report $\mu\sim$1.5 for $R_{p}>1.5\;R_{\oplus}$, where the reduction in $\mu$ likely comes from the exclusion of planets with radii between $1.0\;R_{\oplus}<R_{p}<1.5\;R_{\oplus}$.

Similarly, we find the mode of the Rayleigh distributed mutual inclinations is given by $\sigma_{i}=2.3^{+0.9}_{-0.4}$ degrees, consistent with the bulk of previous results with $\sigma_{i}\sim 2^{\circ}$ (Lissauer et al., 2011; Tremaine & Dong, 2012; Figueira et al., 2012; Fang & Margot, 2012; Johansen et al., 2012; Fabrycky et al., 2014; Ballard & Johnson, 2016; Gaidos et al., 2016).

We are unable to achieve a good match to the Kepler Q1$-$Q16 detections for a flared disk model, contrary to the reported result by Fang & Margot (2012), where a flared disk model reproduced the $\bm{\hat{N_{k}}}$ distribution without the need for a Kepler dichotomy. The discrepancy likely comes from the unique $N_{pl}$ distribution chosen by Fang & Margot (2012), a “bounded uniform” distribution. The bounded uniform distribution is produced by first choosing a maximum number of planets $n_{\text{i,max}}$ from a Poisson distribution, then choosing the number of planets in the system $n_{i}$ from a uniform distribution between 1 and $n_{\text{i,max}}$.

It has previously been shown that the Kepler $\bm{\hat{N_{k}}}$ distribution can be matched without the need for a Kepler dichotomy, when $N_{pl}$ is drawn from an exponential distribution (Gaidos et al., 2016). We find that the Kepler sample is consistent with $\sigma_{i}=2.4^{+0.9}_{-0.5}$ degrees and $\mu=1.6^{+0.3}_{-0.2}$ drawn from an exponential distribution. Here we disagree with Gaidos et al. (2016), who preferred near-coplanar mutual inclinations. While we also achieve good fits for near-coplanar orbital planes, comparing $\bm{\xi}$ distributions strongly rules out mutual inclinations $\lesssim 1.4^{\circ}$. This illustrates the importance of modeling both $\bm{\hat{N_{k}}}$ and $\bm{\xi}$ distributions, where Gaidos et al. (2016) only modeled the $\bm{\hat{N_{k}}}$ distribution.

7.2 Flat Disk Model

For a flat planetary disk model (Dullemond & Monnier, 2010), a good fit to the Kepler candidates can be achieved when the number of planets per star is drawn from both a Poisson or exponential distribution. That is, independent of the $N_{pl}$ distribution chosen, the flat disk model removes the need for a Kepler dichotomy. When $N_{pl}$ is drawn from a Poisson distribution, we find $Z_{0}=1.6^{+0.6}_{-0.4}\;R_{*}$ and $\mu=2.4^{+0.6}_{-0.4}$. Notably, the mean number of planets per star $\mu$ is consistent between the assumed planetary disk models.

We use a flat planetary disk model with a Gaussian disk thickness $Z_{0}$. We can compare this value to the inner Solar System (Figure 2). For the inner Solar System planets, $Z_{\text{max}}\approx 8\;R_{*}$, giving $Z_{0}\approx 5\;R_{*}$, where $Z_{0}\approx 2\,Z_{\text{max}}/\pi$. This is significantly larger than our derived value of $Z_{0}=1.6^{+0.6}_{-0.4}\;R_{*}$ for our sample of closely packed Kepler systems. This may give some indication of the flat disk model’s applicability at larger semi-major axes, or may be reflective of the different parameter spaces probed.

8 Summary and Conclusion

We estimate the inherent orbital period and planet radius distributions for the Kepler Q1$-$Q16 catalog, within the parameter space [3 days $<P<200$ days] and $[1\;R_{\oplus}<R_{p}<5\;R_{\oplus}]$. We find that both distributions are well described by broken power laws, with breaks occurring at $\sim$15 days and \sim$$2.7\;R_{\oplus}. These inherent distributions are used to populate model planetary systems for flat and flared planetary disk models, and for the number of planets per star $N_{\text{pl}}$ drawn from Poisson and exponential distributions.

We confirm that a flared planetary disk model with $N_{\text{pl}}$ drawn from a Poisson distribution, is not consistent with the Kepler detections. We also confirm that Kepler detections are well matched when $N_{\text{pl}}$ is drawn from an exponential distribution, without the need to invoke a dichotomous planetary system population. In this paper we use a flat inner planetary disk model, where planets with smaller periods tend to have larger inclinations. When a flat rather than a flared planetary disk model is assumed, model systems are consistent with Kepler detections, without the requirement of a Kepler dichotomy, and independent of the chosen $N_{\text{pl}}$ distribution.

We find that the mean number of planets per star $\mu$ is largely model independent, $\sim$2.0 when $N_{\text{pl}}$ is drawn from a Poisson distribution, and $\sim$1.6 when $N_{\text{pl}}$ is drawn from an exponential distribution, for [3 days $<P<200$ days] and [$1\;R_{\oplus}<R_{p}<5\;R_{\oplus}$]. This contrasts with the Solar System where there are 0 planets within this parameter space.

Similarly, we find for a flared planetary disk model, mutual inclinations are distributed with a mode \sim$$2.2^{\circ}. For a flat planetary disk model, the Gaussian disk thickness $Z_{0}\sim\,$1.5 $R_{*}$, much lower than the $\sim$5 $R_{*}$ of the inner Solar System.

8.1 The Kepler Dichotomy

The underproduction of model systems with a single detected transiting planet has been well studied. This has lead to the invocation of a dichotomous planetary system population, where one population suppresses the number of detected transiting planets, resulting in a higher likelihood of producing a single detected transiting planet. Many physical explanations for the existence of the dichotomy have been put forward (Johansen et al., 2012; Weissbein & Steinberg, 2012; Hansen & Murray, 2013; Moriarty & Ballard, 2015; Lai & Pu, 2016; Spalding & Batygin, 2016).

Dawson et al. (2016) generated sets of planetary systems with various gas depletion factors using N-body simulations of planetary embryos. No set of simulations was a good match to the period ratio, $\Delta$ (Equation 19), planet multiplicity and $\xi$ distributions of the observed Kepler sample. Some improvement was found when simulated planetary systems were allowed to be a mix of “dynamically hot” and “dynamically cold systems”. However, this improvement becomes less pronounced when taking into account the partial correlations between these distributions, particularly between $\xi$ and $\Delta$.

It has also been shown that the requirement of the dichotomy is not robust to the assumed distribution for the number of planets per star (Fang & Margot, 2012; Gaidos et al., 2016). This is confirmed in this paper, and in addition, we show that a planetary system dichotomy is also not required for a flat inner planetary disk model. This result is independent of the choice of distribution for the number of planets per star $N_{p}$. We emphasize that we apply the flat planetary disk model only to the short period range of Kepler candidates.

Of the sets of model assumptions explored in this paper, the need for a Kepler dichotomy only exists for a flared inner planetary disk, with the number of planets per star drawn from a Poisson distribution.

The Kepler dichotomy describes the apparent need for a dichotomous planetary system population, with respect to a star’s probability of producing multiple transiting planets. We show that the Kepler dichotomy is only required under specific model assumptions. Specifically, when the inner part of a planetary disk is assumed to be flared, while also requiring the number of planets per star to be Poisson distributed. When removing either or both of these assumptions, the need for a Kepler dichotomy disappears.

Appendix A Inclination angles of planetary orbital planets

There are a number of different angles used in the literature which have all been referred to as the planet inclination. Where we have used an inclination angle, we have attempted to be as explicit as possible. The below figure illustrates different inclination angles used throughout the paper.

Appendix B Planet distribution function

Our planet distribution function (PLDF) has 7 free parameters, $F_{0}$, $\beta_{1}$, $\beta_{2}$, $P_{\text{brk}}$, $\alpha_{1}$, $\alpha_{2}$, $R_{\text{brk}}$.

[TABLE]

where $F_{0}$ is the number of planets per star within our parameter space, $R_{\text{brk}}$ and $P_{\text{brk}}$ are the transition points between the two power laws for the planet radius and orbital period respectively. The normalization constant $C$ is calculated from the requirement

[TABLE]

where the integration limits $R_{\text{min}}$, $R_{\text{max}}$, $P_{\text{min}}$ and $P_{\text{max}}$ are given in Equation 3.

We follow Youdin (2011) and Burke et al. (2015) by implementing a Poisson likelihood for our PLDF. By maximizing this likelihood we can obtain best-fit parameters for our model.

[TABLE]

where $N_{\text{exp}}$ is the expected number of planet detections for the set of model parameters, and is given by

[TABLE]

where $\eta_{j}(P,R_{p})=\eta_{\text{detect}}\times\eta_{\text{transit}}$ is the combined transit and pipeline detection efficiency of the $j$th star for the specified period and radius. The pipeline detection efficiency $\eta_{\text{detect}}$ is given by Equation 6 and the transit probability $\eta_{\text{transit}}=R_{*}/a$, where $a$ is the semi-major axis.

We calculate $\left[\sum_{j=1}^{N_{*}}\eta_{j}(P,R_{p})\right]$ for a grid in orbital period and planet radii, in bins of 1.5 days and 0.05 $R_{\oplus}$ respectively. For each grid point, we sum over all stars in our sample. The mean combined transit and pipeline detection efficiency $\langle\eta(P,R_{p})\rangle$ can then be found be dividing this term by the number of stars in our sample, $N_{*}$.

Appendix C Simulated planetary systems

C.1 Testing the stability of sequential planet pairs

The dynamical spacing $\Delta$ describes the separation of two planets in units of their mutual Hill radius. The mutual Hill radius of two planets is given by

[TABLE]

where $m_{i}$ and $m_{j}$ are the planet masses for the inner and outer planets respectively. The dynamical spacing $\Delta$ is the semi-major axis spacing of the two planets, in units of the mutual Hill radius.

[TABLE]

where $a_{i}$ and $a_{j}$ are the semi-major axes of the inner and outer planets respectively. Analytic stability solutions exist for a system which contains exactly two planets, $\Delta_{ij}\gtrsim 3.46$ (Gladman, 1993; Chambers et al., 1996), although it is not possible to ensure this requirement for our simulated systems555 Although our simulated systems may produce exactly two planets within our parameter space, we cannot rule out the possibility of additional planets outside of our parameter space, which would invalidate the analytic solution.. For systems with $n_{i}\geq 3$, we use an empirical stability criteria for two adjacent planet pairs (three sequential planets). A set of three sequential planets with indicies $i$, $j$, and $k$ is deemed unstable when

[TABLE]

where $\Delta_{ij}$ and $\Delta_{jk}$ are the dynamical spacing of the inner and outer planet pair from the three sequential planets (Lissauer et al., 2011). If there are only two simulated planets in a system, $\Delta_{ij}<10$ results in the system being labeled unstable.

Should any set of planets fail the above stability criteria, the system is deemed unstable and new planet periods are redrawn for all $n_{i}$ planets as in step 5. New planetary radii are not redrawn, since passing the stability criteria is biased towards sets of planets with small planetary radii, where stability is more easily achieved. Redrawing planetary radii immediately would result in a simulated $R_{p}$ distribution skewed towards small $R_{p}$, relative to the underlying distribution in Section 4. Should the stability criteria fail $10^{3}$ times for the same set of planetary radii, new $R_{p}$ and periods for all $n_{i}$ planets are redrawn as in step 4.

C.2 Orbital plane inclinations relative to the observer

Once stability has been established for a model system, each planet is then assigned an inclination relative to the observer, the $i$ variable commonly seen in transit and radial velocity detections.

For Rayleigh distributed mutual inclinations (flared disk in Figure 1), $i$ is assigned as follows. An inclination $\phi$ is drawn from a Rayleigh distribution with mode $\sigma_{\phi}/\sqrt{2}$, where $\sigma_{\phi}$ is the mode of the Rayleigh distributed mutual inclinations. The factor of $1/\sqrt{2}$ is a conversion factor between the Rayleigh distributed mutual inclinations, and the Rayleigh distributed planet inclinations around the invariable plane. The orbital plane of the planet is then rotated by a random uniform angle $\Omega$, giving

[TABLE]

For a flat disk (Figure 1), the perpendicular height above the invariable plane $Z_{0}$ is drawn from a Gaussian distribution with a mean of 0 and standard deviation $\sigma_{Z}$, in units of stellar radii. For a flat disk, unlike a flared disk, the assigned inclination $i$ is dependent on the semi-major axis of the planet. Again, the orbital plane of the planet is rotated by a random uniform angle $\Omega$, to account for a random viewing angle.

[TABLE]

resulting in a tendency for larger inclinations for close-in planets and vice-versa (right planel of Figure 1).

Appendix D Orbit-normalised transit duration ratio $\bm{\xi}$

For a planet which transits through the centre of its star:

[TABLE]

where $v_{\text{orb}}$ and $T_{\text{dur}}$ represent the orbital velocity (assuming a circular orbit) and the transit duration of the planet respectively. Note that for the Kepler sample, the simplification of a circular orbit is justified since $\xi$ is weakly dependent on eccentricity (Fabrycky et al., 2014). In addition, eccentricity values for the Kepler sample are generally found to be associated with near-circular orbits (e.g. Fabrycky et al. (2014); Hadden & Lithwick (2014); Dawson et al. (2016)), or with mean values around $\sim$0.1 (Moorhead et al., 2011; Hansen & Murray, 2013).

When the transit is not through the centre of the star $\left(2R_{*}=2\sqrt{R_{*}^{2}-b^{2}}\right)$:

[TABLE]

where $b$ is the impact parameter representing the transiting planet. For a pair of planets which transit the same host star:

[TABLE]

where the "in" and "out" subscripts represent the inner and outer planets respectively. From Kepler’s 3rd law:

[TABLE]

Dividing 25 by 26:

[TABLE]

The RHS ratio is particularly useful for planetary transits as it is composed of well-measured variables. Setting the RHS to $\xi$ (Steffen et al., 2010):

[TABLE]

From 27, a coplanar planetary pair will only give $\xi=1$ if the invariable plane (Fig 10) is exactly edge-on to the observer. For inclined invariable planes, a coplanar planetary pair will give $\xi>1$, as $b_{\text{out}}>b_{\text{in}}$. Values of $\xi<1$ are due to $b_{\text{out}}<b_{\text{in}}$, and are not possible in cases of perfect coplanarity.

Appendix E Previous coplanarity studies

Comparison of exoplanet coplanarity studies

Planet sample

Reference $\Delta\phi$ distribution Observables Dispersiona Sample (quarter, multiplicity) Period (days) Radius ($R_{\oplus}$) Stellar sample Dichotomyg

Lissauer et al. (2011) Rayleigh $\bm{\hat{N_{k}}}^{\,b}$ $\sigma_{\phi}\sim 2.0^{\circ}$ Kepler (Q2, 1-6) $3-125$ $1.5-6$ FGK dwarfs 2.8

Tremaine & Dong (2012) Fisher $\bm{\hat{N_{k}}}$ $\sigma_{\phi}^{c}<4.0^{\circ}$ RV & Kepler (Q2, 1-6) $<200$ $<22$ FGK dwarfs

Figueira et al. (2012) Rayleigh $\bm{\hat{N_{k}}}$ $\sigma_{\phi}^{d}\sim 1.4^{\circ}$ HARPS & Kepler (Q2, 1-3) $<50$ $>2$ FGK dwarfs

Fang & Margot (2012) Rayleigh, R of R $\bm{\hat{N_{k}}}$, $\bm{\xi}^{e}$ $\sigma_{\phi}^{c}\sim 1.4^{\circ}$ Kepler (Q6, 1-6) $<200$ $1.5-30$ FGK dwarfs 1

Johansen et al. (2012) uniform $i$ + rotation${}^{\,f}$ $\bm{\hat{N_{k}}}$ $\sigma_{\phi}<3.5^{\circ}$ Kepler (Q6, 1-3) $<240$ $<22$ FGK dwarfs 3

Weissbein & Steinberg (2012) Rayleigh $\bm{\hat{N_{k}}}$ no fit Kepler (Q6, 1-6) $<75^{h}$ $-$ FGK dwarfs

Hansen & Murray (2013) Rayleigh $\bm{\hat{N_{k}}}$

Kepler (Q6, 1-6) $<1.1$ AU

2

Fabrycky et al. (2014) Rayleigh $\bm{\xi}$ $\sigma_{\phi}\sim 1.8^{\circ}$ Kepler (Q6, 1-6) $<130^{h}$ $-$ FGK dwarfs

Ballard & Johnson (2016) Rayleigh $\bm{\hat{N_{k}}}$ $\sigma_{\phi}=2.0^{\circ\,+4.0}_{\;\;\;-2.0}$ Kepler M-dwarfs (Q16, 1-5) $1-200$ $-$ M stars 3

Gaidos et al. (2016) Rayleigh $\bm{\hat{N_{k}}}$ $\sigma_{\phi}\sim 0^{\circ}$ Kepler M-dwarfs (Q16, 1-5) $<180$ $1-4$ M stars

This paper Rayleigh / Flat disk $\bm{\hat{N_{k}}},\bm{\xi}$ $1.6^{+0.6}_{-0.3}$ Kepler (Q16, 1-6) $3-200$ $1-5$ FGK dwarfs

• a

The mode of the Rayleigh distribution of $\phi$ values (Fig. 10) around the invariable plane.

• b

$\bm{\hat{N_{k}}}$ is the multiplicity vector for the numbers of observed k-planet systems, i.e. $\bm{\hat{N_{k}}}=[\hat{N_{1}},\hat{N_{2}},\hat{N_{3}},...]$.

• c

Converted from the mean $\mu$ of the mutual inclination Rayleigh distribution: $\sigma_{\phi}=\sqrt{2/\pi}\>\>\sigma_{i}$.

• d

Converted from Rayleigh distribution relative to the invariable plane: $\sigma_{\phi}=\sqrt{2}\;\sigma_{\Delta\theta}$.

• e

$\bm{\xi}$ is the normalized transit duration ratio (Appendix D) as given in Steffen et al. (2010).

• f

Each planet is given a random uniform inclination between $0^{\circ}-5^{\circ}$. This orbital plane is then rotated uniformally between $0-2\pi$ to give a random longitude of ascending node.

• g

The factor by which the number of simulated 1-planet systems are lower than observed

• h

Converted from a maximum semi-major axis, assuming a Solar mass star