A Flat Inner Disk Model as an Alternative to the Kepler Dichotomy in the Q1 to Q16 Planet Population
Timothy Bovaird Charles H. Lineweaver
 
firstpage–lastpage 2016

Quantum phase transitions of light in a dissipative Dicke-Bose-Hubbard model
    Wu-Ming Liu
    December 30, 2023
============================================================================




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


 
exoplanets, Kepler, inclinations, multiple-planet systems, invariable plane



§ INTRODUCTION

The Kepler Q1-Q16 catalog <cit.> uses 47 months of Kepler data collected from ∼190,000 stars. 
This has resulted in the detection of over 4,000 planet candidates orbiting ∼3,200 stars. 
While the majority of the ∼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.



 §.§ The Kepler Dichotomy



The mutual inclination distribution between planets around Kepler stars has been well studied
(<cit.>, see Appendix <ref>). 
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 ∼1^∘-5^∘. 

In contrast to the agreement for mutual inclinations, 
some studies report a significant underproduction of simulated systems with a single detected planet <cit.>. 
These studies underproduce the number of simulated stars with a single detected transiting planet by a factor of ∼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” <cit.>, 
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 <cit.>, 
instability or suppressed planet formation caused by stellar binaries <cit.>,
varying surface density profiles and disk masses <cit.>,
varying strengths of gas depletion
or spin-orbit misalignment between the star and planet <cit.>.
<cit.> 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 ∼55 per cent of M dwarfs in the Kepler field.





 §.§ Detected Transiting Planets


We define the true planetary system multiplicity vector  as the number of stars which are host to k planets.
For the Kepler mission (and any transit survey), 
the observed vector  will be significantly lower than the true , 
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

    = [N̂_̂1̂,N̂_̂2̂,N̂_̂3̂,..],

where N̂_̂1̂, N̂_̂2̂ and 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≤ k ≤ 6 and =[2608, 413, 141, 52, 18, 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 <ref>). 
In Figure <ref>, 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 √(Δθ^2+Δθ_(y-z)^2).



 §.§ An Alternative Disk Model


In the studies mentioned in Section <ref>, the true mutual inclinations between simulated planets are drawn from a Rayleigh distribution with mode σ_Δϕ. 
The Rayleigh distribution is composed of two Gaussian distributed components, with standard deviations equal to σ_Δϕ.
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 <ref>. 
Rayleigh distributed mutual inclinations represent a `flared disk' model, where a planet's height above the invariable plane[
The mode of the Rayleigh distribution of inclinations relative to an invariable plane σ_ϕ, 
is related to the Rayleigh distribution of mutual inclinations σ_Δϕ, by σ_ϕ≈σ_Δϕ/√(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 <ref>).
<cit.> 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 <ref>.
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 <ref> we define our stellar and planetary samples, based on minimizing false positives and false negatives.
In Section <ref> we estimate the transit and detection completeness across our parameter space.
In Section <ref> we estimate the underlying orbital period and planet radius distributions.
We outline the process of producing model planetary populations in Section <ref>.
In Section <ref>,
we compare the simulated detections in our model systems to the Kepler Q1-Q16 candidates, for both flat and flared disk models.
In Section <ref> we discuss the results from our model planet populations, 
including estimates for the mean number of planets per star within our parameter space.



§ 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 catalog[
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 <ref>.
We apply the following cuts to the input catalog:

    K <  T_eff    <  K
    σ_CDPP_45   <  ppm
    
    R_*    <   R_⊙
    
    T_baseline   >  days
    
    f_duty   >

where T_eff and R_* are the stellar effective temperature and radius respectively. 
The 4.5 hour CDPP (combined differential photometric precision, <cit.>) of the star, 
a measure of the combined instrumental and stellar noise, is given by σ_CDPP_45.
T_baseline is the timespan of observations for each star and f_duty is the fraction of valid observations over T_baseline.
Note that the combination of T_baseline > 1000 days and f_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  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  Kepler Objects of Interest (KOIs) labeled `candidate' by the Q1-Q16 pipeline[
http://exoplanetarchive.ipac.caltech.edu/cgi-bin/TblView/nph-tblView?app=ExoTbls&config=koi].
An additional  KOIs are labeled `not dispositoned'. We update these dispositions using the Kepler Q17 catalog for reference.
This results in  KOIs changing from `not dispositioned' to `candidate'.

These planets form our initial sample of  planet candidates from the Q1-Q16 catalog <cit.>.
To conform with our input stellar sample, we remove planets around host stars outside of our stellar parameter space defined by Equation <ref>.
This reduces our sample of observed planets from  to . 

We set an upper orbital period limit of  days to avoid the increase in false-positives towards the Kepler orbital period of ∼372 days <cit.>, 
and to remain consistent with the Kepler pipeline completeness calculations <cit.>.

The Kepler pipeline is known to have an increasing false-negative rate with decreasing orbital period for orbital periods ≲ 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 <cit.>.
In addition, a small fraction of the fitted planetary radii for planets with orbital periods ≲ 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  days, in order to retain a sample of Kepler stars with ≥ 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 R_⊕ and a lower limit MES of  are chosen since false-positives are dominated by low MES (≲ 8) detections <cit.>. 
An upper planet radius limit of R_⊕ 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:

    days <    P <  days
     R_⊕ <    R_p <  R_⊕
    
    M   ES >

This results in our observed sample of  candidates in  planetary systems, within the parameter space outlined in Equations <ref> and <ref>.
The observed planetary system multiplicity vector  (Equation <ref>) for our parameter space is given by

    = .




§ TRANSIT AND DETECTION EFFICIENCY



When attempting to estimate the underlying multiplicity vector  given the observed  (Equation <ref>), 
there exists a degeneracy between the underlying multiplicity distribution and the underlying mutual inclination distribution.
For example, an observed  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 <ref> and <ref>).
We produce sets of  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  for a set of simulated systems to the  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.



 §.§ 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 <cit.>, shown in Figure <ref>.

We follow the approximation of the pipeline MES by <cit.>, which includes a limb-darkening approximation and accounts for nonzero impact parameters.

    MES = 0.84 δ(c+s√(δ))/σ_cdpp√(n_tr)

where c=1.0874 and s=1.0187 for G dwarfs, 
δ= (R_p/R_*)^2, and n_tr is the number of transits.
Values of σ_cdpp are reported for 14 different transit durations from 1.5 hours to 15 hours for each Kepler star <cit.>.
The σ_CDPP value chosen for Equation <ref> 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_tr=(T_baseline× f_duty) / P, where P is the planet period,
T_baseline is the total observing time for the Q1-Q16 catalog (∼1426 days) 
and f_duty is the duty cycle for the observed star; the fraction of valid observations over the observing baseline.
Note that T_baseline and f_duty are reported for each star, accounting for systematics such as the differences in CCD detectors and pixels.

We define η_detect as the Kepler pipeline completeness, shown in Figure <ref>.
The pipeline completeness as a function of the multiple event statistic is approximately represented by the Γ cumulative distribution function

    η_detect(MES) = 1/c^b Γ(b)∫_0^MES-β x^b-1e^-x/cdx

where Γ is the Gamma function.
For our sample of FGK dwarfs, b=4.35, c=1.05 and β=4.093 <cit.>.

To calculate the total probability of transit detection η(P,R_p), 
we must also take into account the geometric transit probability η_transit of a planet,

    η_transit = R_*/a,

where a is the semi-major axis of the planet.
The product of these two equations gives the total transit detection probability η(P,R_p),
the probability of the planet transiting η_transit,
and the probability of the transiting planet being detected by the Kepler Q1-Q16 pipeline η_detect,

    η(P,R_p) = η_transit×η_detect


Given an input star from our stellar sample (Equation <ref>) and using Equations <ref> to <ref>, 
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 <ref>),
we estimate the mean total transit and detection completeness ⟨η(P,R_p)⟩,
by taking the mean value of η(P,R_p) at each grid point over all stars in our stellar sample.
This is shown in Figure <ref>,
where ⟨η(P,R_p)⟩ ranges from ∼0 to a maximum of ∼0.1.
Transit and detection probabilities >5 per cent only exist for planets with orbital periods ≲8 days.
It can be seen that the pipeline detection probability becomes important for planetary radii less than 2.5 Earth radii.





§ 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 <cit.>, 
with a logarithmic plateau at ≲ 2.5R_⊕.
This logarithmic plateau is also seen when the pipeline efficiency is probed using transit injection and recovery experiments <cit.>.

For orbital periods between 50 and 300 days, a single power law is sufficient to describe the orbital period distribution <cit.>. 
Our parameter space includes planets with orbital periods less than 50 days, 
where the transit and detection completeness is more dynamic,
particularly for periods ≲15 days (Figure <ref>).
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, β_1, β_2, P_brk, α_1, α_2, R_brk, 
where F_0 is the number of planets per star within our parameter space, 
and R_brk and P_brk are the transition points between the two power laws for the planet radius and orbital period respectively, 

    df/dP dR_p = C F_0  g(P,R_p) 
     =
       
    C F_0 P^β_1 R_P^α_1   if P < P_brk and R_p < R_brk
     
    C F_0 P^β_1 R_P^α_2 R_brk^α_1-α_2   if P < P_brk and R_p ≥ R_brk
    
    C F_0 P^β_2 P_brk^β_1-β_2 R_P^α_1   if P ≥ P_brk and R_p < R_brk
    
    C F_0 P^β_2 P_brk^β_1-β_2 R_P^α_2 R_brk^α_1-α_2   if P ≥ P_brk and R_p ≥ R_brk

where the power law exponents α_1, α_2 and β_1, β_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 ⟨η(P,R_p)⟩ 
(Figure <ref>).
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 <ref>) is shown in Appendix <ref>, 
with maximum likelihood parameters of

Maximum likelihood PLDF parameters 

F_0     β_1     β_2     P_brk     α_1     α_2     R_brk 

0.852     1.007     -0.932     15.332     -1.168     -4.906     2.740 





indicating breaks in the power law distributions at ∼15 days and ∼2.7 R_⊕ 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≲2.7 R_⊕.
These results appear to be consistent with <cit.>, 
where breaks in the logarithmic orbital period and planet radius rates are indicated at ∼20 days and ∼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 <ref> and <ref> 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 ⟨η(P,R_p)⟩ for the stars in our sample.

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



 §.§ 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, . 
In this paper we trial two different  distributions. 

The first trial distribution is a modified Poisson distribution N_pl,Poi <cit.>. 
Each star is assigned a random number of planets, drawn from a Poisson distribution with mean . 
Stars drawn with zero planets are redrawn from the same Poisson distribution, 
until all  model planetary systems are populated with planets,
resulting in a final mean ≥. 

The second trial distribution is a modified Exponential distribution N_pl,exp <cit.> and is produced in the same way as N_pl,Poi, 
except that stars are assigned a number of planets drawn from an exponential distribution with mean .
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 <cit.>.
We include this trial distribution for comparative purposes.



§ 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 , 
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  ranges from 0.5 to 3.5 in steps of .
For the flared disk model, the mode of the Rayleigh distributed mutual inclinations ranges from 0 to 5 degrees in steps of  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  R_*.

This results in a total of  grid points for each set of model assumptions, 
with  model planetary systems generated at each grid point.
Each model planetary system is produced as follows:



(Space*=0.1cm)
    A random star is chosen from our sample of  Kepler stars outlined in Section <ref>, 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  (Figure <ref>),  is chosen from a random point on a sphere: 0 ≤cos[] ≤ 1.

    The number of planets in the system  is drawn randomly, according to the assumed  distribution from Section <ref>, with a mean value based on the current grid point.

    The radii of the  planets are drawn from the underlying distribution in Section <ref> (R_⊕ < R_p <  R_⊕), and converted to their corresponding masses[M_p≈(R_p/a)^b, where a∼1.11 and b∼2.41.].

    The periods of the  planets are drawn randomly from the underlying distribution in Section <ref> (days< P<  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 =2. 
If any set of planets in the system is deemed unstable, the system is labeled unstable and new planet periods for all  planets are redrawn as in step 5.
See Section <ref> 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  (step 2) and the assumed planetary disk model (Section <ref>),
and the parameter value at the current grid point.
See Appendix <ref> for a complete description of how planet inclinations are assigned for flat and flared disk models.


The above steps generate the  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.



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

    b =a/R_*cos i ≤ 1

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 <ref>).
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 η_detect (Equation <ref>).
For each simulated transiting planet, a uniform random number Y_m is drawn between 0 and 1.
A simulated planet is labeled as detected if it transits, 
and if its pipeline detection efficiency η_detect > Y_m.
All simulated planets which meet this criteria are added to the detected planet sample for the grid point, X_ij, 
where i and j represent the current grid point.



§ COMPARING SIMULATED AND OBSERVED PLANET DETECTIONS



The simulations outlined in Section <ref> were performed across a grid for the 4 sets of model assumptions. 
For each grid point, the simulated planet detections X_ij are used to generate two distributions,
the system multiplicity vector  (Equation <ref>),
and the distribution of orbit normalized transit duration ratios  <cit.>.
Unlike the  distribution, 
the  distribution only consists of model systems with two or more detected transiting planets.
For a pair of planets orbiting the same star,

    ξ=T_dur,in/P_in^1/3/T_out,in/P_out^1/3,

where T_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, ξ is calculated, 
giving (+1)/2 values of ξ for a star with  planets.
For each grid point we generate the ensemble _ij distribution
by calculating the ξ value for each unique pair of simulated transit detections orbiting the same star, 
across all  model systems. 
For a deeper discussion of ξ, see Appendix <ref>.

The  and  distributions are compared to the  and  distributions of the Kepler Q1-Q16 candidates, 
and are used to assess the goodness of fit at each grid point.
We perform a χ^2 goodness of fit test (Equation <ref>) 
comparing the simulated  to the observed  for our parameter space (Equation  <ref>).
We scale  such that ∑=∑.
To compensate for the poor quality of the χ^2 test with low cell counts, values less than 5 are merged into their adjacent cells.

    χ^2 = ∑^n_k=1( - )^2/.


Similarly, we perform a two-sample Kolmogorov-Smirnov (KS) test between the simulated _ij distribution at each grid point, and the  distribution of the observed Q1-Q16 Kepler candidates within our parameter space. 



 §.§ Flared Disk and Poisson distributed planets per star


In the top panel of Figure <ref>,  is compared to  at each grid point,
under the assumption of a flared planetary disk and a Poisson distributed number of planets per star.
The χ^2 values are represented by the 1σ, 2σ and 3σ values relative to the best fit.
As expected, no good fit is found to the  distribution, 
as is the case in the majority of previous studies
<cit.>.

The bottom panel displays the resulting p-values from the KS test between the _ij and  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 <ref>. 
The mean number of planet per star cannot be determined from the  distribution alone.

It is clear from Figure <ref> 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 ≲1^∘ (top panel).
However, modes ≲ 1.5^∘ are ruled out by comparing orbit normalized transit duration ratios (bottom panel).







 §.§ 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 <ref>).
Unlike for the flared disk model, the best fit  is a good match to ,
giving a χ^2/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_*≲ Z_0 ≲ 2 R_*. 
There is significant overlap between the two tests within this region. 





 §.§ Combining independent tests


The results from comparing the  and  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,

    χ^2_combined≈ -2 ∑_m=1^M ln p_m

where M is the number of tests combined and p_m is the p-value of the mth 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_combined, across the parameter space for each set of model assumptions.
P_combined is derived from the likelihood P_ij∝exp (-χ_combined,ij^2/2) and the requirement ∑ P_ij = 1.

Figure <ref> displays the probability grids for each set of model assumptions, along with the best-fit point and the 1σ and 2σ probability contours.
Panel a) of Figure <ref> combines the tests of Figure <ref>,
and Panel b) combines the tests of Figure <ref>.
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. 



§ RESULTS AND DISCUSSION





 §.§ 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 =2.0^+0.3_-0.2,
over our parameter space, [3 days <P<200 days] and [1 R_⊕<R_p<5 R_⊕].
For a similar orbital period and planet radius parameter space, <cit.> find =2.2±0.3 for Kepler M dwarfs.
<cit.> report ∼1.5 for R_p>1.5 R_⊕,
where the reduction in  likely comes from the exclusion of planets with radii between 1.0 R_⊕ < R_p < 1.5 R_⊕.

Similarly, we find the mode of the Rayleigh distributed mutual inclinations is given by σ_i = 2.3^+0.9_-0.4 degrees, 
consistent with the bulk of previous results with σ_i∼2^∘
<cit.>.

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 <cit.>,
where a flared disk model reproduced the  distribution without the need for a Kepler dichotomy.
The discrepancy likely comes from the unique N_pl distribution chosen by <cit.>, a “bounded uniform” distribution.
The bounded uniform distribution is produced by first choosing a maximum number  of planets n_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_i,max.

It has previously been shown that the Kepler  distribution can be matched without the need for a Kepler dichotomy,
when N_pl is drawn from an exponential distribution <cit.>.
We find that the Kepler sample is consistent with σ_i = 2.4^+0.9_-0.5 degrees and = 1.6^+0.3_-0.2
drawn from an exponential distribution.
Here we disagree with <cit.>, who preferred near-coplanar mutual inclinations.
While we also achieve good fits for near-coplanar orbital planes, 
comparing  distributions strongly rules out mutual inclinations ≲ 1.4^∘.
This illustrates the importance of modeling both  and  distributions,
where <cit.> only modeled the  distribution.



 §.§ Flat Disk Model


For a flat planetary disk model
<cit.>, 
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 = 2.4^+0.6_-0.4.
Notably, the mean number of planets per star  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 <ref>).
For the inner Solar System planets, Z_max≈8 R_*, giving Z_0≈5 R_*, where Z_0≈2  Z_max/π. 
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.



§ 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_⊕<R_p<5 R_⊕].
We find that both distributions are well described by broken power laws, with breaks occurring at ∼15 days and ∼2.7 R_⊕.
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_pl drawn from Poisson and exponential distributions.

We confirm that a flared planetary disk model with N_pl drawn from a Poisson distribution, 
is not consistent with the Kepler detections. 
We also confirm that Kepler detections are well matched when N_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_pl distribution.

We find that the mean number of planets per star  is largely model independent,
∼2.0 when N_pl is drawn from a Poisson distribution,
and ∼1.6 when N_pl is drawn from an exponential distribution,
for [3 days <P<200 days] and [1 R_⊕<R_p<5 R_⊕].
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 ∼2.2^∘.
For a flat planetary disk model, the Gaussian disk thickness Z_0∼1.5 R_*,
much lower than the ∼5 R_* of the inner Solar System.



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

<cit.> 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, Δ (Equation <ref>), planet multiplicity and ξ 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 ξ and Δ.














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



mn2e 







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






§ PLANET DISTRIBUTION FUNCTION



Our planet distribution function (PLDF) has 7 free parameters, F_0, β_1, β_2, P_brk, α_1, α_2, R_brk.

    df/dP dR_p = C F_0  g(P,R_p) 
     =
       
    C F_0 P^β_1 R_P^α_1   if P < P_brk and R_p < R_brk
     
    C F_0 P^β_1 R_P^α_2 R_brk^α_1-α_2   if P < P_brk and R_p ≥ R_brk
    
    C F_0 P^β_2 P_brk^β_1-β_2 R_P^α_1   if P ≥ P_brk and R_p < R_brk
    
    C F_0 P^β_2 P_brk^β_1-β_2 R_P^α_2 R_brk^α_1-α_2   if P ≥ P_brk and R_p ≥ R_brk

where F_0 is the number of planets per star within our parameter space, R_brk and P_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

    ∫_R_min^R_max∫_P_min^P_max C  g(P, R_p)  dP  dR_p = 1

where the integration limits R_min, R_max, P_min and P_max are given in Equation <ref>.

We follow <cit.> and <cit.> by implementing a Poisson likelihood for our PLDF. By maximizing this likelihood we can obtain best-fit parameters for our model. 

    ln(L) ∝[ ∑^N_pl_i=1ln(C  F_0  g(P,R_P) )]- N_exp

where N_exp is the expected number of planet detections for the set of model parameters, and is given by

    N_exp = C  F_0 ∫_R_min^R_max∫_P_min^P_max[∑_j=1^N_*η_j(P, R_p)] g(P,R_P)  dP dR_p

where η_j(P, R_p)=η_detect×η_transit is the combined transit and pipeline detection efficiency of the jth star for the specified period and radius. 
The pipeline detection efficiency η_detect is given by Equation <ref> and the transit probability η_transit=R_*/a, where a is the semi-major axis.

We calculate [∑_j=1^N_*η_j(P, R_p)] for a grid in orbital period and planet radii, in bins of 1.5 days and 0.05 R_⊕ respectively.
For each grid point, we sum over all stars in our sample. 
The mean combined transit and pipeline detection efficiency ⟨η(P,R_p)⟩ can then be found be dividing this term by the number of stars in our sample, N_*. 





§ SIMULATED PLANETARY SYSTEMS



 §.§ Testing the stability of sequential planet pairs


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

    R_H,ij=(m_i+m_j/3M_*)^1/3a_i+a_j/2

where m_i and m_j are the planet masses for the inner and outer planets respectively. 
The dynamical spacing Δ is the semi-major axis spacing of the two planets, in units of the mutual Hill radius.

    Δ_ij=a_j-a_i/R_H,ij

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, Δ_ij≳3.46 <cit.>, 
although it is not possible to ensure this requirement for our simulated systems[
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 ≥ 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

    Δ_ij+Δ_jk < 18

where Δ_ij and Δ_jk are the dynamical spacing of the inner and outer planet pair from the three sequential planets <cit.>.
If there are only two simulated planets in a system, Δ_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  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 <ref>. 
Should the stability criteria fail 10^3 times for the same set of planetary radii, new R_p and periods for all  planets are redrawn as in step 4.



 §.§ 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 <ref>), i is assigned as follows. 
An inclination ϕ is drawn from a Rayleigh distribution with mode σ_ϕ /√(2), where σ_ϕ is the mode of the Rayleigh distributed mutual inclinations. 
The factor of 1/√(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 Ω, giving

    i_flare=+ϕ cos(Ω).


For a flat disk (Figure <ref>), the perpendicular height above the invariable plane Z_0 is drawn from a Gaussian distribution with a mean of 0 and standard deviation σ_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 Ω, to account for a random viewing angle.

    i_flat=  + arcsin(Z_0/a) cos(Ω)

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



§ ORBIT-NORMALISED TRANSIT DURATION RATIO 


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

    2R_*≈ v_orbT_dur

where v_orb and T_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 ξ is weakly dependent on eccentricity <cit.>.
In addition, eccentricity values for the Kepler sample are generally found to be associated with near-circular orbits (e.g. <cit.>), or with mean values around ∼0.1 <cit.>.

When the transit is not through the centre of the star (2R_* = 2√(R_*^2-b^2)):

    2√(R_*^2-b^2) = v_orbT_dur

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

    2√(R_*^2-b_in^2)=T_dur,in v_orb,in
       2√(R_*^2-b_out^2)=T_dur,out v_orb,out

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

    v_orb,in∝a_in/P_in∝ P_in^-1/3

Dividing <ref> by <ref>:

    √(R_*^2-b_in^2)/√(R_*^2-b_out^2) = T_dur,in/P_in^1/3/T_dur,out/P_out^1/3

The RHS ratio is particularly useful for planetary transits as it is composed of well-measured variables. Setting the RHS to ξ <cit.>: 

    ξ=T_dur,in/P_in^1/3/T_out,in/P_out^1/3

From <ref>, a coplanar planetary pair will only give ξ=1 if the invariable plane (Fig <ref>) is exactly edge-on to the observer. For inclined invariable planes, a coplanar planetary pair will give ξ>1, as b_out > b_in. Values of ξ<1 are due to b_out<b_in, and are not possible in cases of perfect coplanarity. 




§ PREVIOUS COPLANARITY STUDIES



	Comparison of exoplanet coplanarity studies
	
	
	
		                    2|c|Planet sample 

		(lr)6-7
		Reference     Δϕ distribution     Observables     Dispersion^a     Sample (quarter, multiplicity)     Period (days)     Radius (R_⊕)     Stellar sample     Dichotomy^g

		<cit.>     Rayleigh     ^ b     σ_ϕ∼2.0^∘     Kepler (Q2, 1-6)     3-125     1.5-6     FGK dwarfs     2.8

		<cit.>     Fisher          σ_ϕ^c <4.0^∘     RV & Kepler (Q2, 1-6)     <200     < 22     FGK dwarfs     - 

		<cit.>     Rayleigh          σ_ϕ^d ∼ 1.4^∘     HARPS & Kepler (Q2, 1-3)     <50     >2     FGK dwarfs     -  

		<cit.>     Rayleigh, R of R     , ^e     σ_ϕ^c∼1.4^∘     Kepler (Q6, 1-6)     <200     1.5-30     FGK dwarfs     1 

		<cit.>     uniform i + rotation^ f          σ_ϕ<3.5^∘     Kepler (Q6, 1-3)      <240     <22     FGK dwarfs     3

		<cit.>     Rayleigh          no fit     Kepler (Q6, 1-6)      <75^h     -     FGK dwarfs     -

		<cit.>     Rayleigh          -     Kepler (Q6, 1-6)     <1.1 AU     -     -     2

		<cit.>     Rayleigh          σ_ϕ∼1.8^∘     Kepler (Q6, 1-6)     <130^h     -     FGK dwarfs     - 

		<cit.>     Rayleigh          σ_ϕ = 2.0^∘ +4.0_   -2.0     Kepler M-dwarfs (Q16, 1-5)     1-200     -     M stars     3

		<cit.>     Rayleigh          σ_ϕ∼ 0^∘     Kepler M-dwarfs (Q16, 1-5)     <180     1-4     M stars     - 

		This paper     Rayleigh / Flat disk     ,     1.6^+0.6_-0.3     Kepler (Q16, 1-6)     3-200     1-5     FGK dwarfs     - 

		
	[]
		
  a The mode of the Rayleigh distribution of ϕ values (Fig. <ref>) around the invariable plane.
		
  b  is the multiplicity vector for the numbers of observed k-planet systems, 
		i.e. = [N̂_̂1̂, N̂_̂2̂, N̂_̂3̂,...].
		
  c Converted from the mean μ of the mutual inclination Rayleigh distribution: σ_ϕ=√(2/π)  σ_i.
		
  d Converted from Rayleigh distribution relative to the invariable plane:  σ_ϕ=√(2) σ_Δθ.
		
  e  is the normalized transit duration ratio (Appendix <ref>) as given in <cit.>.
		
  f Each planet is given a random uniform inclination between 0^∘-5^∘. This orbital plane is then rotated uniformally between 0-2π 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
         




