Architectures of Exoplanetary Systems. I: A Clustered Forward Model for Exoplanetary Systems around Kepler's FGK Stars
MMNRAS , 1–31 (2019) Preprint 18 November 2019 Compiled using MNRAS L A TEX style file v3.0
Architectures of Exoplanetary Systems. I: A Clustered ForwardModel for Exoplanetary Systems around
Kepler ’s FGK Stars
Matthias Y. He , , , (cid:63) , Eric B. Ford , , , , and Darin Ragozzine Department of Astronomy and Astrophysics, 525 Davey Laboratory, The Pennsylvania State University, University Park, PA 16802, USA Center for Exoplanets and Habitable Worlds, 525 Davey Laboratory, The Pennsylvania State University, University Park, PA 16802, USA Center for Astrostatistics, 525 Davey Laboratory, The Pennsylvania State University, University Park, PA 16802, USA Institute for CyberScience, The Pennsylvania State University, University Park, PA 16802, USA Department of Physics and Astronomy, N283 ESC, Brigham Young University, Provo, UT 84602, USA
Accepted ?. Received ?; in original form ?
ABSTRACT
Observations of exoplanetary systems provide clues about the intrinsic distribution of planetarysystems, their architectures, and how they formed. We develop a forward modelling frameworkfor generating populations of planetary systems and “observed” catalogues by simulating the
Kepler detection pipeline (
SysSim ). We compare our simulated catalogues to the
Kepler
DR25catalogue of planet candidates, updated to include revised stellar radii from Gaia DR2. Weconstrain our models based on the observed 1D marginal distributions of orbital periods,period ratios, transit depths, transit depth ratios, transit durations, transit duration ratios, andtransit multiplicities. Models assuming planets with independent periods and sizes do notadequately account for the properties of the multiplanet systems. Instead, a clustered pointprocess model for exoplanet periods and sizes provides a significantly better description of the
Kepler population, particularly the observed multiplicity and period ratio distributions. Wefind that 0 . + . − . of FGK stars have at least one planet larger than 0 . R ⊕ between 3 and 300 d.Most of these planetary systems ( ∼ Kepler dichotomy is evidence for a population of highlyinclined planetary systems and is unlikely to be solely due to a population of intrinsically singleplanet systems. We provide a large ensemble of simulated physical and observed catalogues ofplanetary systems from our models, as well as publicly available code for generating similarcatalogues given user-defined parameters.
Key words: methods: statistical – planetary systems – planets and satellites: detection,fundamental parameters, terrestrial planets – stars: statistics
Within the past decade, NASA’s
Kepler mission (Borucki et al. 2010,2011a,b; Batalha et al. 2013) has discovered thousands of exoplanetsand hundreds of multiplanet systems around Sun-like stars. Withinthese systems, there is an abundance of short period planets (i.e. withorbits much smaller than that of Earth) and tightly packed multipleplanets (Latham et al. 2011; Lissauer et al. 2011b, 2014; Roweet al. 2014). The wealth of transit detections generated by
Kepler and the relative homogeneity of the sample allows for the study ofexoplanet systems as a whole, enabling statistical exploration of thepopulation of exoplanets detectable by
Kepler .Multitransiting planetary systems are especially valuable be-cause they serve as crucial tests for our models of planetary for-mation, their resulting architectures, and their subsequent evolution (cid:63)
Contact e-mail: [email protected] and stability (Ragozzine & Holman 2010; Ford 2014; Winn &Fabrycky 2015). The abundance of systems with many transitingplanets indicates that systems with small mutual inclinations arecommon, as small mutual inclination angles between planets inthe same system are required to explain how a single observer cansee so many planets in a transiting configuration. These observa-tions contribute to our theories of planet formation, supporting thepicture that planets likely formed in relatively flat, gaseous discs.Previous studies have explored the degree of coplanarity requiredto explain the multitude of multitransiting systems, focusing almostexclusively on the observed multiplicity and the transit duration ra-tio distributions (Latham et al. 2011; Lissauer et al. 2011b; Fang &Margot 2012; Figueira et al. 2012; Johansen et al. 2012; Tremaine& Dong 2012; Weissbein, Steinberg, & Sari 2012; Fabrycky et al.2014). However, one emergent puzzle from the studies of observedtransiting multiplicities is the apparent excess of single transitingsystems, which has led some authors to speculate the existence of a © 2019 The Authors a r X i v : . [ a s t r o - ph . E P ] N ov He, Ford, and Ragozzine second population of intrinsic singles or highly inclined multiplanetsystems, a so-called “
Kepler dichotomy” (Lissauer et al. 2011b; Jo-hansen et al. 2012; Hansen & Murray 2013; Ballard & Johnson2016).It is unclear how the planet multiplicities and their orbital in-clinations mesh with other observed properties of the transitingpopulation, such as the distributions of their orbital periods and pe-riod ratios. The period ratios of adjacent planet pairs in multiplanetsystems, which is a direct measure of their physical separations andthus stability, have been studied (Fabrycky et al. 2014; Steffen &Hwang 2015). The period ratios range from close to unity to overseveral dozen, with 1.172 being the smallest observed period ratio inthe KOI-1665 system (Lissauer et al. 2011a). Smaller period ratiosof 1.038 (in KOI-284) and 1.065 (in KOI-2248) are now interpretedas due to transiting planets with very similar periods around dif-ferent stars in a binary system (Lissauer et al. 2011a, 2014). Thesestudies also find relative excesses of (apparently) adjacent planetpairs with period ratios near (slightly larger than) first-order mean-motion resonances (MMRs), namely near the 3:2 and 2:1 ratios(Fabrycky et al. 2014; Steffen & Hwang 2015). Inferring the truerate of near resonant systems is difficult, due to the interaction ofthe geometric joint transit probability, detection efficiency, and theunknown distribution of mutual inclinations and eccentricities.A few studies have also attempted to understand the size dis-tribution of planets in multiplanet systems, by probing their relativeradii using transit depth ratios in order to minimize the effects of ouruncertainties in the stellar radii. For example, Ciardi et al. (2013)found that for planet pairs with planets larger than ∼ R ⊕ , the outerplanet tends to be larger than the inner one, although they did notobserve this trend for planets (cid:46) R ⊕ . In addition, planet sizes ap-pear to be highly correlated, as evidenced by the peaked nature ofthe adjacent radii ratio distribution (Weiss et al. 2018a), and thisclustering extends to planet masses (Millholland, Wang, & Laugh-lin 2017). While the mechanisms of photoevaporation (Owen & Wu2013; Fulton et al. 2017; Owen & Wu 2017; Van Eylen et al. 2017;Carrera et al. 2018) or core-powered mass-loss from formation(Ginzburg, Schlichting, & Sari 2016, 2018; Gupta & Schlichting2018) have been proposed to explain these features in the distribu-tions of radii and radii ratios, complex observational biases limit ourability to distinguish models or understand the relative contributionof physical and observational effects. To further illustrate this point,most recently Zhu (2019) has challenged the statistical significanceof the correlations reported by Ciardi et al. (2013) and Weiss etal. (2018a). Zhu (2019) found that testing the significance of cor-relations by bootstrap re-sampling with cuts on the signal-to-noiseratios (as opposed to planet size) resulted in smaller and less sig-nificant correlations for sizes of planets in one system, uniformityof spacing, and preference for the outer planet to be larger than theinner planet. They suggest that a better way to study such correla-tions would be via forward modelling the detection and selectionprocesses using a model for the intrinsic distribution of planetaryproperties. This study takes that approach and provides a fresh per-spective on the significance of correlations in orbital period, planetsize, and uniformity of spacing.There is an overwhelming wealth of information in the ob-served Kepler population of exoplanets, especially encoded in themultiplicity, period ratio, and transit depth and duration ratio dis-tributions. Exploratory analyses are difficult to interpret due tocomplex geometric and sensitivity biases in the
Kepler detectionpipeline. Earlier works such as those discussed above attemptedto study these features but had fewer detections to work with, amore limited understanding of the
Kepler pipeline’s detection effi- ciency, greater uncertainties in the stellar properties, and intractablebiases in the planet catalogue due to human-influenced vetting ofplanet candidates. In this work, we are able to address these con-cerns by making use of several recent analyses and updates to the
Kepler stellar and planetary properties. We use the final
Kepler
DR25 catalogue (Thompson et al. 2018), which was generated us-ing a completely automated procedure for the vetting of planetcandidates using the
Kepler Robovetter . We rely and build on theExoplanets Systems Simulator (“SysSim”), which makes use ofseveral additional DR25 data products describing
Kepler ’s detec-tion pipeline (Burke & Catanzarite 2017a,b,c; Christiansen 2017;Coughlin 2017), and is described in Hsu et al. (2018, 2019). Finally,we adopt the improved stellar properties from the European SpaceAgency’s Gaia mission (Gaia Collaboration et al. 2018), and planetcandidates around a clean sample of FGK dwarfs as defined by Hsuet al. (2019).In this paper, we outline a new framework for simulating cata-logues of observed transiting planets in §2, which we use to exploreseveral statistical, yet physically motivated, models for the intrin-sic distribution of planetary systems. In particular, we explore amarked, clustered point process for generating planetary systems.For each planetary system, we attempt to generate the planets byfirst drawing a period and radius scale for each cluster and thendrawing orbital periods and planet sizes for each planet in a clusterconditioned on the cluster’s period and radius scales. These modelsare described in §2.1-§2.3. We define a set of summary statistics,a distance function to compare models to data, and a subset of the
Kepler
DR25 catalogue we use to constrain our models, in §2.4.We describe the optimization procedure used to explore the mul-tidimensional parameter space of each model in §2.5. In §2.6, weprovide a brief review of the statistical machinery of Gaussian pro-cesses, and describe how we adapt it to serve as an “emulator” forour models and compute the credible regions for each model pa-rameter using Approximate Bayesian Computation. The results forthe credible regions of the parameters for each model are presentedand discussed in §3. The results of each of our models, for boththe observed and physical (intrinsic) distributions of planetary sys-tems, and a direct comparison between them are described in §4.We compute some estimates for the fraction of stars with planetsand the number of planets per star or planetary system in §4.3. Adiscussion of the
Kepler dichotomy is presented in §4.4. We discusshow our models can be improved in §4.6, and summarize our mainresults and encourage the use of our code and model catalogues in§5.
The models described in this paper are developed in the frameworkof the Exoplanets Systems Simulator, which we refer to as “SysSim”.The SysSim codebase is written in the Julia language (Bezanson etal. 2014) and can be installed as the ExoplanetsSysSim.jl package(Ford et al. 2018b). Our models can be accessed at https://github.com/ExoJulia/SysSimExClusters .Our general framework for performing an (approximate)Bayesian analysis can be summarized with the following steps:
Step 0: Define a statistical description for the intrinsic distri-bution of exoplanetary systems
We present three separate statistical models in this paper (§2.1).Each model is described by a set of physically motivated modelparameters.
MNRAS000
MNRAS000 , 1–31 (2019) lustered Planetary Systems Step 1: Generate an underlying population of exoplanetarysystems ( physical catalogue ) from a given model
Each model involves randomly drawing stars from the
Kepler stellarcatalogue and populating them with planets. Planet radii and periodsare drawn. Then, masses are assigned from the radii and a rejection-sampling algorithm is applied to only keep planetary systems thatare physical.
Step 2: Generate an observed population of exoplanetary sys-tems ( observed catalogue ) from the physical catalogue
The systems are oriented by assigning orbits to each planet and theSysSim forward model simulating the
Kepler detection pipeline isapplied to generate a catalogue of observed, transiting planets alongwith their “measured” properties.
Step 3: Compare the simulated observed catalogue with the
Kepler data
A set of summary statistics are computed for both the simulatedobserved and
Kepler data, and a distance function is computed tocompare the two catalogues.
Step 4: Optimize the distance function to find the best-fittingmodel parameters
Steps 0–3 are repeated by passing the distance function into anoptimization algorithm in order to explore the parameter space andattempt to find the best-fitting model parameters.
Step 5: Explore the posterior distribution of model parame-ters using a Gaussian Process (GP) emulator
For each model, a GP model is trained on a set of points (model pa-rameters with computed distances) in order to form an emulator forthe full forward model. The emulator is used to quickly characterizethe parameter space.
Step 6: Compute credible intervals for model parameters andsimulated catalogues using Approximate Bayesian Computing
Model parameters are drawn from prior distributions, and used todraw an emulated distance from the GP emulator. If the emulateddistance is less than a specified distance threshold, we accept the setof model parameters. The credible regions are reported based on theaccepted sample of model parameters, and we present and analyzethe simulated physical and observed catalogues that resulted in thedistance from the full forward model being less than the distancethreshold.We describe each of these steps in full detail in the followingsubsections.
We describe three separate models, starting with a baseline non-clustered model, extending to a clustered periods model, and thena clustered periods and sizes model. First, we give a brief overviewof all three models. Then we provide details about the method fordrawing each of the properties under the different models. Resultsfrom the non-clustered model shown in §4, will provide empiricalsupport for motivating use of the clustered models.In the non-clustered model, we first draw a number of planetsfor each star and then draw orbital periods and planet sizes inde-pendently for each planet in a system, using a simple power lawfor period and a broken power law for radius. In both of the clus-tered models, we first draw a number of “clusters” of planets foreach star. For each cluster, we draw a number of planets and a pe-riod scale. In the clustered periods and sizes model, we also drawa radius scale for each cluster. Then, the periods and sizes of theplanets are drawn from distributions centred on the period scale andthe radius scale, respectively (i.e. conditioned on the properties of
Cluster 2Cluster 1
Orbital period
Star Cluster 2Cluster 1
Clustered Periods and Sizes modelClustered Periods modelNon-clustered model
Figure 1.
Cartoon illustration of our three models, drawn along an arbitraryorbital period axis. In the non-clustered model, the number of planets isdrawn from a Poisson distribution, and the orbital periods and planet radii aredrawn independently from a power law and a broken power law, respectively.In the clustered models, the number of clusters is drawn from a Poissondistribution (two clusters are shown here for illustrative purposes). Similarly,the number of planets per cluster is drawn from a zero-truncated Poissondistribution. For the clustered periods model, the orbital periods of theplanets are clustered together, but their sizes are still drawn independently.For the clustered periods and sizes model, the radii of planets in the samecluster are also similar. The models do not enforce any correlation betweenperiods and planet radii. This paper provides a forward model for drawingplanetary systems from each of the three models and simulating
Kepler observations of such planetary systems. In §4, we show that the two clusteredmodels provide a significantly better match to the
Kepler observations thanthe non-clustered model. their parent cluster). Thus, the properties of planets are explicitlycorrelated with those of other planets from the same cluster. Thisleads to closely spaced planets often having strong correlations, butmore widely spaced planets having weaker correlations. We show acartoon illustration of our three models in Figure 1.
We will constrain our models based on the observed multiplicitydistribution (as described in §2.4.1) and use the results to addressboth the overall rate of planets per star and the fraction of stars withat least one planet in this study. Therefore, it is important to considerthe process for assigning the number of planets to each star.For the non-clustered model, we draw the number of planets ineach system from a Poisson distribution, N p ∼ Poisson ( λ p ) . In theclustered models we first draw a number of clusters from a Poissondistribution, N c ∼ Poisson ( λ c ) . Then, we draw a number of planetsfor each cluster from a zero-truncated Poisson (ZTP) distribution, N p ∼ ZTP ( λ p ) , where the probability mass function is given by: g ( k | λ ) = λ k ( e λ − ) k ! . (1) MNRAS , 1–31 (2019)
He, Ford, and Ragozzine
The choice of the ZTP is simply to avoid drawing empty clusterswith no planets. The way in which zero-planet systems are drawn is an im-portant consideration for our models. In our non-clustered model,assigning a star no planets occurs when N p = N c = N c (cid:54) N c , max and N p (cid:54) N p , max ,respectively. We set N c , max =
10 and N p , max =
15 in our clusteredmodels, and N p , max =
20 in our non-clustered model, unless statedotherwise in order to avoid generating systems with extremely largenumbers of planets, as these are computationally expensive.In addition to truncation effects, the true mean rates of planetsper system (non-clustered model), clusters per system and planetsper cluster (clustered models) are somewhat less than the valuessuggested by the parameters λ p and λ c due to our stability criteriaand rejection sampling algorithm (described in §2.1.7 and §2.2).Thus, while the values of λ p (non-clustered) and λ c λ p (clustered)serve as approximations for the mean number of planets per star,they should not be overinterpreted. We use the true rates of planetsper star for more detailed calculations in §4. Non-clustered model:
Orbital periods are drawn independently from a simple power-lawdistribution: f ( P ) ∝ P α P , P min (cid:54) P (cid:54) P max , (2)where f ( P ) is the probability density function (PDF) and α P is thepower-law index for the period distribution. We choose P min = P max =
300 d.
Clustered periods and clustered periods and sizes models:
The orbital period of each planet is drawn conditionally on theperiod scale for its parent cluster, P c . For each cluster we draw atrial P c from a simple power law, and the planets in each clusterhave trial periods drawn from a log-normal distribution conditionedon P c with a cluster width that is scaled to the number of planets: f ( P c ) ∝ P c α P (3) P (cid:48) i ∼ Lognormal ( , N p σ P ) (4) P i = P c P (cid:48) i , P min (cid:54) P c (cid:54) P max (5)where P (cid:48) i are the unscaled periods, P i are the true periods, N p is thenumber of planets in the cluster, and σ P is the cluster width scalefactor, per planet in the cluster. The trial periods are accepted orrejected based on a heuristic for orbital instability, as described in§2.2. The mean of a ZTP-distributed random variable X with parameter λ isgiven by E [ X ] = λ − e − λ . In practice, we first draw unscaled periods P (cid:48) i ∼ Lognormal ( , N p σ P ) ,truncated between P (cid:48) min = (cid:112) P min / P max and P (cid:48) max = (cid:112) P max / P min (thisis to avoid drawing unscaled periods such that maximum period ratio isgreater than P max / P min , which would have no chance fitting between P min and P max regardless of P c ). We then draw a cluster period scale P c fromEquation (3) (truncated between P min / min { P (cid:48) i } and P max / max { P (cid:48) i } ) togive P i = P (cid:48) i P c , checking that the drawn value of P c allows the planetsin the cluster to be stable with all other planets in the system; see §2.2 for adetailed outline of this process. Non-clustered model and clustered periods model:
Planet radii are drawn independently from a broken power law: f ( R p ) ∝ (cid:40) R p α R , R p , min (cid:54) R p (cid:54) R p , break R p α R , R p , break < R p (cid:54) R p , max , (6)where f ( R p ) is the PDF, and α R and α R are the broken power-lawindices for the planet radius distribution below and above R p , break ,the break radius, respectively. Our choice of a broken power law forthe radius distribution is motivated by previous studies suggestingthat there is an observed break at around 2–3 R ⊕ , with a rise inoccurrence down to ∼ R ⊕ and a plateau below that (Youdin 2011;Howard et al. 2012; Petigura, Marcy, & Howard 2013b). Clustered periods and sizes models:
The radius of each planet is drawn conditionally on the radius scalefor its parent cluster, R p , c . The cluster radius scales R p , c are drawnfrom a broken power law, while the planets in each cluster have theirradii drawn from a log-normal distribution conditioned on R p , c witha fixed scale for the cluster width σ R : f ( R p , c ) ∝ (cid:40) R p , c α R , R p , min (cid:54) R p , c (cid:54) R p , break R p , c α R , R p , break < R p , c (cid:54) R p , max . (7) R p , i ∼ Lognormal ( R p , c , σ R ) . (8)For all three models, we limit our analyses to R p , min = . R ⊕ and R p , max = R ⊕ , as the distribution of larger or smaller planetsis not well constrained by Kepler observations.
Planets vary significantly in structure and composition (Wolfgang& Laughlin 2012; Lopez & Fortney 2014; Chen & Kipping 2016;Wolfgang, Rogers, & Ford 2016). We adopt a non-parametric mass–radius relation from Ning, Wolfgang, & Ghosh (2018) in order toassign planet masses probabilistically given their drawn planet radii.This mass-radius relation involves a series of Bernstein polynomialsused to model the joint distribution of mass and radius based on asample of 127
Kepler exoplanets with RV or TTV masses, and ismore flexible than simpler, power law mass–radius relations. Toaccelerate computations, we pre-compute and use a 1001 × We draw the orbital eccentricities from a Rayleigh distribution, e ∼ Rayleigh ( σ e ) , where the Rayleigh parameter ( σ e ) defines thescale for eccentricities, following previous works (e.g. Moorhead etal. 2011; Fabrycky et al. 2014). Finally, we allow for two populations of planetary systems withseparate distributions of mutual inclinations between planets. Weuse Rayleigh distributions for both populations, as the Rayleigh dis-tribution has been used to describe mutual inclinations in previousworks (e.g. Fabrycky & Winn 2009; Lissauer et al. 2011b; Fang &Margot 2012; Dawson, Lee & Chiang 2016). Thus, for a fraction f σ i , high of systems we draw from a broad distribution of mutual MNRAS000
Kepler exoplanets with RV or TTV masses, and ismore flexible than simpler, power law mass–radius relations. Toaccelerate computations, we pre-compute and use a 1001 × We draw the orbital eccentricities from a Rayleigh distribution, e ∼ Rayleigh ( σ e ) , where the Rayleigh parameter ( σ e ) defines thescale for eccentricities, following previous works (e.g. Moorhead etal. 2011; Fabrycky et al. 2014). Finally, we allow for two populations of planetary systems withseparate distributions of mutual inclinations between planets. Weuse Rayleigh distributions for both populations, as the Rayleigh dis-tribution has been used to describe mutual inclinations in previousworks (e.g. Fabrycky & Winn 2009; Lissauer et al. 2011b; Fang &Margot 2012; Dawson, Lee & Chiang 2016). Thus, for a fraction f σ i , high of systems we draw from a broad distribution of mutual MNRAS000 , 1–31 (2019) lustered Planetary Systems inclinations while for the remaining 1 − f σ i , high fraction of systemswe draw from a narrower distribution: i m ∼ (cid:40) Rayleigh ( σ i , high ) , u < f σ i , high Rayleigh ( σ i , low ) , u (cid:62) f σ i , high , (9)where u ∼ Unif ( , ) and σ i , low (cid:54) σ i , high .In addition, we check whether each planet is near a 2:1, 3:2, 4:3,or 5:4 MMR with another planet, and draw i m ∼ Rayleigh ( σ i , low ) for planets that are. For our purposes, we regard adjacent planet pairsas “near an MMR” if their period ratio is within 5% exterior to theMMR. This amounts to period ratios between 1.5 and 1.575 as nearthe 3:2 MMR and period ratios between 2 and 2.1 as near the 2:1MMR, for example. Our motivation for forcing planets near MMRsto have more coplanar orbits than other planets is to explore whetherthis alone can produce the apparent relative excesses of planets withperiod ratios just wide of MMRs (Fabrycky et al. 2014; Steffen &Hwang 2015), or whether it is essential to generate more planetsnear MMRs than what is drawn naturally from our model for orbitalperiods (see §2.1.2). Since planets with more coplanar orbits aremore likely to transit together than highly mutually inclined planets,our model has the effect of producing slightly more observed planetpairs near MMRs than planet pairs with arbitrary period ratios.To test the robustness of our conclusions to this special treat-ment of the near-MMR planets, we also explore a model without thelowering of mutual inclinations for such planets, using our clusteredperiods and sizes model. We refer to this model as the “alternativeMMR inclinations” model for the remainder of the paper. The resultsof this model will be primarily discussed in §3 and §4.4.We leave f σ i , high , σ i , high , and σ i , low as free parameters of themodels. We test whether planetary systems are likely to be long-term Hillstable after drawing their periods. When a planetary system is iden-tified as likely unstable, then it is discarded and we attempt to redraworbital periods. Our instability criterion is based on the spacing be-tween adjacent planets ( ∆ ) normalized by the mutual Hill radius( R H ), which is given by: R H = (cid:18) a in + a out (cid:19) (cid:20) m in + m out M (cid:63) (cid:21) / , (10)where a in , a out are the semimajor axes and m in , m out are the massesof the inner and outer planets, respectively, and M (cid:63) is the mass of thestellar host. We define an instability criterion that is parametrizedin terms of R H : ∆ = a out ( − e out ) − a in ( + e in ) R H . (11)To avoid generating planetary systems likely to be unstable, we re-quire ∆ (cid:62) ∆ c where ∆ c is the minimum separation (held fixed). Fortwo circular and coplanar orbits, the minimum separation required tobe Hill stable (i.e., no close encounters) is given by ∆ > √ (cid:39) . ∆ <
10 are virtually always unstable. Pu & Wu (2015)find that ∆ (cid:38)
12 is required for long-term stability of planets inmultiplanet systems given certain assumptions about the massesand eccentricity distribution for the
Kepler exoplanets. In our earlyexploratory analyses where ∆ c was treated as a free parameter, wefound that our models preferred smaller values, ∆ (cid:46)
10. Therefore, we set ∆ c = Our non-clustered model has a total of 10 free parameters: λ p , α P , α R , α R , R p , break , σ e , ∆ c , f σ i , high , σ i , high , and σ i , low . Theclustered periods model has an additional two parameters λ c and σ P (the cluster width in log-period per planet in the cluster) to give12 free parameters, while the clustered periods and sizes modeladds yet another parameter σ R (the cluster width in log-radius) fora total of 13 free parameters. In our early exploratory analysis wefound that the break radius R p , break and the minimum separation ∆ c were not well constrained by Kepler observations (see §2.5 formore details). In particular, R p , break can take on a wide range ofvalues and the models seem to prefer smaller values of ∆ c (cid:46) R p , break = R ⊕ and ∆ c = physical catalogue We outline the procedure for generating an underlying population ofplanetary systems (a physical catalogue ) from the clustered periodsand sizes model as follows. The procedure for the other two modelsare analogous (we explain how to modify the procedure after thesesteps).(i) Set a number of target stars N stars , sim for each simulated cat-alogue, typically equal to the number of Kepler targets being usedas observational constraints.(ii) Set a value for each of the model parameters.(iii) For each target, assign stellar properties drawn from the
Ke-pler stellar catalogue (updated based on Gaia DR2, as describedin Hsu et al. 2019, and allowing for the uncertainties in the stellarproperties). Stellar properties include radius and mass which affectthe observed transit properties, as well as the one-sigma depth func-tion, window function, and contamination that are necessary for theplanet detection model (Hsu et al. 2019).(iv) Draw a number of clusters in the system, N c ∼ Poisson ( λ c ) .Re-sample until N c (cid:54) N c , max .(v) For each cluster:(a) Draw a number of planets in the cluster from a zero-truncated Poisson (ZTP) distribution, N p ∼ ZTP ( λ p ) . Re-sampleuntil N p (cid:54) N p , max .(b) Draw a radius for each planet in the cluster: first, draw acharacteristic radius R p , c for the cluster according to Equation(7). If N p =
1, the radius of the cluster’s one planet is simply R p = R p , c . If N p >
1, draw a radius for each of the cluster’splanets, R p , i ∼ Lognormal ( R p , c , σ R ) , where i = , ..., N p (thelog is base- e ). Draw their masses using the mass–radius relationfrom the non-parametric model in Ning, Wolfgang, & Ghosh(2018).(c) Draw an orbital eccentricity e ∼ Rayleigh ( σ e ) and argu-ment of periastron ω ∼ Unif ( , π ) for each planet in the cluster.(d) Draw orbital periods for each planet in the cluster: If N p =
1, assign an unscaled period of P (cid:48) =
1. If N p >
1, draw their
MNRAS , 1–31 (2019)
He, Ford, and Ragozzine
Table 1.
List describing the important model parameters, and how they are used, in this paper.Parameter Definition of parameter Equation Relevant quantity/distribution N stars , sim Number of simulated stars - - λ c Mean number of clusters per star (before rejection sampling) - N c ∼ Poisson ( λ c ) N c , max Maximum number of clusters per star - N c (cid:54) N c , max λ p Mean number of planets per system (non-clustered model) or cluster (clusteredmodels) (before rejection-sampling) - N p ∼ Poisson ( λ p ) or ZTP ( λ p ) N p , max Maximum number of planets per cluster - N p (cid:54) N p , max α P Power law index for distribution of periods (non-clustered model) or period scales(clustered models) (2), (3) P c ∼ f ( P c ) ∝ P c α P P min Minimum period (days) (2), (3) P = P (cid:48) P c (cid:62) P min P max Maximum period (days) (2), (3) P = P (cid:48) P c (cid:54) P max α R Power law index for planetary radius distribution for R p (cid:54) R p , break (6), (7) R p , c ∼ f ( R p ) ∝ R p α R , R p (cid:54) R p , break α R Power law index for planetary radius distribution for R p > R p , break (6), (7) R p , c ∼ f ( R p ) ∝ R p α R , R p > R p , break R p , break Break radius for planetary radii ( R ⊕ ) (6), (7) R p , min < R p , break < R p , max R p , min Minimum planetary radius ( R ⊕ ) (6), (7) R p (cid:62) R p , min R p , max Maximum planetary radius ( R ⊕ ) (6), (7) R p (cid:54) R p , max σ P Scale factor, per planet in the cluster, for the cluster (unscaled) period distribution (4) P (cid:48) ∼ Lognormal ( , N p σ P ) σ R Scale factor for the cluster radius distribution (8) R p ∼ Lognormal ( R p , c , σ R ) σ e Scale for orbital eccentricities - e ∼ Rayleigh ( σ e ) ∆ c Minimum allowed separation between adjacent planets in mutual Hill radii - ∆ (cid:62) ∆ c f σ i , high Fraction of systems with relatively high mutual inclinations, σ i , high (9) - σ i , high Scale for mutual inclinations for systems with high mutual inclinations (deg) (9) i m ∼ Rayleigh ( σ i , high ) σ i , low Scale for mutual inclinations for systems with low mutual inclinations (deg) (9) i m ∼ Rayleigh ( σ i , low ) N attempts Maximum number of attempts for re-sampling periods and period scales for eachcluster - - unscaled periods P (cid:48) i ∼ Lognormal ( , N p σ P ) , where i = , ..., N p (the log is base- e ), and sort them in increasing order. Check if ∆ = [ a i + ( − e i + )− a i ( + e i )]/ R H ( i , i + ) (cid:62) ∆ c for all i in the cluster.Re-sample the unscaled periods P (cid:48) i until this condition is satisfiedor the maximum number of attempts N attempts is reached. If thecase is the latter, discard the cluster. Draw a period scale factor P c (d) according to Equation (2) and multiply each planet’s unscaledperiods by the period scale for its parent cluster: P i = P (cid:48) i P c .(e) Test for stability: Check if ∆ (cid:62) ∆ c for all adjacent planetpairs in the entire system, including planets from previouslydrawn clusters. If the system is identified as likely unstable, re-draw P c for the current cluster until this condition is satisfied oruntil the maximum number of attempts N attempts is reached. Ifthe latter case occurs, discard the cluster.(vi) For each system, draw an inclination angle (relative to theplane of the sky) for the reference plane of the system isotropically,cos i ref ∼ Unif ( , ) . For each system, draw a number u ∼ Unif ( , ) .If u < f σ i , high , set σ i = σ i , high ; otherwise, set σ i = σ i , low . Assignan orbit to each planet in the system:(a) Compute the period ratios P = P i + / P i of all adjacentplanet pairs in the system. For each planet, check if it is near anyMMRs with any adjacent planet by checking if P mmr (cid:54) P (cid:54) We adopt a value of N attempts = N p with a small set value of σ P ) clusters from significantlyincreasing the computational time for simulating a physical catalogue . Onthe other hand, if we attempt the drawing of planets in each cluster only a fewtimes (or even just once), the rejection-sampling algorithm tends to producetoo few planetary systems. Thus, a modest number of attempts per clusterprovides a compromise between these two extremes. Also, a consequence ofthis procedure is that the true mean rates of clusters and planets per clusterin our samples are somewhat less than λ c and E [ N p ] , respectively. . P mmr for any P mmr in { , . , / , . } (aside from theinner- and outer-most planet, each planet is part of two adjacentplanet pairs). If the planet is not near any MMRs with anotheradjacent planet, draw a mutual inclination angle (relative to thereference plane) i m ∼ Rayleigh ( σ i ) ; otherwise, draw a mutualinclination angle i m ∼ Rayleigh ( σ i , low ) regardless of whether σ i , high or σ i , low was assigned to the non-resonant planets of thesystem.(b) For each planet’s orbit, draw an angle of ascending node Ω ∼ Unif ( , π ) in the reference plane.(c) Compute the inclination angle i (relative to the plane ofthe sky) for each planet’s orbit using the spherical law of cosines,cos ( i ) = cos ( i ref ) cos ( i m ) + sin ( i ref ) sin ( i m ) cos ( Ω ) .The procedure is very similar for the other models. For theclustered periods model, instead of drawing a characteristic radius R p , c for each cluster as in Step 5(b), we simply assign planet radiiby drawing them directly from Equation (6). In the case of the non-clustered model, Steps 4–5 are simplified to drawing a number ofplanets N p ∼ Poisson ( λ p ) and directly drawing their periods andradii by Equations (2) and (6). The non-clustered model could be considered as a special case of theclustered models by setting the number of planets per cluster to one, N p = N c . In practice, for the non-clustered model, wedraw periods for all the N p planets simultaneously and accept or reject allperiods at once, so as to improve computational efficiency and to minimizeartefacts near P min and P max for systems with many planets.MNRAS000
List describing the important model parameters, and how they are used, in this paper.Parameter Definition of parameter Equation Relevant quantity/distribution N stars , sim Number of simulated stars - - λ c Mean number of clusters per star (before rejection sampling) - N c ∼ Poisson ( λ c ) N c , max Maximum number of clusters per star - N c (cid:54) N c , max λ p Mean number of planets per system (non-clustered model) or cluster (clusteredmodels) (before rejection-sampling) - N p ∼ Poisson ( λ p ) or ZTP ( λ p ) N p , max Maximum number of planets per cluster - N p (cid:54) N p , max α P Power law index for distribution of periods (non-clustered model) or period scales(clustered models) (2), (3) P c ∼ f ( P c ) ∝ P c α P P min Minimum period (days) (2), (3) P = P (cid:48) P c (cid:62) P min P max Maximum period (days) (2), (3) P = P (cid:48) P c (cid:54) P max α R Power law index for planetary radius distribution for R p (cid:54) R p , break (6), (7) R p , c ∼ f ( R p ) ∝ R p α R , R p (cid:54) R p , break α R Power law index for planetary radius distribution for R p > R p , break (6), (7) R p , c ∼ f ( R p ) ∝ R p α R , R p > R p , break R p , break Break radius for planetary radii ( R ⊕ ) (6), (7) R p , min < R p , break < R p , max R p , min Minimum planetary radius ( R ⊕ ) (6), (7) R p (cid:62) R p , min R p , max Maximum planetary radius ( R ⊕ ) (6), (7) R p (cid:54) R p , max σ P Scale factor, per planet in the cluster, for the cluster (unscaled) period distribution (4) P (cid:48) ∼ Lognormal ( , N p σ P ) σ R Scale factor for the cluster radius distribution (8) R p ∼ Lognormal ( R p , c , σ R ) σ e Scale for orbital eccentricities - e ∼ Rayleigh ( σ e ) ∆ c Minimum allowed separation between adjacent planets in mutual Hill radii - ∆ (cid:62) ∆ c f σ i , high Fraction of systems with relatively high mutual inclinations, σ i , high (9) - σ i , high Scale for mutual inclinations for systems with high mutual inclinations (deg) (9) i m ∼ Rayleigh ( σ i , high ) σ i , low Scale for mutual inclinations for systems with low mutual inclinations (deg) (9) i m ∼ Rayleigh ( σ i , low ) N attempts Maximum number of attempts for re-sampling periods and period scales for eachcluster - - unscaled periods P (cid:48) i ∼ Lognormal ( , N p σ P ) , where i = , ..., N p (the log is base- e ), and sort them in increasing order. Check if ∆ = [ a i + ( − e i + )− a i ( + e i )]/ R H ( i , i + ) (cid:62) ∆ c for all i in the cluster.Re-sample the unscaled periods P (cid:48) i until this condition is satisfiedor the maximum number of attempts N attempts is reached. If thecase is the latter, discard the cluster. Draw a period scale factor P c (d) according to Equation (2) and multiply each planet’s unscaledperiods by the period scale for its parent cluster: P i = P (cid:48) i P c .(e) Test for stability: Check if ∆ (cid:62) ∆ c for all adjacent planetpairs in the entire system, including planets from previouslydrawn clusters. If the system is identified as likely unstable, re-draw P c for the current cluster until this condition is satisfied oruntil the maximum number of attempts N attempts is reached. Ifthe latter case occurs, discard the cluster.(vi) For each system, draw an inclination angle (relative to theplane of the sky) for the reference plane of the system isotropically,cos i ref ∼ Unif ( , ) . For each system, draw a number u ∼ Unif ( , ) .If u < f σ i , high , set σ i = σ i , high ; otherwise, set σ i = σ i , low . Assignan orbit to each planet in the system:(a) Compute the period ratios P = P i + / P i of all adjacentplanet pairs in the system. For each planet, check if it is near anyMMRs with any adjacent planet by checking if P mmr (cid:54) P (cid:54) We adopt a value of N attempts = N p with a small set value of σ P ) clusters from significantlyincreasing the computational time for simulating a physical catalogue . Onthe other hand, if we attempt the drawing of planets in each cluster only a fewtimes (or even just once), the rejection-sampling algorithm tends to producetoo few planetary systems. Thus, a modest number of attempts per clusterprovides a compromise between these two extremes. Also, a consequence ofthis procedure is that the true mean rates of clusters and planets per clusterin our samples are somewhat less than λ c and E [ N p ] , respectively. . P mmr for any P mmr in { , . , / , . } (aside from theinner- and outer-most planet, each planet is part of two adjacentplanet pairs). If the planet is not near any MMRs with anotheradjacent planet, draw a mutual inclination angle (relative to thereference plane) i m ∼ Rayleigh ( σ i ) ; otherwise, draw a mutualinclination angle i m ∼ Rayleigh ( σ i , low ) regardless of whether σ i , high or σ i , low was assigned to the non-resonant planets of thesystem.(b) For each planet’s orbit, draw an angle of ascending node Ω ∼ Unif ( , π ) in the reference plane.(c) Compute the inclination angle i (relative to the plane ofthe sky) for each planet’s orbit using the spherical law of cosines,cos ( i ) = cos ( i ref ) cos ( i m ) + sin ( i ref ) sin ( i m ) cos ( Ω ) .The procedure is very similar for the other models. For theclustered periods model, instead of drawing a characteristic radius R p , c for each cluster as in Step 5(b), we simply assign planet radiiby drawing them directly from Equation (6). In the case of the non-clustered model, Steps 4–5 are simplified to drawing a number ofplanets N p ∼ Poisson ( λ p ) and directly drawing their periods andradii by Equations (2) and (6). The non-clustered model could be considered as a special case of theclustered models by setting the number of planets per cluster to one, N p = N c . In practice, for the non-clustered model, wedraw periods for all the N p planets simultaneously and accept or reject allperiods at once, so as to improve computational efficiency and to minimizeartefacts near P min and P max for systems with many planets.MNRAS000 , 1–31 (2019) lustered Planetary Systems observed catalogue The underlying population of planetary systems is then used to gen-erate an observed catalogue of exoplanets using a procedure whichsimulates the observational pipeline employed by
Kepler . This pro-cedure is the same regardless of the choice of model used for gener-ating the physical catalogue . The full pipeline that is implementedin SysSim as described in Hsu et al. (2018, 2019) incorporates many
Kepler
DR25 data products, including tabulated window functions,one sigma depth functions, and a detection efficiency derived fromanalysing pixel-level transit injections. We briefly summarize themain considerations here.First, we reduce the number of planets by only keeping (in theobserved catalogue) planets that transit. This amounts to requiring that the impact parameter b is lessthan 1 + R p / R (cid:63) , b = a cos iR (cid:63) − e + e sin ω < + R p / R (cid:63) . (12)Next, we select a subset weighted by their detection probabilities.We require at least three transits to be observed by the Kepler mission based on window functions provided as part of
Kepler
DR25, as described in Hsu et al. (2019). A planet is labelled asdetected if a random number u ∼ Unif ( , ) is less than the planet’sdetection probability (conditioned on it transiting), as calculated byassuming the joint detection and vetting efficiency model describedin Hsu et al. (2019).For each of the detected transiting planets, we compute thetrue transit depths accounting for limb darkening and the true transitdurations according to Kipping (2010), equation (15) therein: t dur = P π ρ c √ − e sin − (cid:18) √ − b a R ρ c sin i (cid:19) (13) (cid:39) P π a R √ − e + e sin ω (cid:112) − b , (14) ρ c ≡ − e + e sin ω, (15)where t dur is the full-width half-maximum transit duration and a R ≡ a / R (cid:63) is the semimajor axis in units of the stellar radius. The one-term analytic expression in Equation (13) only assumes that theplanet–star separation does not change during the entire transit eventand neglects the planet size (Kipping 2010). Using this definitiongrazing transits (i.e., 1 (cid:54) b < + R p / R (cid:63) ) have zero duration.Finally, we add measurement noise to the true period, transitdepth, and transit durations and compile the “observed” proper-ties of the detected transiting planets into a simulated observedcatalogue. Some of the simulations used for exploring parameterspace and training the emulator (see §2.5-§2.6) use a simplifiednoise model where the fractional uncertainties in orbital periods,transit depths, and transit durations are held fixed. However, thefinal simulations for constraining model parameters and inferringthe distributions of observed and physical properties use a diagonalversion of the transit noise model of Price & Rogers (2014), whichis based on a trapezoidal transit model and accounts for the finiteintegration time. We use SysSim in the single-observer mode, and do not make use ofCORBITS for sky-averaging (Brakensiek & Ragozzine 2016). Working insingle-observer model allows our forward model to reproduce the variationsin the observed
Kepler catalogue due to the finite number of targets (seeHsu et al. 2019) and avoids issues with two-planet correlations that are notcalculated by CORBITS.
Next, we compare the simulated observed catalogues of exoplanetsderived from our models to the actual observed population of ex-oplanets by the
Kepler mission. In principle, one could attempt togenerate simulated catalogues that precisely match all planet prop-erties to those of the actual
Kepler catalogue. Of course, the oddsof generating such a catalogue are minuscule, even if one coulduse the perfect model for the underling exoplanet population dueto the stochastic nature of the model. Instead, we identify a set ofsummary statistics that encode the most physically important prop-erties of the distributions of planetary systems observed by
Kepler in §2.4.1. Next, we define a distance function that quantifies thedegree of similarity between the summary statistics for the two cat-alogues in §2.4.2. We describe a procedure for identifying sets ofparameters for our physical models that approximately minimizethe distance between simulated and observed catalogues (allowingfor inevitable shot noise due to the stochastic nature of the model) in§2.5. Finally, we constrain the model parameters using approximateBayesian Computing (ABC). Given the cost of the full model, wetrain a Gaussian process emulator to approximate the distribution ofthe distances computed from our forward model in §2.6. We obtainsamples from the ABC posterior by drawing trial model parametersbased on the GP emulator, computing the distance with the fullSysSim forward model, and accepting parameter values that resultin a small distance between the simulated observed catalogue andthe
Kepler
DR25 catalogue.
Our procedure for generating observed catalogues yields an ob-served catalogue with a “measured” orbital period P , transit duration t dur , and transit depth δ for each observed planet. A good forwardmodel must result in a similar number of detected planets, as wellas a similar number of systems with m detected planets. Addition-ally, we want our models to reproduce the observed distributions oforbital periods, transit depths and transit durations. Finally, we wantour forward models to produce planetary systems that have realisticcorrelations between the orbital periods and sizes of planets withinthe same system. Therefore, we compute the following summarystatistics for each observed catalogue:(i) the overall rate of observed planets relative to the number oftarget stars, f = N p , tot / N stars , where N p , tot and N stars are the totalnumbers of observed planets and stars in the catalogue, respectively,(ii) the observed multiplicity distribution, { N m } , where N m isthe number of systems with m observed planets and m = , , , ... ,(iii) the observed orbital period distribution, { P } ,(iv) the observed distribution of period ratios of apparently ad-jacent planets, {P} , where P = P i + / P i ,(v) the observed transit depth distribution, { δ } ,(vi) the observed distribution of the transit depth ratios of appar-ently adjacent planets, { δ i + / δ i } ,(vii) the observed transit duration distribution, { t dur } ,(viii) the observed distribution of (period-normalized) transit du-ration ratios of apparently adjacent planets near mean motion reso-nances, { ξ res } , and(ix) the observed distribution of (period-normalized) transit du-ration ratios of apparently adjacent planets not near mean motionresonances, { ξ non − res } .Previously published studies have attempted to match a subsetof these summary statistics, but not all at once. Each summary statis- MNRAS , 1–31 (2019)
He, Ford, and Ragozzine tic is most sensitive to one or two model parameters, but typicallyhave weaker dependencies on other model parameters.For example, the transit duration distribution is useful forcharacterizing the orbital eccentricities (Ford, Quinn, & Veras2008). Moreover, the period-normalized transit duration ratio ξ (also known as the orbital-velocity normalized transit duration ra-tio; Fabrycky et al. 2014) is a useful probe of the orbital propertiesof exoplanets in multitransiting systems. It can be shown that theratio of the period to the transit duration cubed is proportional to thestellar density, P / t ∝ ρ (cid:63) (Seager & Mallén-Ornelas 2003). Thus,for two planets transiting the same star, it is useful to consider theratio of their period-normalized transit durations. This quantity hasbeen used to test whether two planets really transit the same star; ifthey do, ξ should be near unity (Lissauer et al. 2012; Fabrycky et al.2014). Following Steffen et al. (2010) and Fabrycky et al. (2014),we define ξ as: ξ = (cid:18) t dur , in t dur , out (cid:19) (cid:18) P out P in (cid:19) / . (16)More importantly for our purposes, the distribution of ξ encodesinformation about the orbital eccentricities and impact parameters(i.e., inclination angles) of the transiting planets (Fabrycky et al.2014; Morehead 2016). It is useful to transform this quantity bytaking the logarithm, log ξ , so that negative values imply ξ < ξ >
1. For planets in circular, coplanar orbitsaround the same star, log ξ must be non-negative (if they both transitat the equator, b = ξ =
1; otherwise, the inner planet musthave a smaller impact parameter b than the outer planet, and thus alonger period-normalized transit duration, ξ > ξ , because the impact parameters become more randomizedas the outer transiting planets need not necessarily have larger valuesof b than the inner planets. For eccentric orbits, the distribution oflog ξ becomes more spread out (i.e., more extreme values of ξ aremore common) with increasing eccentricities because the velocitiesof the planets during transit become more randomized, dependingon whether the transits occur near periastron or apastron.Our forward model involves assigning systems to one oftwo separate mutual inclination scales and assigning planets nearMMRs to follow the smaller mutual inclination distribution ( i m ∼ Rayleigh ( σ i , low ) ; see Step 6(a) in §2.2). In order to constrain theinclinations of both populations, we include two summary statisticsbased on the ξ distribution: { ξ res } calculated using only observedplanet pairs apparently near a MMR (i.e. based on the observedperiods) and { ξ non − res } calculated using the remaining planet pairs(i.e., planets pairs not apparently near any MMRs). ABC is a powerful method for performing inference on modelswhere it is impractical to write an explicit likelihood, such as thecase for studying multiplanet systems observed by
Kepler (Hsu et al.2018). In ABC, one must define a distance function to quantify howdifferent two catalogues are. This distance is a function of only thesummary statistics for each catalogue and goes to zero when the twosummary statistics are identical. For simple analytic distributions,one can identify sufficient summary statistics for which one canrigorously prove that the ABC posterior approaches the true poste-rior as a distance threshold goes to zero. In practice, ABC is mostuseful for complex problems like ours, for which it is impractical toidentify sufficient statistics. In these cases, the ABC posterior will be broader than the true posterior, since some draws may reasonablyreproduce the summary statistics, but differ in some way that didnot contribute to the distance function.Having identified summary statistics that are both observableand physically significant in §2.4.1, we proceed to define a dis-tance function for each summary statistic. For most of our summarystatistics, we use either the two-sample Kolmogorov–Smirnov (KS;Kolmogorov 1933; Smirnov 1948) distance or a rescaled variant ofthe two-sample Anderson–Darling (AD; Anderson & Darling 1952;Pettitt 1976) distance. A component distance of zero would meanthat the marginal distributions for each of these summary statisticsare identical. The use of KS distances rewards distributions thatagree best near the median of the marginal distribution, while theuse of AD distances is more sensitive to differences in the tails ofthe marginal distributions. We perform calculations using both, pri-marily as a sensitivity test. By comparing results using two differentdistance functions, we can check whether any of our conclusionsare sensitive to the choice of distance function.For our total distance, we take a linear combination of individ-ual distance terms for each summary statistic: D W = (cid:34)(cid:213) i (cid:48) w i (cid:48) (cid:107)D i (cid:48) (cid:107) α D (cid:35) / α D , (17)where D W is the total weighted distance, we set α D =
1, and w i (cid:48) = / ˆ σ (D i (cid:48) ) is the weight (defined below) of the i (cid:48) th distanceterm, D i (cid:48) . In principle, one could have chosen another value of α D ,corresponding to the Euclidean normal or maximum norm. Eitherof these would result in a total distance that is more sensitive tothe most discrepant summary statistic than our choice of α D = Kepler populationdescribed in §2.4.1, as well as the overall rate of planets per observedstar, f Kepler . We aim to assign weights to the individual distancesthat reflect the precision with which they were measured by
Kepler .Therefore, we specialize Equation 17 to: D W = D f ˆ σ ( D f ) + D mult ˆ σ ( D mult ) + (cid:213) i = D i ˆ σ (D i ) , (18)where D f is the distance between the two rates of observed planets, D mult is the distance of the multiplicity distributions, D i is thedistance between the distributions of the i th observable, and ˆ σ (D) is the estimated root-mean-square (RMS) of the distance D givena “perfect” model (i.e., comparisons of simulated catalogues toa reference simulated catalogue with the same model parameters).The indices for the observables run from i = i =
7, denoting thefollowing distributions of measured properties in order: period { P } ,period ratio {P} , transit depth { δ } , transit depth ratio { δ i + / δ i } ,transit duration { t dur } , and period-normalized transit duration ratiofor near-resonance { ξ res } and not-near-resonance { ξ non − res } planetpairs. Rate of Observed Planets:
Since the rate of observed planetsis a scalar, the distance function is simply, D f = | f sim − f Kepler | . Multiplicity:
The observed multiplicities { N m } is a vector ofintegers. If each system had the same probability to be observedas an m -planet system, then it would be a multinomial distribution.Formally, { N m } is not drawn from a multinomial since each systemhas a different set of probabilities for being observed as an m -planetsystem. Nevertheless, the theory of multinomial distribution is use-ful for identifying an appropriate distance function for the observedmultiplicities. The Cressie–Read power divergence (CRPD) statistic MNRAS000
The observed multiplicities { N m } is a vector ofintegers. If each system had the same probability to be observedas an m -planet system, then it would be a multinomial distribution.Formally, { N m } is not drawn from a multinomial since each systemhas a different set of probabilities for being observed as an m -planetsystem. Nevertheless, the theory of multinomial distribution is use-ful for identifying an appropriate distance function for the observedmultiplicities. The Cressie–Read power divergence (CRPD) statistic MNRAS000 , 1–31 (2019) lustered Planetary Systems (Cressie & Read 1984) is commonly used in comparing multino-mial distributions like the multiplicity distribution and is knownto be more robust than other choices like χ for cases like ourswhere there is a large dynamic range between the values in eachcategory and one of the categories ( O ) often has values less than5. Therefore, for D mult , we adopt the CRPD statistic which is givenby: ρ CRPD = (cid:213) j O j (cid:20)(cid:18) O j E j (cid:19) / − (cid:21) , (19)where O j is the number of “observed” systems according to one ofour models and E j is the number of expected systems in the j th bin based on the rate of such systems in the Kepler data set. Theindices j label the bins corresponding to 1, 2, 3, 4, and 5 + observedplanets. Note that by definition, this statistic only sums over bins j for which E j (cid:44) E j = O j > D mult meets the primary goal of clearly favouringbetter matches to the Kepler data. Continuous Distributions:
For remaining D i terms, we per-form the full analysis with two different distance functions for com-paring samples of continuous random variables: the two-sampleKolmogorov–Smirnov (KS) distance, and a rescaled variant of thetwo-sample Anderson–Darling (AD) distance. For two finite sam-ples of sizes n and m , described by empirical distribution functions F n ( x ) and G m ( x ) , the two-sample KS and AD distances are definedas: D KS = max | F n ( x ) − G m ( x )| , (20) D AD ≡ nmN ∫ ∞−∞ [ F n ( x ) − G m ( x )] H N ( x )[ − H N ( x )] dH N ( x ) , (21) = nm N − (cid:213) i = ( M i N − ni ) i ( N − i ) , (22)where N = n + m is the combined sample size, H N ( x ) = [ nF n ( x ) + mG m ( x )]/ N is the empirical distribution function for the combinedsample, and M i is the number of observations less than or equal tothe i th smallest in the combined sample (Pettitt 1976).The KS distance is simply the maximum absolute differencein the cumulative distribution functions (CDFs) and is thus mostsensitive to the difference in the bulk locations of the two distribu-tions. The AD distance on the other hand, more heavily weights thetails of the distributions and is thus more sensitive to differencesat the extremes of the distributions. In practice, we find that thestandard AD distance given by Equation (22) does not sufficientlyaccount for the differences in the sample sizes n and m ; in other Lissauer et al. (2011a) use the exact multinomial statistic, but using thisstatistic is intractable for our case. After significant testing of different statis-tics, CRPD was shown to be the best approximation of the exact statistic fordistributions similar to ours. We also considered switching the interpretation of O j and E j , i.e. let O j denote the number of actual observed Kepler systems and E j denote thenumber of expected observed systems given a model. During exploratoryanalyses, we found that ρ CRPD would occasionally give infinite values ifa simulated catalogue included zero 5+ planet systems. Eliminating thatterm from the sum is not viable, as this can result in negative values, whichviolates the non-negative property of a distance function and would leadto favouring models with fewer multiplanet systems than were observed by
Kepler . words, two samples can give a very small AD distance even when n and m are vastly different. In the context of our models, this hasthe consequence that very small simulated observed catalogues (i.e.with only a handful of observed planets) can still result in small ADdistances, even enough to counteract larger distances of D f for theoverall rate of planets per star. This is clearly undesirable as suchmodels are terrible fits to the
Kepler observations. In order to coun-teract this unintended consequence, we modify the AD distance topenalize such models (i.e. those that result in almost no observedplanets) by dividing out the nm / N constant in front of the integral,giving: D AD (cid:48) = ∫ ∞−∞ [ F n ( x ) − G m ( x )] H N ( x )[ − H N ( x )] dH N ( x ) (23) = N ( nm ) N − (cid:213) i = ( M i N − ni ) i ( N − i ) . (24)For the remainder of the paper, we will refer to our modified ADdistance given by Equation (24) as simply “the AD distance” unlessotherwise noted. We refer to the total weighted distance functionsinvolving the KS or rescaled AD distance terms as D W , KS and D W , AD (cid:48) , respectively.In order to compute the weights w i (cid:48) for the individual dis-tance terms, we first generate a simulated observed catalogue usingthe clustered periods and sizes model which serves as a referencecatalogue, and then repeatedly simulate catalogues using the samemodel (i.e. with identical model parameters) which would be a “per-fect” model. The RMS ( ˆ σ (D) ) of each individual distance term andthe weights (the inverse of the RMS, w = / ˆ σ ) computed this wayare listed in Table 2, while the model parameters used to gener-ate this reference catalogue are given in Table 3. The distances foreach term are not zero because the simulations involve Monte Carlonoise and only a finite number of planets are drawn per iteration.Indeed, the true Kepler catalogue of observed planets is finite insize. Thus, we simulate a reference catalogue that contains a similarnumber of observed planets given the same number of target starsas the
Kepler mission. However, for the purposes of computing theweights and optimizing the distance function to find the best-fittingmodel parameters (see §2.5), we wish to reduce stochastic noise inorder to improve our power to distinguish between different models.Therefore, we simulate catalogues from the “perfect” model withfive times as many stars to give observed catalogues that are fivetimes as large as that of our
Kepler sample, balancing the desire toreduce stochastic noise and computational time. We compute theweights with 1000 repeated simulations from the “perfect” model.By summing the distances weighted by their variations for a“perfect” model, the total weighted distance emphasizes the dis-tances of observables that are well characterized and prevents anysingle distance term with high stochastic noise from dominatingthe total distance function. Also, this implies that each individ-ual distance term included in Equation (18) typically contributes aweighted value of roughly 1 to the total distance in the case of usingthe perfect model. This is apparent if one considers the effect of the nmN term in Equation(21): let m be the number of observed planets in the Kepler data and n bethe number of observed planets given by the model. If n (cid:39) m , this term isroughly ∼ m /
2; but if n ∼ (cid:28) m , this term becomes roughly mm + ∼ m is relatively large, the nmN term is thus much smaller when n (cid:28) m than when n (cid:39) m .MNRAS , 1–31 (2019) He, Ford, and Ragozzine
Table 2.
Table of weights for the individual distance terms as computed froma reference clustered periods and sizes model. The weights w are computedfrom the root mean squares of the distances, ˆ σ (D) , and are shown here asrounded whole numbers simply for guidance purposes.Distance term ˆ σ (D) w = / ˆ σ D f D mult σ (D) w = / ˆ σ ˆ σ (D) w = / ˆ σ D (for { P } ) 0.0173 58 0.000448 2230 D (for {P } ) 0.0288 35 0.00113 882 D (for { δ } ) 0.0179 56 0.000480 2085 D (for { δ i + / δ i } ) 0.0299 33 0.000911 1098 D (for { t dur } ) 0.0213 47 0.000524 1907 D (for { ξ res } ) 0.0678 15 0.00652 153 D (for { ξ non − res } ) 0.0328 30 0.00121 827 We couple the
Kepler
DR25 stellar propertiescatalogue with the results of the second
Gaia data release (DR2)(Andrae et al. 2018; Gaia Collaboration et al. 2018) in order to takeadvantage of its significantly improved stellar parameters for a largefraction of the
Kepler target stars. These improvements (of primaryinterest, the stellar radii) result from the refined parallax and thusdistance measurements of the Gaia mission. Furthermore, the GaiaDR2 parallaxes allow for a cleaner sample of main–sequence targetstars thanks to a more precise positioning of the targets in color-luminosity space. Additionally, the astrometric information allowsfor identification of targets likely consisting of multiple stars withcomparable luminosity.We use the same
Kepler target list as in Hsu et al. (2019), whoidentified a clean sample of FGK main-sequence (M-S) stars foroccurrence rate studies. This involved performing a series of cutson the
Kepler
DR25 stellar table based on measurements reportedin the Gaia DR2 and updating the stellar radii with values reportedin the Gaia DR2. For full details, see §3.1 of Hsu et al. (2019).In order to minimize sensitivity to uncertainties in stellar radii,impact parameters, and the limb darkening model, we have chosen adistance function based on the distribution of transit depths insteadof the distribution of planet radii or planet–star radius ratios, sincethe measured transit depth does not depend on our knowledge ofthe stellar radius (like the planet radius) and is better modelled as aGaussian distribution than the planet–star radius ratio (due to effectsof limb darkening and covariance with impact parameter). Never-theless, the uncertainties in stellar radii still affect our simulationsvia the transit depths of the simulated planets.
Planet Catalogue:
We start from the
Kepler
Data Release 25(DR25) (Thompson et al. 2018)
Kepler
Objects of Interest (KOI)catalogue as the basis for our study. This table is most suitablefor a population study of the
Kepler exoplanets because it involvesuniform vetting of the Q1-Q17 light curves obtained by
Kepler ,and thus does not involve human biases across individual systems.It is derived from processing using the SOC pipeline release 9.3,and involves fully automated dispositioning of the Threshold Cross-ing Events (TCEs) using the
Kepler Robovetter (Thompson et al.2015; Twicken et al. 2016). Specifically, the automated procedureinvolves determining whether the TCEs are transit-like, and if so,tests whether there is any evidence of an eclipsing binary, shift in the T o t a l s y s t e m s Kepler P (days)10 N u m b e r Kepler P i + 1 / P i N u m b e r N u m b e r i + 1 / i N u m b e r t dur (mins)050100 N u m b e r N u m b e r Near MMRNot near MMR
Figure 2.
The
Kepler
Q1-Q17 DR25 population of KOIs vetted as planetcandidates by the robovetter, given our cuts as described in §2.4.3. Fromtop to bottom: histograms of the observed planet multiplicities, periods( P ), period ratios ( P ), transit depths ( δ ), transit depth ratios ( δ i + / δ i ),transit durations ( t dur ), and period-normalized transit duration ratios ( ξ ). Weseparate the ξ distribution into two, for planets near- and not-near-MMRs asdefined in §2.1.6 (most planets are not near an MMR). The period, periodratio, transit depth, and transit depth ratio bins are uniformly spaced in log,while the bins in the other histograms are linearly spaced.MNRAS000
Q1-Q17 DR25 population of KOIs vetted as planetcandidates by the robovetter, given our cuts as described in §2.4.3. Fromtop to bottom: histograms of the observed planet multiplicities, periods( P ), period ratios ( P ), transit depths ( δ ), transit depth ratios ( δ i + / δ i ),transit durations ( t dur ), and period-normalized transit duration ratios ( ξ ). Weseparate the ξ distribution into two, for planets near- and not-near-MMRs asdefined in §2.1.6 (most planets are not near an MMR). The period, periodratio, transit depth, and transit depth ratio bins are uniformly spaced in log,while the bins in the other histograms are linearly spaced.MNRAS000 , 1–31 (2019) lustered Planetary Systems in-transit centroid position, or contamination from another source(Mullally et al. 2015; Coughlin et al. 2016). TCEs that pass all thetests are dispositioned as planetary candidates, and their planetaryand orbital parameters are computed using a Markov Chain MonteCarlo (MCMC) fitting algorithm (Rowe et al. 2015).Starting from the Kepler
DR25 KOI catalogue, we removeplanet candidates around targets that were excluded from the stellarcatalogue, as described above. Next, we keep only KOIs designatedas planet candidates by the
Kepler Robovetter . Then, we replace thetransit depths and durations in the
Kepler
DR25 catalogue (whichwere maximum likelihood estimators) with the median values fromthe MCMC-based posterior samples described in Rowe et al. (2015).We update the planet radii based on the observed transit depths, theupdated stellar radii from
Gaia
DR2, and the limb darkening modelfrom the DR25 stellar catalogue. Finally, we limit the planetarycatalogue to include planet candidates with periods between P min = P max =
300 d (see §2) and with updated planet radii between R p , min = . R ⊕ and R p , max = R ⊕ .These cuts result in a total of 79 935 Kepler targets (hereafterdenoted by N stars , Kep ), with 2137 total planet candidates in 1561systems (with periods between 3 and 300 d and planet radii be-tween 0.5 and 10 R ⊕ ). Of these, a total of 390 multiplanet systemswith 966 planets are included, with the remaining 1171 planets be-ing in single systems. The resulting population of KOIs is shownin Figure 2, where we plot histograms of the number of planetsper system ( N m , planet multiplicity), periods ( P ), period ratios ofapparently adjacent planet pairs ( P i + / P i ), transit durations ( t dur ),period-normalized transit duration ratios ( ξ ), transit depths ( δ ), andtransit depth ratios ( δ i + / δ i ). These distributions serve as the targetobserved distributions for our models. In particular, the distribu-tions of the period ratios, transit duration ratios, and transit depthratios are especially insightful to model because they probe the ar-chitectures of planetary systems, yet are insensitive to the stellarparameters. In order to compare our forward models to the
Kepler observa-tions, we need to find model parameters that result in simulatedcatalogues that are similar to the
Kepler
DR25 catalogue. Giventhe complexity and computational expense of the model, we take amultistage approach. First, we use an optimizer to identify good re-gions of parameter space. Results from the optimizer are then usedto train a Gaussian Process emulator (described in §2.6). Finally,we draw samples from prior distributions for model parameters anduse rejection sampling to construct the ABC posterior for inference.Here we describe the optimization stage that seeks the set ofmodel parameters that minimize the distance function (for given amodel, observed data set, and distance function). This is a challeng-ing problem for several reasons. First, evaluating the distance func-tion is computationally expensive, primarily due to the cataloguesimulation. Secondly, the parameter space is large. Our models in-volve several free parameters: 8, 10, and 11 for the non-clustered,clustered periods, and clustered periods and sizes models, respec-tively (even after fixing the break radius R p , break and minimumseparation for stability ∆ c ). There can be correlations or potentiallycomplicated interplay between model parameters. The third andmost challenging factor is that the forward model is stochastic dueto sampling variance and the finite number of targets. Even for fixedmodel parameters, the computed distance varies from one realiza-tion to the next due to Monte Carlo randomness in drawing targetproperties, physical properties of planetary systems, and simulated measurement noise. This means that traditional optimization algo-rithms that assume a deterministic function are not appropriate forour problem. In theory, one could reduce the variance in the dis-tances drawn from our forward model by simulating significantlymore targets than Kepler observed. While this could be a useful(but expensive) way to find the “best-fitting” model parameters, it isnot appropriate for accurately characterizing the uncertainties in themodel parameters. The summary statistics for the
Kepler catalogueinclude features that might be real or merely the result of small num-ber statistics. In ABC, our forward model should also have variancedue to the finite number of targets observed, in order for the ABCposterior to properly weight model parameters accounting for theextent of variations due to the finite number of targets. Therefore,we use the same number of targets as in the
Kepler catalogue duringboth the optimization stage (this section) and the emulator stage (asdescribed in the next section).We use an adaptive Differential Evolution optimizer func-tion with radius limited sampling, implemented by “BlackBoxOp-tim” ( https://github.com/robertfeldt/BlackBoxOptim.jl ). This package includes a general purpose optimization function“bboptimize()” which provides various algorithms, some of whichare designed to deal with stochastic noise. For our purposes, we usethe optimizer “adaptive_de_rand_1_bin_radiuslimited”. This opti-mizer uses a population-based algorithm that iteratively searches aspecified region in parameter space in order to try and minimize atarget fitness function without assuming that the function is differ-entiable. Thus it is well suited to high dimensional problems withstochastic noise. We use the total weighted distance given by Equa-tion (18) as the target fitness function, and leave most of the keymodel parameters of each model as free parameters with an allowedsearch range (minimum, maximum) for each parameter. Table 3 liststhe search ranges we used for each free parameter, in each model.We repeat 50 runs of the optimizer on each of our models,each with a different set of initial values for the free parametersdrawn randomly (uniformly, while for λ c and λ p , uniformly in log)within each of the search ranges for each parameter. For each run,we set the “population size” parameter equal to 4 n params , where n params is the number of free parameters in the model. We simulate N stars , sim = σ i , low (cid:54) σ i , high , we transform these twoparameters into two dummy variables r , r using a mapping fromthe unit square to a triangle (Osada et al. 2002). We also set acondition to avoid simulating the model if ln ( λ c ) + ln ( λ p ) > . Osada et al. (2002) prescribe the equation P = ( − √ r ) A + √ r ( − r ) B + r √ r C in their Equation (1) of Section 4.2 to uniformly map theunit square to any arbitrary triangle of vertices ( A , B , C ) , i.e. by drawing r and r randomly between 0 and 1. We adopt this transformation andset the vertices to A = ( ◦ , ◦ ) , B = ( ◦ , ◦ ) , and C = ( ◦ , ◦ ) for ( σ i , high , σ i , low ) .MNRAS , 1–31 (2019) He, Ford, and Ragozzine
Table 3.
Table of the free parameters and their search ranges (min, max)explored, of each model. The last column lists the parameter values of theclustered periods and sizes model used to generate the reference catalogueas described in §2.4.2. We set R p , break = R ⊕ and ∆ c = f σ i , high ( , ) ( , ) ( , ) λ c - ( . , ) ( . , ) λ p ( , ) ( . , ) ( . , ) α P (− , ) (− , ) (− , ) α R (− , ) (− , ) (− , ) − α R (− , ) (− , ) (− , ) − σ e ( , . ) ( , . ) ( , . ) σ i , high ( ◦ ) ( , ) ( , ) ( , ) σ i , low ( ◦ ) ( , σ i , high ) ( , σ i , high ) ( , σ i , high ) σ R - - ( , . ) σ P - ( , . ) ( , . ) Computing the full forward model as detailed in §2.2-§2.3 iscomputationally expensive. Generating a physical catalogue with N stars , sim = N stars , Kep = ∼ λ c , λ p , and thecluster width in log-period per planet σ P due to repeated draws forstability. (The procedure to simulate an observed catalogue froma pre-existing physical catalogue is faster, taking just a few sec-onds.) Thus, it would be prohibitively time-consuming to simulatethe full model for millions of iterations. Fortunately, for the purposeof finding the region(s) of parameter space that result in observedcatalogues similar to the Kepler catalogue, we only need to store thetotal weighted distance for each proposed set of model parametersand do not necessarily need to save all the information in the cata-logues. Based on our initial exploratory analyses, the total weighteddistance has a single dominant mode for each physical model, whichis identified by the optimizer from §2.5. In the vicinity of that mode,the mean of distance function varies smoothly, but with considerablevariance due to Monte Carlo noise from the finite number of stars.Thus, we can dramatically accelerate the exploration of parame-ter space by approximating the total weighted distance of our fullmodel using a fast emulator. We adopt the machinery of Gaussianprocesses (GPs) to train an emulator for our distance function anduse the GP to explore the model parameter space in a computation-ally feasible manner (see the textbook by Rasmussen & Williams2006 for an extensive guide to and discussion of GPs, and O’Hagan2004 for an introduction to the GP emulator approach).A Gaussian process emulator is a statistical model that aims tomimic a more complicated and expensive function by “emulating”the outputs of the expensive function given the same inputs. Acovariance function specifies the correlation between draws fromthe GP for any pair of input values. For any set of inputs (i.e., modelparameters), the GP emulator returns a Gaussian distribution for itsprediction of the model output (i.e., distance) that depends on theobserved values of the function at a set of training points. Whilethe detailed outputs of our physical (clustered and non-clustered)models are complex, we only require the GP emulator to providea good approximation to the total weighted distance in a region ofparameter space that results in simulated catalogues that are a good match to the
Kepler data. For each set of model parameters, the GPemulator returns a distribution, which approximates the mean andvariance of the distribution of distances that would be returned bythe full model (including the effects of the finite number of targets).This allows us to explore the parameter space very quickly in orderto efficiently estimate how often realizations with a given set ofmodel parameters would result in a weighted distance less than themaximum acceptable distance for our ABC posterior sample.Here we provide an overview of the GP emulator before pro-viding more specific details below. For the remainder of this paper,we let θ denote the free parameters of our physical models (e.g. λ p , α P , etc. as listed in Table 3, of our non-clustered and clusteredmodels) and φ denote the (hyper) parameters of the GP emulator.We have a distance function D W ( θ ) (this is either D W , KS ( θ ) or D W , AD (cid:48) ( θ ) ) that we evaluated at a large number of points θ dur-ing the optimization stage. As discussed below, we use a subset ofpoints that yielded low distances during the optimization stage astraining points for the GP emulator (see §2.6.2). For a given modeland set of training points, we find the “best-fitting” values of thehyperparameters φ b f which maximize the log likelihood for theGP. Then, we use the emulator with φ b f to predict the distance D W ( θ ) for a much larger number of points in the model parameterspace. We draw trial values of θ from a prior and reject draws thatresult in a predicted distance greater than a distance threshold. Thedistribution of the accepted points provides a sample from the ABCposterior (see §2.6.3). We compute credible regions for each modelparameter θ from the quantiles of the accepted points. GPs are especially useful in our context due to their flexibility inmodelling stochastic processes with intractable functional forms.This is exactly the case for our distance function.Mathematically, a GP is described by a prior mean function m ( x ) and a covariance (i.e. kernel) function k ( x , x (cid:48) ; φ ) : f ( x ) ∼ GP (cid:0) m ( x ) , k ( x , x (cid:48) ; φ ) (cid:1) , (25)where f ( x ) is the function we wish to model ( f ( x ) = D W ( θ ) forour purposes), and φ are the hyperparameters of the kernel.The prior mean function can be used to model an underlyingdeterministic process if one is known (such as the periodic motionof a star in a series of radial velocity measurements, e.g. in Rajpaulet al. 2015). For our problem, an underlying functional dependence(i.e. D W as a function of the model parameters θ ) is not known;thus we use a constant prior mean function. When evaluating theemulator near training points, the predicted values will be stronglyaffected by the training points and only minimally affected by theprior mean. The choice of our mean function will be important whenevaluating the emulator far from training points.In an ideal world, we would supply enough training points toadequately characterize the model behavior over the full parameterspace, so the emulator would return accurate predictions at anygiven point. However, this is not practical for the entirety of our8–11 dimensional parameter space, as the number of training pointsis limited to just a few thousand, due to the computational cost oftraining the GP emulator. Therefore, we select training points to be inthe vicinity of the best-fitting region as found during optimization.We verified that the distribution of distances returned by the GPemulator is accurate in the region of interest. In order to ensure thatthe GP emulator consistently returns values well above the distancethreshold in other regions, we set the prior GP mean to a largevalue, i.e., near the high end of the distance for the training points MNRAS , 1–31 (2019) lustered Planetary Systems (e.g., m ( x ) =
75 for emulating the distance function involving KSdistance terms, D W , KS , and m ( x ) =
250 for emulating the distancefunction involving AD distance terms, D W , AD (cid:48) ). Since the prior GPmean is significantly larger than the distance threshold to be used byABC, the emulator will almost certainly return emulated distanceswell above distance threshold, when it is given parameter valuesfar from the training points. For the emulator to have a reasonablechance of returning a distance that would be accepted, the modelparameters must be near a sizeable number of training points thatcause the mean of the prediction to drop well below the prior mean.For the covariance function, we choose the squared exponentialkernel with a separate length scale, λ i , for each dimension (i.e.model parameter θ i , with i = , , ..., d where we let d denote thenumber of dimensions/model parameters), k ( x , x (cid:48) ; φ ) = σ f exp (cid:34) − (cid:213) i ( x i − x i (cid:48) ) λ i (cid:35) , (26)where φ = ( σ f , λ , λ , ..., λ d ) are the hyperparameters of the kernel.The hyperparameter σ f controls the strength of correlation betweenpoints, and also serves as the standard deviation of the Gaussianprior (i.e. far away from any training data). Intuitively, the lengthscales λ i govern how far points can be from one another, in eachdimension, before they become uncorrelated with each other.We find the best values for the hyperparameters φ by attempt-ing to maximizing the log marginal likelihood:log p ( y | X , φ ) = − y T k − y y −
12 log | k y | − n π, (27)where y = D W ( θ ) are the function values at the training points, k y = k ( θ , θ (cid:48) ; φ ) + σ n I is the covariance matrix for y (and σ n = ˆ σ (D W ) are the estimated uncertainties at the training points), and n is the total number of training points ( X = θ ). The choice oftraining points is described in §2.6.2.In practice, optimizing a 8–11D function is computationallyexpensive. Rather than simultaneously optimizing all of the hyper-parameters, we set σ f = λ i (cid:48) (listed in Table 4) informed by inspection of the dis-tribution of training points. Then, we maximize the log-marginallikelihood by varying an overall length scale factor λ global , where λ i = λ global λ i (cid:48) . Note that for the purposes of training the emula-tor for our clustered models, we also transform the parameters λ c and λ p into a sum and difference of their log-values, ln ( λ c λ p ) andln ( λ p λ c ) , as the region of best values for these transformed parame-ters are more Gaussian than that of the rates of clusters and planetsper cluster parameters themselves. For each model, we choose a set of training points { θ } for theGP emulator by taking a subset of the model evaluations from theoptimization runs as described in §2.5. Since we ran 50 individualoptimizers with 5000 total model evaluations per run, we have a poolof 2 . × model evaluations scattered around the d -dimensionalparameter space. We rank order these points by D W and choose thetop 10 points, keeping every 10 th point for a total of 10 points.The choice of keeping every 10 th point instead of all points is due to(1) the computational limits in calculating and inverting the kernelmatrix k ( θ , θ (cid:48) ; φ ) , which scales as n and thus prevents us fromusing more than a few thousand training points, and (2) the desireto include some points far enough away from the minimum so thatthe GP emulator makes reasonable predictions for a wider region ofparameter space. Table 4.
Table of length scale hyperparameters λ i and bounds for comput-ing credible regions for the parameters of each model. The same sets ofhyperparameters and box bounds were used for both the analyses involvingthe KS and the AD distances. Trial sets of model parameters were drawnuniformly in each bounded interval.Parameter Non-clustered Clustered periods Clustered periodsand sizes λ i Bounds λ i Bounds λ i Bounds f σ i , high ( , . ) ( . , . ) ( . , . ) ln ( λ c λ p ) * 0.03 ( . , . ) ( , ) ( , . ) ln ( λ p λ c ) - - 1 ( , ) (− , ) α P (− . , . ) (− , . ) (− , ) α R (− , ) (− , ) (− , ) α R (− , − ) (− , − ) (− , − ) σ e ( , . ) ( , . ) ( , . ) σ i , high ( ◦ ) 30 ( , ) ( , ) ( , ) σ i , low ( ◦ ) 0.25 ( , . ) ( , . ) ( , ) σ R - - - - 0.15 ( . , . ) σ P - - 0.1 ( . , . ) ( . , . ) Note. *This is ln ( λ p ) for the non-clustered model. The combination of stochastic noise in D W due to the MonteCarlo noise of our simulations and keeping the best rank-orderedpoints would introduce a bias towards smaller than average distancesat these points. In order to avoid this bias, we recompute the values D W by simulating a new catalogue with the full SysSim model ateach of these points. We train the GP emulator using updated valuesof D W for a random subset of 2000 points of the 10 availablepoints. The prior mean function, kernel function, and training points fullydefine a GP model. For each model and distance function ( D W , KS or D W , AD (cid:48) ) combination we train a GP emulator and use it topredict the distance function at a large number of points θ in the d -dimensional model parameter space. In order to improve com-putational efficiency, we draw samples from a reduced range foreach parameter (based on inspecting the results of the optimiza-tion stage). Effectively, we assume a uniform prior for each modelparameter by drawing points uniformly in the d -dimensional box,which the bounds for each parameter are specified in Table 4.For each draw of model parameters from the prior distribution,we draw an emulated distance from the GP emulator. Finally, weonly accept draws if the emulated distance is less than a distancethreshold. We repeat this procedure until we have accumulated 10 draws from the ABC posterior.We adopt distance thresholds of D W , KS =
35 and D W , AD (cid:48) =
140 for both the clustered models, and a distance threshold of D W , KS =
55 for the non-clustered model. The distance thresholdsare less than the medians of the training points, so that the mean andvariance of the GP emulators are constrained both in and around the We do not train an emulator or compute credible regions for the non-clustered model using the AD distances, because upon inspection of theobserved catalogues generated using the best model parameters resultingfrom the optimization scheme in Section §2.5, we find that the overall rateof planets is such a poor fit with this distance function and model that it isnot worthy of more detailed investigation.MNRAS , 1–31 (2019) He, Ford, and Ragozzine perimeter of the regions of parameter space that are plausibly goodfits. Given the stochastic nature of our forward model, even a per-fect model would typically return a distance of ∼ × √ (cid:39)
20, i.e.,the number of component distance functions in our total weighteddistance times the square root of the ratio of the number of targetstars used when generating weights to the number of target starsused during the inference step). The best distances found during theoptimization stage ( ∼
20 for the clustered models and ∼
45 for thenon-clustered model, when using the KS distance) set the smallestdistance thresholds we could have chosen. The use of larger distancethresholds, means that credible regions based on our ABC poste-rior samples will be somewhat larger than the ideal posterior if wehad been able to use a smaller distance threshold. In practice, thesize of the ABC posterior credible intervals appear to be shrinkingonly slowly as we decrease the distance threshold further. Since thecredible regions can only shrink with a smaller distance threshold,the credible regions that we report are conservative, i.e., larger thanif we would have obtained with infinite computational resources.In Figure 3, we display a “corner” plot (plotted using corner.py ; Foreman-Mackey 2016) showing the ABC posteriorbased on a sample of 10 for the clustered periods and sizes modelwith the KS distance function. Similar figures for the non-clustered,clustered periods, and alternative MMR inclinations models areavailable in supplementary online material. Similarly, supplementalFigure 4 shows the analogous plot using the AD distance function.We interpret and discuss the meaning of these results in detail inthe next section. In Table 5, we report the best-fitting values of the free model param-eters for each of the models. We discuss the results and the effectof each of the free model parameters on the simulated observedpopulation in this section.
First, we briefly summarize how our models compare with eachother, before reporting the results for each parameter of each modelin detail. Our clustered models are clearly preferred over the non-clustered model, as evidenced by the significantly smaller best-fitting weighted distances as described in §2.6. The improvement inthe clustered models can be traced to improvements in the compo-nent distances for the multiplicities, period ratios, transit durationratios for planet pairs not near resonance, and to a lesser extent peri-ods and transit durations. Going from clustered periods to clusteredperiods and radii results in some component distance improving(based on period ratios and radius ratios), but the change in totaldistance is more modest.All three models are able to match the overall observed rateof planets (via our analysis involving KS distances; the analysis in-volving AD distances fails for the non-clustered model). The periodand radius (power/broken-power law) distributions are broadly con-sistent across all the models and distances, suggesting a shallowlyincreasing number of planets in log-period, a roughly flat distribu-tion of planet sizes below the break radius, and a sharply decreasingoccurrence of planets above the break radius. Both clustered modelsprefer: (1) a sizable fraction ( ∼ ∼ ◦ . These findings are consistent across analyses using either the KS and AD distances, corroborating previous stud-ies (e.g., Mulders et al. 2018). All three models imply that mostplanets have small orbital eccentricities. Finally, the clustered mod-els suggest that both periods and planet sizes are highly clustered,with consistent results for the cluster widths in log-period and radiususing both KS and AD distances.Our results are also relatively unchanged when consideringthe alternative MMR inclinations model, i.e. our fully clusteredmodel without the reduction of mutual inclinations for the near-MMR planets. All the parameters with the exception of f σ i , high remain consistent within statistical uncertainties. While f σ i , high issomewhat reduced ( ∼ . f σ i , high ∼
40% of highly mutually inclined systems includesplanets near an MMR that have low mutual inclinations. Overall,the rate of planets per star, the period and radius distributions, theeccentricity and mutual inclination scales, and the period and radiuscluster widths all remain the same. λ c , λ p ) In our two clustered models, the parameters λ c and λ p parameterizethe mean numbers of clusters and planets per cluster, respectively,before rejection-sampling and any truncation. The true mean rateof clusters in our simulations is somewhat less than λ c , while thetrue mean rate of planets per cluster is typically greater than λ p forsmall values, and always greater than one, due to the draws fromthe zero-truncated Poisson distribution. These parameters largelycontrol the overall rates of planets and how they are distributedwithin and between clusters.It is clear that with increasing λ c and λ p , the overall numberof observed planets increases, due to there being intrinsically moreplanets around each star. This effect is perhaps slightly more sensi-tive to λ c , likely because increasing λ p is more likely to result inrejected systems which slows the increase in the number of planetswith increasing λ p (since it is harder to fit more planets into the samecluster than to increase the number of clusters, due to stability). Theproduct of these two parameters is most directly constrained by thetotal rate of planets per star (i.e. the distance term D f in equation18). In our clustered periods and sizes model, we find that λ c (cid:39) . + . − . indicating that ∼
80% of systems with planets have onecluster. The number of planets per cluster peaks at four with λ p (cid:39) . + . − . . The combination of these gives a typical λ c λ p ∼ . σ credible regions overlap. The overall rate of planetsper star ( ∼ .
6) is similar.For the clustered periods model, we find similar values forthe rates of clusters and planets per cluster using both KS and ADdistances: λ c (cid:39) . + . − . and λ p (cid:39) . + . − . using KS distances.Again, while the KS analysis results in slightly more planets percluster and fewer clusters than the AD analysis, the difference isnot significant given the credible regions. These results are verysimilar to those of the clustered periods and sizes model using KSdistances.For our non-clustered model, the number of planets per system N p is described by a single parameter, N p ∼ Poisson ( λ p ) . Wefind that using KS distances, λ p (cid:39) . + . − . , implying a medianvalue of ∼ . MNRAS000
6) is similar.For the clustered periods model, we find similar values forthe rates of clusters and planets per cluster using both KS and ADdistances: λ c (cid:39) . + . − . and λ p (cid:39) . + . − . using KS distances.Again, while the KS analysis results in slightly more planets percluster and fewer clusters than the AD analysis, the difference isnot significant given the credible regions. These results are verysimilar to those of the clustered periods and sizes model using KSdistances.For our non-clustered model, the number of planets per system N p is described by a single parameter, N p ∼ Poisson ( λ p ) . Wefind that using KS distances, λ p (cid:39) . + . − . , implying a medianvalue of ∼ . MNRAS000 , 1–31 (2019) lustered Planetary Systems f i ,high = 0.42 +0.080.07 . . . . . l n ( c ) ln( c )=0.10 +0.440.35 . . . . . l n ( p ) ln( p )=1.35 +0.360.44 . . . . . P P = 0.40 +0.640.56 . . . . . R R =1.02 +0.640.70 . . . . . R R =4.41 +1.360.79 . . . . . e e =0.020 +0.0140.010 i , h i g h () i ,high ( )=47.74 +17.2816.74 . . . . . i , l o w () i ,low ( )=1.40 +0.540.39 . . . . . R R = 0.31 +0.070.07 . . . . . f i ,high . . . . . P . . . . . ln( c ) . . . . . ln( p ) . . . . . P . . . . . R . . . . . R .
01 0 .
02 0 .
03 0 .
04 0 . e
20 40 60 80 i ,high ( ) . . . . . i ,low ( ) .
16 0 .
24 0 .
32 0 .
40 0 . R .
12 0 .
16 0 .
20 0 .
24 0 . PP = 0.21 +0.040.04 Figure 3.
Marginal ABC posterior distributions of the model parameters for the clustered periods and sizes model, showing the projections of the points thatpass a distance threshold of D W , KS =
35 as drawn from the GP emulator. The prior mean function was set to a constant value of 75. is comparable to the value in our best-fitting clustered models. Asdiscussed before, we do not recover meaningful results for thismodel using AD distances.Returning to our fully clustered model, the values of λ c and λ p are somewhat shifted in our alternative MMR inclinations model.There are slightly more clusters per system and fewer planets withineach cluster, using both KS and AD distances, although these dif-ferences are within the statistical uncertainties. The overall rate ofplanets per system effectively remains the same.We conclude that the Kepler observations robustly constrainthe mean number of planets per system and the results are insensitiveto the choice of forward model or distance function. While there is more variance across models and distance functions for the inferredratio of planets per cluster to clusters per system, the differencesare less than the statistical uncertainties for three of the four casesconsidered (i.e., clustered periods and fully clustered models, eachwith KS and AD distances). α P ) For all our models, the overall period distribution is described by asingle power law between 3 and 300 d with index α P . We allowed α P to vary between − α P = − MNRAS , 1–31 (2019) He, Ford, and Ragozzine
Table 5.
Table of the best-fitting values for the free parameters of each model, with 1 σ credible intervals. The last two columns (under the label “Alt. MMRinclinations”) show the results for the clustered periods and sizes model without the lowering of mutual inclinations for near-MMR planets. The rows with λ c and λ p are equivalent to those with ln ( λ c ) and ln ( λ p ) , respectively, and are shown for interpretability. We set R p , break = R ⊕ and ∆ c = f σ i , high . + . − . . + . − . . + . − . . + . − . . + . − . . + . − . . + . − . ln ( λ c ) - - − . − . + . − . − . + . − . − . − . + . − . . + . − . . + . − . . + . − . λ c - - 0 . . + . − . . + . − . . . + . − . . + . − . . + . − . . + . − . ln ( λ p ) .
25 1 . + . − . .
44 1 . + . − . . + . − . .
31 1 . + . − . . + . − . . + . − . . + . − . λ p . . + . − . . . + . − . . + . − . . . + . − . . + . − . . + . − . . + . − . α P − . − . + . − . . + . − . . + . − . . + . − . . + . − . . + . − . − . + . − . α R − . − . + . − . − . − . + . − . − . + . − . − . − . + . − . − . + . − . − . + . − . − . + . − . α R − . − . + . − . − . − . + . − . − . + . − . − . − . + . − . − . + . − . − . + . − . − . + . − . σ e . + . − . . + . − . . + . − . . + . − . . + . − . . + . − . . + . − . σ i , high ( ◦ ) 60 62 + −
50 49 + − + −
50 48 + − + − + − + − σ i , low ( ◦ ) 0.3 0 . + . − . . + . − . . + . − . . + . − . . + . − . . + . − . . + . − . σ R - - - - - 0.3 0 . + . − . . + . − . . + . − . . + . − . σ P - - 0.2 0 . + . − . . + . − . . + . − . . + . − . . + . − . . + . − . three of our models result in an increasing occurrence of planets withlog-period, strongly disfavouring a flat distribution ( α P = − α P are consistent within statistical uncertainties acrossall of our clustered models. However, assuming a non-clusteredmodel results in both a shallower slope and significantly smalleruncertainties on α P .The clustered periods and sizes model results in α P (cid:39) . + . − . using KS distances, and favours a shallower slope whenusing AD distances, α P (cid:39) . + . − . . The results for the clusteredperiods models are very similar for the KS and AD analyses, giving α P (cid:39) . + . − . and α P (cid:39) . + . − . , respectively. The constrainton α P appears to be much tighter for the non-clustered model, with α P (cid:39) − . + . − . . We conclude that one should be cautious aboutoverinterpreting measurements of α P that have assumed orbital pe-riods are drawn independently for planets within the same system. α R , α R ) We allowed α R and α R , the power law indices for the brokenpower-law as given in Equation (6), to vary from − − R p , break control the overall distribution of planetaryradii; we fix R p , break = R ⊕ here, so α R controls the distributionof Earth-sized planets between 0.5 and 3 R ⊕ and α R controls thedistribution of larger planets between 3 and 10 R ⊕ . The continuityof p ( R p ) induces a correlation between α R and α R .In our distance function, the most direct constraint on the radiusdistribution comes from fitting the transit depth and transit depthratio distributions. The sizes of planets also affect the observed pe-riod ratio distribution, since small planets can be packed closer toeach other but may be harder to detect. Formally, the radius distri-bution must be inferred simultaneously with the period distributionand other model parameters (Youdin 2011). Nevertheless, it appears that one could infer the overall radius distribution reasonably welleven without simultaneously modelling the planetary architectures.There is a clear need for a broken power law for the radiusdistribution, given the joint posterior for α R and α R stronglyexcludes equal values. Moreover, the radius distribution falls muchmore quickly above R p , break than below it; α R (cid:39) − . + . − . while α R (cid:39) − . + . − . for our clustered periods and sizes model with KSdistances (and slightly steeper values using AD distances). Similarly,the clustered periods and non-clustered models also prefer α R (cid:46) − α R (cid:39) −
1. This is expected given that the
Kepler populationdoes not have a high occurrence of larger, Jupiter-sized planets. Thedifferences in the constraints on α R and α R across models anddistances functions are smaller than the statistical uncertainties.The 68.3% credible interval for α R includes a flat distributionof α R = − α R , we caution that our radiusmodel cannot include a local minima, i.e., the radius valley in theunderlying distribution (Fulton et al. 2017; Van Eylen et al. 2017;Hsu et al. 2019). Therefore, we leave the generalization of our modelto allow for a more flexible radius distribution for future studies. σ e ) The orbital eccentricity distribution primarily affects the transitduration and duration ratio distributions, with a slight influence onthe period ratio distribution due to the stability criteria. In all ourmodels, the eccentricities are very low; for the clustered models,both KS and AD analyses result in σ e (cid:39) . σ e ∼ . MNRAS000
Kepler populationdoes not have a high occurrence of larger, Jupiter-sized planets. Thedifferences in the constraints on α R and α R across models anddistances functions are smaller than the statistical uncertainties.The 68.3% credible interval for α R includes a flat distributionof α R = − α R , we caution that our radiusmodel cannot include a local minima, i.e., the radius valley in theunderlying distribution (Fulton et al. 2017; Van Eylen et al. 2017;Hsu et al. 2019). Therefore, we leave the generalization of our modelto allow for a more flexible radius distribution for future studies. σ e ) The orbital eccentricity distribution primarily affects the transitduration and duration ratio distributions, with a slight influence onthe period ratio distribution due to the stability criteria. In all ourmodels, the eccentricities are very low; for the clustered models,both KS and AD analyses result in σ e (cid:39) . σ e ∼ . MNRAS000 , 1–31 (2019) lustered Planetary Systems §4.1.1): this model is unable to produce enough observed planetswith small period ratios, and thus desires near circular orbits inorder to space planets as close as possible while remaining stable.As such, the fits to the transit duration and ξ distributions are alsoworse in the non-clustered model (see supplemental Figure 7).Our results are consistent with those of several previous stud-ies based on properties of systems with multiple transiting planets.Fabrycky et al. (2014) used an early sample of 899 planets in mul-titransiting Kepler systems and fit to the ξ distribution to constrainboth eccentricity and mutual inclination dispersions. They found abest fit of σ e (cid:39) .
032 (also assuming a Rayleigh distribution of e ), but caution that a wide range of σ e were plausible since the ξ distribution is not as sensitive to e as it is to mutual inclinations. Wu& Lithwick (2013), who found e ∼ .
01, and Hadden & Lithwick(2014), who reported an rms of σ e (cid:39) . + . − . , both performedanalyses of TTVs to arrive at these results. Van Eylen & Albrecht(2015) assumed a Rayleigh distribution and arrived at a slightlylarger value of σ e = . ± . σ e = . Kepler discoveries. They found a larger mean eccentricityin the range of 0.1–0.25 for their population, but cautioned that themeasurement uncertainty and potential systematic effects in the hoststellar densities could bias the derived eccentricities. This concernhas now been mitigated thanks to improved stellar properties and theability to select a cleaner sample of host stars thanks to Gaia DR2.Shabram et al. (2015) investigated the eccentricity distribution ofa small sample of mostly single, large, short-period planets usingboth transit and occultation ( h = e cos ω and k = e sin ω ) measure-ments. When they assumed a one-component Gaussian distributionfor h and k (equivalent to a Rayleigh distribution for e ), they found avalue of σ e = . + . − . . However, they also show that adopting atwo-component Gaussian (very similar to a mixture of two Rayleighdistributions) is favoured. That model results in about 90% of plan-ets with σ e = . + . − . and 10% with σ e = . + . − . .Some studies have found evidence suggesting that the eccen-tricity distribution may differ between systems with a single transit-ing planet and multiple transiting planets. Moorhead et al. (2011)found a difference in the transit duration distribution (normalizedby the duration estimated for a transit over the diameter of the starwith a circular orbit of the same period) that was statistically sig-nificant, but cautioned about the potential impact of uncertaintiesin the stellar parameters. Xie et al. (2016) also found a differencein the eccentricity distribution for Kepler single transiting planets( ¯ e ≈ .
3) and systems with multiple detected transiting planets( ¯ e = . + . − . ), based on a larger catalogue and stellar proper-ties updated by LAMOST. A similar finding is reported by VanEylen et al. (2019), who find σ e = . ± .
06 for singles and σ e = . + . − . for multis, using a smaller sample of planets butwith much more precisely measured stellar properties, thanks to as-teroseismology. Mills et al. (2019) updated stellar properties basedon Gaia and spectra from the California–Kepler Survey for a largesample of planets. They found that singles have higher mean eccen-tricities, with ¯ e = . + . − . for singles and ¯ e = . ± . (cid:39)
69% of singlesbeing drawn from a low-eccentricity population with σ e , low < . σ e , high > .
3. Formally, our low-inclination population contributesa small number of systems with a single transiting planet detected.However, we find that this fraction is so small that the transit durationdistribution of singles will remain a robust probe of the propertiesof the eccentricity distribution of the highly excited population.Several of these studies also find correlations of eccentricities withplanet radii (e.g., Hadden & Lithwick 2014; Shabram et al. 2015;Dawson, Lee & Chiang 2016) and stellar metallicity (e.g., Shabramet al. 2015; Mills et al. 2019).In this study, we adopted a single Rayleigh distribution of e for simplicity and to avoid increasing the number of model param-eters. Thus, at some level, our results are likely averaging over asmall population with larger eccentricities and/or any correlationsbetween eccentricity and planetary or stellar properties. However,we expect this to have a relatively small effect on our results for theeccentricity distribution, since the ξ distribution is only observablefor systems with multiple transiting planets. While the transit du-ration distribution is affected by planets regardless of multiplicity,this provides a weaker constraint than the ξ distribution as almosthalf of the known transiting planets are in systems with multipletransiting planets. We encourage future work to explore the impactof allowing the high-inclination population to have a broader dis-tribution of eccentricities and potentially correlations with stellarproperties. f σ i , high , σ i , high , σ i , low ) The mutual inclination angles between planets have pronouncedeffects on the multiplicity distribution and the period-normalizedtransit duration ratio, as discussed in §2.4.1. As listed in Table 3,we allow the fraction of systems with broad mutual inclinations f σ i , high to vary in the entire range between 0 and 1, and the mutualinclination scales σ i , high and σ i , low for the Rayleigh distributionsto vary between 0 ◦ (coplanar orbits) and 90 ◦ (effectively isotropic),with the constraint that σ i , low (cid:54) σ i , high . As discussed in §2.1.6, σ i , low also governs the mutual inclinations for planets near MMRs.For both of our clustered models, the fraction of systems with σ i , high is close to 40%; f σ i , high (cid:39) . + . − . (0 . + . − . ) using KS(AD) distances for the clustered periods and sizes model. Resultswith the clustered periods model are very similar. The broad mutualinclination scale σ i , high is clearly greater than ∼ ◦ , but poorly un-constrained at significantly larger inclinations, since these systemscontribute almost solely to the number of observed single-transitingsystems.Our result of f σ i , high ∼ . f iso (cid:39) .
38 reported in Mulders et al. (2018) and consistent withprevious studies on the
Kepler dichotomy (e.g. Ballard & Johnson2016 for M dwarfs). For the rest of the systems, the mutual incli-nations between planets are only a few degrees, σ i , low (cid:39) . + . − . (1 . + . − . ) degrees using KS (AD) distances. Additionally, (cid:39) (cid:39)
40% of systems labelled high mutualinclination, were actually assigned a low inclination to the sys-tem’s reference plane, due to being near a first-order MMR with anadjacent planet. While σ i , low is small, strictly coplanar orbits arestrongly prohibited primarily due to the log ξ distribution.We note that in the clustered periods and sizes model, if wedo not include the lowering of mutual inclinations for planets near MNRAS , 1–31 (2019) He, Ford, and Ragozzine an MMR, f σ i , high decreases to ∼ .
32. This is expected becausemore systems need to host multiple planets with relatively smallmutual inclinations in order to make up for the lost contributionto the multitransiting (2+ observed planets) systems from the near-MMR planets. Thus, this value is comparable to the true fraction ofsystems with highly mutually inclined planets in our usual clusteredmodels. Notably, the values of σ i , high and σ i , low do not change inthis model, highlighting their robustness to our treatment of thenear-MMR planets.The non-clustered model is unable to adequately match theobserved multiplicity distribution as shown in Figure 4 and dis-cussed later in §4.1.1. This leads the non-clustered model to favourvery small values of f σ i , high and σ i , low . Since only a few percent ofsystems have high inclinations, the value of σ i , high is only weaklyconstrained. This demonstrates the importance of adopting a clus-tered model, not just for inferring the distribution of orbital periods,but also for inference about the distribution of mutual inclinations.One reason for this is that the correlation between the true multi-plicity and inclination distribution noticed by earlier authors (e.g.,Lissauer et al. 2011b; Tremaine & Dong 2012) is partially brokenwhen including the period ratio distribution, which sets the typ-ical spacing of planets and affects the observed multiplicity in anon-clustered model.Our results for the mutual inclination distribution are in strongagreement with many previous studies. The first study on the mutualinclinations of the Kepler exoplanets was done by Lissauer et al.(2011b), who also simulated planetary systems in order to matchthe observed multiplicities of the first few months of
Kepler data.They found an excess of single transiting systems, and by fitting themultis alone with a mixture of 3- and 4-planet systems found a bestfit of σ i = ◦ using a Rayleigh distribution of mutual inclinations.Tremaine & Dong (2012) combined both Kepler and RV surveyresults to estimate mean inclinations in the range 0–5 ◦ ; a similarresult of (cid:54) ◦ was reported by Johansen et al. (2012), who tried tomatch the number of single, double, and triple transiting systemsby assuming all systems are intrinsic triples. Figueira et al. (2012)also used Kepler and RV (HARPS) data to constrain σ i (cid:54) ◦ as-suming a Rayleigh distribution, although they included only planetslarger than 2 R ⊕ and noted that larger inclinations of ∼ ◦ are pos-sible if the mass–radius relationship is more extreme. Fang & Mar-got (2012) explored combinations of bounded-uniform and zero-truncated Poisson distributions of multiplicities, with Rayleigh andRayleigh of Rayleigh distributions of mutual inclinations, in orderto fit the observed multiplicity and ξ distributions. They found thebest fits with the bounded-uniform and either Rayleigh or Rayleighof Rayleigh distributions with σ = ◦ (or σ σ = ◦ ). Fabrycky etal. (2014) also used the observed multiplicity and ξ distributions toconstrain mutual inclinations in the range 1–2 . ◦ . A similar rangeof 0.3–2 . ◦ (for mean mutual inclinations) was reported by Xie etal. (2016). Thus, our results from the clustered models corroboratethese previous studies and also serve as the tightest constraints onthe mutual inclination distribution to date, with consistent resultsvia two independent analyses involving KS and AD distances.Zhu et al. (2018) took a different approach to modelling themutual inclination distribution compared to these previous works,by parametrizing the mutual inclination dispersion as a functionof planet multiplicity k , with σ i , k ∝ k α . By matching to boththe observed transiting multiplicity and TTV-inferred multiplicitydistributions, they found a steep inverse relation, − (cid:46) α < − σ i , = . ◦ . The results of our clusteredmodels, considering only the low mutual inclination population, are generally consistent with their mutual inclination dispersion forhigh-multiplicity systems. σ P , σ R ) The parameters σ P and σ R refer to the scale factors for the clusterperiod and radius distributions, respectively. While σ P is effectivelythe width of the cluster in log-period, per planet in the cluster, σ R isthe fixed width of the cluster in log-radius regardless of the numberof planets (recall Equations 4 and 8). These parameters are largelyconstrained by the period ratio and transit depth (i.e. radius) ratiodistributions, respectively.We find that σ P (cid:39) . + . − . (0 . + . − . ) using KS (AD) dis-tances for our model with just clustered periods and very simi-lar values for the model that includes clustering in radii, where σ R (cid:39) . + . − . using KS (0 . + . − . using AD) distances. We findthat extremely small values of σ P (cid:46) .
05 are highly excluded,likely due to the effects of our stability criteria which causes someclusters to be discarded after many attempts to fit planets in and thussignificantly decreasing the overall rate of observed planets.These values suggest that the periods and radii of planets arehighly clustered, as smaller values imply a greater degree of cluster-ing (i.e. smaller cluster widths). For example, if we assume σ P (cid:39) . P c =
20 d, thisimplies that (cid:39)
68% of draws ( ∼ (cid:39)
68% of planetsdrawn for a cluster of super-Earth sized planets around 2 R ⊕ willhave radii between 1 .
48 and 2 . R ⊕ , assuming σ R (cid:39) . Kepler data
As outlined in §2, our forward modelling procedure produces asimulated observed catalogue of transiting exoplanets, includingmeasured planet multiplicities per system, observed properties foreach planet (period, transit duration, and transit depth), and theratios of period, transit depth, and transit duration for pairs of ap-parently adjacent planets. In this section, we display simulated ob-served catalogues from each of our three models, non-clustered,clustered periods, and clustered periods and sizes, in Figures 4–6 respectively. For each model, we generate one catalogue with N stars , sim = N stars , Kep = θ used to generate these catalogues are listedin Table 5. In order to compute the 16 and 84 percentiles for eachbin, we also generate 100 additional catalogues with model pa-rameters drawn from our ABC posterior and that pass our (KS)distance thresholds after generating an observed catalogue usingthe full forward model. Appendix Figure 7 shows how each of thesecatalogues compare to the Kepler
DR25 catalogue in terms of thetotal and individual (weighted and unweighted) distance terms. Wediscuss how each of our models compare to the
Kepler data below.
In Figure 4, we show the results of our non-clustered model (themodel parameters for the population plotted in black are listed inTable 5). The panels from top to bottom display the marginal distri-bution of each summary statistic while the left-hand and right-hand
MNRAS000
MNRAS000 , 1–31 (2019) lustered Planetary Systems F r a c t i o n Simulated16% and 84%Kepler C D F SimulatedKepler D f = 0.0008 CRPD = 0.016 P (days)10 F r a c t i o n SimulatedKepler P (days)0.00.51.0 C D F KS = 0.099 AD = 0.017 P i + 1 / P i F r a c t i o n P i + 1 / P i C D F KS = 0.245 AD = 0.107 F r a c t i o n C D F KS = 0.047 AD = 0.006 i + 1 / i F r a c t i o n i + 1 / i C D F KS = 0.138 AD = 0.035 t dur (mins)0.0000.0250.050 F r a c t i o n t dur (mins)0.00.51.0 C D F KS = 0.137 AD = 0.045 F r a c t i o n Near MMRNot near MMR C D F Near MMRNot near MMR KS = 0.108 AD = 0.017 KS = 0.235 AD = 0.099 Figure 4. Non-clustered model: a simulated observed catalogue of exoplanet systems generated from our non-clustered model.
Left-hand panels: histogramsof the summary statistics, with the same panels as the ones in Figure 2. The black hollow histograms show one simulated catalogue (model parameters listed inTable 5) while the
Kepler
DR25 exoplanets are plotted as grey shaded histograms for comparison. The red dashed histograms show the 16 and 84 percentilesof each bin based on 100 simulated catalogues with parameters drawn from our emulator with D W , KS < Right-hand panels: the corresponding CDFs tothe left-hand panels. The distances for each summary statistic are shown and labelled as arrows between the CDFs. The observed multiplicities are poorly fit,with too few multiplanet systems. While the period distribution is reasonably reproduced, the period ratio distribution is very poorly modelled. Similarly, thetransit depth ratio distribution of the model does not peak as sharply around unity as the data.MNRAS , 1–31 (2019) He, Ford, and Ragozzine sides show the empirical PDFs and CDFs, respectively. The
Ke-pler
DR25 population (as shown in Figure 2) is overplotted as greyhistograms for comparison. The lower right corners of each CDFpanel display the relevant distance terms between the model and the
Kepler population.The non-clustered model is able to produce a population ofsimulated observed exoplanets that resembles the
Kepler
DR25population in some regards. The overall rate of observed planetsper target can be matched almost exactly (the D f term is almostzero). The bulk of the period distribution appears well modelledby a simple power law with α P ∼ − . δ is slightly left-skewed (i.e., modeis at larger value than the median) whereas our model producesa symmetric distribution. The transit duration distribution is alsoreasonably well modelled, although we produce an excess of grazingtransits with zero duration. The distribution is also somewhat shiftedto longer durations due to the extremely small eccentricities. The ξ distribution of our simulated catalogue appears to be slightlymore skewed towards values log ξ >
0, suggesting that the mutualinclinations of the
Kepler planets are larger than σ i , low ∼ . ◦ .However, there are several major, clear differences betweenthe non-clustered model and the Kepler planet catalogues. First,the planet multiplicity distribution is poorly modelled and producesfar too few multiplanet systems. This shortcoming of the modelis persistent for many realizations near the best-fitting parameters.Secondly, the period ratio distribution is clearly not well capturedby this model, motivating our use of clustered models. While mostof the
Kepler adjacent planet pairs have small period ratios andthe distribution falls rapidly above
P ∼
3, the non-clustered modelproduces a much more gradual decline in the tail of the distribution.A stability criteria of ∆ c = Kepler radii ra-tios are highly peaked around unity. In summary, the non-clusteredmodel described in the previous section performed well in somerespects, but there are at least three obvious differences betweenthe simulated observed and
Kepler
DR25 populations: the observedmultiplicity, period ratio, and transit depth ratio distributions. Theseshortcomings motivate our clustered models.
Before we consider our full clustered periods and sizes model, wefirst explore the effect of adding just a clustering of orbital periods,i.e., our clustered periods model (the planetary radii are drawn inthe same manner as in our non-clustered model). In Figure 5, weshow the results of this model with the
Kepler
DR25 population(again, the model parameters used are listed in Table 5).Our clustered periods model is a significantly better descrip-tion of the
Kepler data than our non-clustered model. Notably, theobserved multiplicity and period ratio distributions bear a muchcloser resemblance to the data, and have significantly reduced dis-tances (CRPD distance for the observed multiplicity counts; KSand AD distances for the period ratio distribution). This is perhapsnot surprising given that the main difference between this modeland the previous model is the introduction of period clusters. The clustering in periods allows for some systems to contain planets thatare closely packed more often than in the case of independentlydrawn periods, while also allowing for a more gradually falling tailto large orbital period ratios. The clustering in periods also providea more flexible model for the numbers of planets per system sincethere are two parameters ( λ c and λ p ) instead of just one ( λ p ofthe non-clustered model). The fits to the period and transit durationdistributions are also slightly improved, while the transit durationratio distribution (for planets not near resonances; { ξ non − res } ) issignificantly improved, in both KS and AD distances. This sug-gests that the distribution of orbital eccentricities is not as low aswhat is implied by the non-clustered model, and is instead welldescribed by a Rayleigh distribution with σ e (cid:39) .
01. These resultsalso imply that the mutual inclination distribution is well describedby σ i , low (cid:39) . ◦ for ∼
60% of systems and planets near resonances.This also appears to produce noticeable peaks near the first-orderMMRs in the period ratio distribution (most apparent near
P (cid:39) . P (cid:39) ξ res dis-tribution. There is no significant change in the fit to the transit depthdistribution. The transit depth ratio distribution, on the other hand,is modelled just as poorly as and perhaps even worse than before.These results show that clustered periods alone cannot reproducethe highly peaked nature of the planet radii ratios observed in the Kepler multiplanet systems.The clustering in orbital periods is able to substantially improvethe modelling of the period ratio distribution, but the transit depthratio distribution remains poorly modelled by both the non-clusteredand clustered periods models. The transit depth ratios are morehighly peaked around one in the actual
Kepler population than inthe simulated catalogues of these models (see Figures 4 and 5),motivating our next investigation of the clustered periods and sizesmodel.
Several previous studies have found that the radii of exoplanetsobserved by
Kepler around a single target star are correlated, causingplanets to be more similar in size within each system compared tobetween systems (Ciardi et al. 2013; Millholland, Wang, & Laughlin2017; Weiss et al. 2018a). However, most recently these results havebeen called into question by Zhu (2019), who show that the clusteredradii can be reproduced by a re-sampling of signal-to-noise ratiosand conclude that these correlations are largely due to detectionbiases. Thus, we extend the clustering point process to includeclustering of the planetary radii (and implicitly planet masses). Thisis the full implementation of our clustered periods and sizes modelas detailed in §2.2.We plot the results of this model in Figure 6 and list the modelparameters in Table 5. As is the case in our previous clustered peri-ods model, the distances for the observed multiplicity, period, andperiod ratio distributions are small, and are significantly improvedcompared to the non-clustered model. There are no significant dif-ferences in the distances (KS and AD) for these observables betweenthe two clustered models. A similar conclusion can be drawn forthe transit durations and transit duration ratios; while both clus-tering in periods and clustering in both periods and sizes fit thesemarginal distributions equally well, both clustered models providesubstantially better descriptions of { t dur } and { ξ non − res } and per-haps slightly better fits to { ξ res } as compared to the non-clusteredmodel. However, the transit depth distribution appears to be mod-elled slightly worse than in the previous models. Also, while the MNRAS000
Kepler around a single target star are correlated, causingplanets to be more similar in size within each system compared tobetween systems (Ciardi et al. 2013; Millholland, Wang, & Laughlin2017; Weiss et al. 2018a). However, most recently these results havebeen called into question by Zhu (2019), who show that the clusteredradii can be reproduced by a re-sampling of signal-to-noise ratiosand conclude that these correlations are largely due to detectionbiases. Thus, we extend the clustering point process to includeclustering of the planetary radii (and implicitly planet masses). Thisis the full implementation of our clustered periods and sizes modelas detailed in §2.2.We plot the results of this model in Figure 6 and list the modelparameters in Table 5. As is the case in our previous clustered peri-ods model, the distances for the observed multiplicity, period, andperiod ratio distributions are small, and are significantly improvedcompared to the non-clustered model. There are no significant dif-ferences in the distances (KS and AD) for these observables betweenthe two clustered models. A similar conclusion can be drawn forthe transit durations and transit duration ratios; while both clus-tering in periods and clustering in both periods and sizes fit thesemarginal distributions equally well, both clustered models providesubstantially better descriptions of { t dur } and { ξ non − res } and per-haps slightly better fits to { ξ res } as compared to the non-clusteredmodel. However, the transit depth distribution appears to be mod-elled slightly worse than in the previous models. Also, while the MNRAS000 , 1–31 (2019) lustered Planetary Systems F r a c t i o n Simulated16% and 84%Kepler C D F SimulatedKepler D f = 0.0011 CRPD = 0.001 P (days)10 F r a c t i o n SimulatedKepler P (days)0.00.51.0 C D F KS = 0.056 AD = 0.008 P i + 1 / P i F r a c t i o n P i + 1 / P i C D F KS = 0.096 AD = 0.011 F r a c t i o n C D F KS = 0.042 AD = 0.005 i + 1 / i F r a c t i o n i + 1 / i C D F KS = 0.159 AD = 0.038 t dur (mins)0.0000.0250.050 F r a c t i o n t dur (mins)0.00.51.0 C D F KS = 0.067 AD = 0.008 F r a c t i o n Near MMRNot near MMR C D F Near MMRNot near MMR KS = 0.091 AD = 0.011 KS = 0.059 AD = 0.006 Figure 5. Clustered periods model: a simulated observed population of exoplanet systems generated from our clustered periods model. The panels are thesame as the ones in Figure 4. The black hollow histograms show one simulated catalogue (model parameters listed in Table 5) while the
Kepler
DR25 exoplanetsare plotted as grey shaded histograms for comparison. The red dashed histograms show the 16 and 84 percentiles of each bin based on 100 simulated catalogueswith parameters drawn from our emulator with D W , KS <
35. The observed multiplicity and period ratio distributions are remarkably improved compared tothose of the non-clustered model. However, the transit depth ratio distribution looks similar as before, since no clustering in sizes is present in this model, andthus still poorly fits the data.MNRAS , 1–31 (2019) He, Ford, and Ragozzine F r a c t i o n Simulated16% and 84%Kepler C D F SimulatedKepler D f = 0.0004 CRPD = 0.001 P (days)10 F r a c t i o n SimulatedKepler P (days)0.00.51.0 C D F KS = 0.037 AD = 0.003 P i + 1 / P i F r a c t i o n P i + 1 / P i C D F KS = 0.077 AD = 0.015 F r a c t i o n C D F KS = 0.067 AD = 0.006 i + 1 / i F r a c t i o n i + 1 / i C D F KS = 0.15 AD = 0.044 t dur (mins)0.0000.0250.050 F r a c t i o n t dur (mins)0.00.51.0 C D F KS = 0.087 AD = 0.014 F r a c t i o n Near MMRNot near MMR C D F Near MMRNot near MMR KS = 0.125 AD = 0.025 KS = 0.05 AD = 0.005 Figure 6. Clustered periods and sizes model: a simulated observed population of exoplanet systems generated from our fully clustered model. The panelsare the same as the ones in Figures 4 and 5. The black hollow histograms show one simulated catalogue (model parameters listed in Table 5) while the
Kepler
DR25 exoplanets are plotted as grey shaded histograms for comparison. The red dashed histograms show the 16 and 84 percentiles of each bin based on 100simulated catalogues with parameters drawn from our emulator with D W , KS <
35. This model produces a great fit to the observed summary statistics of thedata. While the transit depth ratio distribution is improved and peaks more closely around unity due to the clustering in planet sizes, the distances ares notreduced compared to that of the non-clustered model, suggesting that there are still differences yet to be modelled. In particular, the observed
Kepler transitdepth ratio distribution appears slightly more asymmetric than what our models can produce. MNRAS000
Kepler transitdepth ratio distribution appears slightly more asymmetric than what our models can produce. MNRAS000 , 1–31 (2019) lustered Planetary Systems transit depth ratios qualitatively appear to be better fit with thismodel, there is almost no improvement in the KS or AD distance.A closer examination reveals that there is a slight offset in theCDFs, due to the observed distribution of transit depth ratios of ad-jacent Kepler planet pairs being asymmetric (in log). This suggeststhat while models with non-clustered planet sizes fail to predictthe highly peaked nature of the transit depth ratio distribution, amodel with clustered sizes still does not adequately reproduce thisproperty; there are additional features shaping the distribution ofadjacent-planet radii ratios that require a more complex model toexplain. In particular, modelling the
Kepler depth and depth ratiodistributions would likely be improved by allowing for a valley inthe planet radius distribution (perhaps due to photoevaporation, coreheating, or some other process, see §1) that is not included in ourmodel. We discuss our speculations for how future models can begeneralized to better match these features in §4.6. physicalcatalogues ) between the models
A major benefit of forward modelling the
Kepler mission is thatwe can directly analyze the predictions of our models for the true,underlying exoplanetary systems. While the model parameters de-scribe the underlying system properties including the rate of planetsand clusters per system, the period and radius distributions, andthe orbital architectures, the physical catalogues generated by ourmodels can be directly examined to see Monte Carlo realizations ofpopulations of planetary systems. In Figures 7 and 8, we plot theunderlying distributions of planetary systems as generated by oneinstance of each model and 100 realizations of the clustered peri-ods and sizes model, respectively. The physical catalogues shownhere are the same simulated populations generating the observedcatalogues as shown in Figures 4–6. We make the following obser-vations: • Planets and clusters per system:
There is a huge differencein the way planets are distributed across different systems betweenour clustered and non-clustered models. The non-clustered model(dotted red) produces very few ( (cid:46) N p ∼ Poisson ( λ p ) . In this model, the most common planetarysystem consists of three planets (between 3 and 300 d and 0.5–10 R ⊕ ). The clustered models (dashed green and solid blue) producemany more systems with no planets due to draws of zero-clustersystems from N c ∼ Poisson ( λ c ) (since the number of planets percluster is a zero-truncated Poisson distribution, N p ∼ ZTP ( λ p ) in the clustered models, and λ c tends to be small). For systemswith planets, the multiplicity tends to be higher; the most commonplanetary system is a four-planet system. While there is a clearpreference for clustered planetary systems, the number of actualclusters tends to be small. For planet-harbouring stars, ∼
80% ofsystems have just a single cluster consisting of three or four planets(over the range of periods considered, 3–300 d).Thus, assuming that planetary systems are highly clustered, ourmodel predicts a large fraction of stars with no planets larger than 0.5 R ⊕ with periods between 3 and 300 d, while the remaining systemsto have many planets in this range. We discuss the occurrence ratesof planets in more detail in §4.3. • Orbital period distribution:
All three models exclude flator falling power laws in log-period (i.e, α P (cid:54) − P max =
300 d, likely due to edge effectsof our rejection-sampling algorithm; this pile-up is most severe inour non-clustered model. However, this has a minimal effect on ourobserved catalogues, as the probability of transits decreases withlonger periods. Additionally, the period distributions produced byour clustered models appear to be “rounded” and deviate from apower law, due to the draws of period clusters. • Spacings of adjacent planets:
Given the same stability criteriaof ∆ (cid:62) ∆ c =
8, the clustered models produce slightly narrowerseparation ( ∆ ) and underlying period ratio distributions, suggestingthat planetary systems are more tightly spaced than one would inferfrom a model with periods drawn independently. All three modelsproduce distributions of separations in mutual Hill radii that aresharply truncated at ∆ c , suggesting that many period or period scaledraws are attempted more than once. The distribution of ∆ beginsto fall out around 20, which is similar to the findings of Weisset al. (2018a) although they assumed a much simpler mass–radiusrelationship (Weiss & Marcy 2014) and only analysed the spacingsof observed planets. The underlying period ratio distributions arehighly peaked around ∼ . • Planet radius distribution:
As listed in Table 5, we chose α R = − , − .
8, and − R p , break = R ⊕ ) and sharply falls above it forthis non-clustered model. However, the catalogue generated fromthe clustered periods and sizes model exhibits a rounded distribu-tion for small radii and the break radius is not as clear (despite alsohaving set α R = − α R are consistent with a flat distribution for all three models. • Planet radius ratios:
The clustered periods and sizes (solidblue) model produces a radius ratio distribution strongly peakedaround unity. The effect of this peak increases for decreasing valuesof σ R due to more highly clustered planet radii. The radius ratiodistributions of the other two models reflect what the distributionlooks like if planet pairs with radii drawn independently from abroken power law are randomly paired. Since we also fit our simulated results to match the overall rate ofobserved planets per surveyed (FGK) star in the
Kepler mission,we can directly estimate the true fraction of stars with planets andthe mean rate of planets per star, given various ranges of periodand planetary radius. We use the same 100 simulated cataloguesgenerated for Figure 8 to calculate the following planet occurrencerates (and their 16% and 84% quantiles) from our clustered periodsand sizes model.
Planets between 0.5 and 10 R ⊕ : The fraction of stars with plan-ets (FSWP) with 3 d < P <
300 d and 0 . R ⊕ < R p < R ⊕ is0 . + . − . . Roughly half of all FGK stars have at least one planet inthis range based on our clustered periods and sizes model. In con-trast, our non-clustered model gives a remarkably different result,producing just a few percent ( ∼ MNRAS , 1–31 (2019) He, Ford, and Ragozzine F r a c t i o n Clustered P+RClustered PNon-clustered N c F r a c t i o n N p F r a c t i o n F r a c t i o n P (days)0.000.020.04 F r a c t i o n P i + 1 / P i F r a c t i o n R p ( R )0.000.010.020.03 F r a c t i o n R p , i + 1 / R p , i F r a c t i o n Figure 7.
Simulated physical catalogues generated by each of our three models. These physical populations correspond to the same populations used togenerate the observed catalogues shown as black histograms in Figures 4–6. From left to right, top to bottom: histograms of the total planet multiplicities,numbers of clusters per system ( N c ), numbers of planets per cluster ( N p ), separations in mutual Hill radii ( ∆ ), periods ( P ), period ratios ( P ), true planetradii ( R p ), and planet radii ratios ( R p , i + / R p , i ). The non-clustered model (dotted red lines) produces fewer stars with no planets due to the single Poissondistribution describing the number of planets N p , while the clustered models (dashed green and solid blue lines) produce many more such systems due todraws of zero-cluster systems from N c ∼ Poisson ( λ c ) . All three models produce shallowly rising period power laws, with a slight pile-up near P max = ∆ (cid:62) ∆ c =
8, the clustered models produce narrower ∆ and period ratio distributions, suggesting that planetary systems are tightly spaced. All three models shown here exhibit a radius power law that is relativelyflat for small sizes towards the break radius ( R p , break = R ⊕ ) and sharply falls above it. Only the clustered periods and sizes (solid blue) model produces anintrinsic radius ratio distribution strongly peaked around unity; the radius ratio distributions of the other two models reflect what happens if planet pairs withradii drawn from a broken power law are randomly paired. ilar result was drawn by Weissbein, Steinberg, & Sari (2012), whostudied the intrinsic multiplicity and mutual inclination distributionsand concluded that either planet occurrence is correlated betweenplanets in the same system and/or some stars are significantly moreplanet-rich than others. Returning to our clustered periods and sizesmodel, the mean number of planets per star is 2 . + . − . . This rateis somewhat less than the credible region values of λ c λ p , due to therejection-sampling procedure. Excluding the stars with no planets(in this range of periods and sizes), the mean number of such planetsper system rises to 4 . + . − . .Zhu et al. (2018) find that the fraction of Sun-like stars withat least one “ Kepler -like” planet is 30 ±
3% based on multiplic-ity and rate of TTVs, although they define “
Kepler -like” to onlyinclude planets larger than 1 R ⊕ (and periods less than 400 d). Incontrast, Zink, Christiansen, & Hansen (2019) find that 0 . + . − . of stars have at least one planet (with orbital periods between 0.5and 500 d and radii between 0.5 and 16 R ⊕ ), while modelling just the observed periods, radii, and multiplicity. Zink, Christiansen, &Hansen (2019) report 8 . ± .
31 planets per star hosting a planetarysystem (but excluding stars that host a single transiting planet with R p (cid:62) . R ⊕ ). This implies that most stars have at least seven plan-ets over the range they consider and only 5.5% of planets have oneto four such planets.Our finding of 0 . + . − . FGK dwarf stars hosting at least oneplanet is greater than that suggested by Zhu et al. (2018) and lessthan that suggested by Zink, Christiansen, & Hansen (2019). Whileeach study considers a slightly different range of orbital periods andsizes, we attribute the bulk of the differences relative to our modelas due to our use of a clustered model for the distribution of orbitalperiods within a planetary system. When we applied a non-clusteredmodel, we find ∼
97% of stars host at least one planet. Switching toa model which accounts for clustering in orbital period and planetsize (or just orbital periods) dramatically increases the fraction ofstars with no planets in our model. Similarly the mean number of
MNRAS000
MNRAS000 , 1–31 (2019) lustered Planetary Systems F r a c t i o n Simulated catalog16% and 84% N c F r a c t i o n +0.110.20 +0.130.10 +0.070.01 N p F r a c t i o n F r a c t i o n P (days)0.000.010.02 F r a c t i o n P i + 1 / P i F r a c t i o n R p ( R )0.000.010.02 F r a c t i o n R p , i + 1 / R p , i F r a c t i o n Figure 8.
Simulated physical catalogues generated from our clustered periods and sizes model. The panels are the same as the ones in Figure 7. The solidblack histograms correspond to the same population used to generate the observed catalogue in Figure 6, while the dashed red histograms correspond to the16 and 84 percentiles of each bin based on 100 simulated catalogues with parameters drawn from our emulator with D W , KS < planets per star with at least one planet (within the range of periodsand sizes considered) increases from 3 . + . − . to 4 . + . − . whenswitching from a non-clustered to a clustered model (see Figure 8).This demonstrates the importance of allowing for clustering whenperforming inference on the fraction of stars with planets or meannumber of planet per stars with at least one planet.Another important way in which our model differs from thoseof Zhu et al. (2018) and Zink, Christiansen, & Hansen (2019) isthat we allow for the mutual inclination distribution to have botha low and high inclination component. Since we treat the fractionof systems from the high inclination component as a free parame-ter, our simulations could have resulted in f σ i , high consistent withzero. Instead, values of zero were strongly excluded in both clus-tered models, regardless of which distance function was chosen.Interestingly, the simulations using a non-clustered model did re-sult in small values of f σ i , high and did not exclude f σ i , high =
0. Thisfinding also underscores the importance of allowing for clusteringwhen performing inference on the mutual inclination distributionof planetary systems.
Earth-sized planets:
We provide results for the rate of Earth-sized planets (here defined to be 0.75–1.25 R ⊕ ) around FGK stars,using our clustered periods and sizes model. The fraction of starswith at least one Earth-sized planet (between 3 −
300 days) is0 . + . − . . The mean number of these planets per star is 0 . + . − . . Considering systems with at least one planet (with sizes between0.5 and 10 R ⊕ and periods between 3 and 300 d), the probabilitythat a planetary system contains an Earth-sized planet (in the sameperiod range) is 0 . ± .
14. If we focus on planetary systems with atleast one planet (not necessarily Earth-size), then the mean numberof Earth-sized planets per system is 1 . + . − . (both with periods3–300 d). Broadly, about half of stars have an Earth-sized planetand most inner planetary systems have Earth-sized planets. Due tothe increasing frequency of planets at longer periods, most of theseEarth-like planets are in 100–300 d orbits. These conclusions followdirectly from combining the inferred distributions for orbital periodand planet radius (that are consistent with previous results) with theinferred fraction of stars with planetary systems. Kepler dichotomy
Previous studies have shown that many simple parametrized modelswith a single population significantly underpredict the number ofsystems with only one transiting planet. Since this implies two pop-ulations, this result is known as the
Kepler dichotomy. This paperfocuses on models that allow for two populations for the mutualinclinations, each with similar distributions for the remaining prop-erties of the system (i.e., no explicit dichotomy in the intrinsic distri-butions of multiplicity, orbital period, planet size, or eccentricity).
MNRAS , 1–31 (2019) He, Ford, and Ragozzine
A model with a single inclination population is nested inside ourmodel. Therefore, we can examine whether our posterior distribu-tion for the fraction of stars with a high mutual inclination planetaryhas significant weight near zero. For each of the clustered modelsand distance functions we considered, we find that the predictedfraction of high-inclination systems (40 ± Kepler dichotomy.Another proposed solution to the
Kepler dichotomy is an al-ternative model consisting of one population of high-multiplicitysystems (i.e., very similar to our low-inclination population) and asecond population of systems with a single planet (both within 3d < P <
300 d and 0 . R ⊕ < R p < R ⊕ ) . As shown in §4.3,the average number of planets per planetary system (i.e., star withat least one planet) is (cid:104) NPPS (cid:105) ∼ . ∼ .
56 in our fully clustered model.Therefore, the average star coming from the high-inclination pop-ulation in our model actually hosts multiple detectable transitingplanets, greatly increasing the fraction of viewing angles for whichthe star would be perceived as hosting a single transiting planet.Using values from our fully clustered model, simply replacingall the stars with at least one planet ( (cid:39) .
56) and assigned to the highmutual inclination population ( f σ i , high = .
42) with stars hosting asingle planet, would require FSWP × f σ i , high × (cid:104) NPPS (cid:105) (cid:39)
Kepler dichotomy, since the onlyway to fit enough systems with one transiting planet around starsthat are not already spoken for is to have more than one (highlyinclined) planet per star.Before completely discarding an abundance of single-planetsystems as an explanation for the
Kepler dichotomy, it is worth con-sidering this argument in more detail. In the fully clustered model, (cid:39)
29% of planets around stars nominally in the high-inclinationpopulation were actually assigned a low inclination (relative to thesystem mid-plane), due to the planet being near a first-order MMRwith another planet in the same system. This results in a subset ofthe planet pairs from the high-inclination population contributingto the number of systems observed to have two transiting planetsand detracting from the fraction of viewing geometries that the sys-tem would be observed as a single planet system. We consider thepossibility of replacing the planets with truly high-inclinations witha population of intrinsic single-planet systems by considering thefraction of stars hosting each of three types of planetary systems:(1) the fraction of stars with planetary systems (i.e. with at leastone planet) initially assigned to the low-inclination population isFSWP × ( − f σ i , high ) (cid:39) × f σ i , high × . (cid:39) × ( − . ) (cid:39)
71% of starsin order to fully replace the planets that remain assigned to highmutual inclinations. Thus, in the fully clustered model, the sum ofthe three populations comes to (cid:39) While our range of periods and radii includes many Hot Jupiters whichare thought to be relatively isolated, the fraction of stars with such systemsis known to be quite low ( ∼ Kepler dichotomy. low inclinations is to reduce, but not eliminate the tension in havingenough stars to accommodate the three populations. Repeating thiscalculation over 100 realizations of our fully clustered model resultsin the total number of stars required to be 110% ± f σ i , high ∼ .
31, as discussed in §3.6),but none of the planets in these systems are reassigned to a lowermutual inclination. The fraction of stars with planets is higher,FSWP ∼ .
67. The mean number of planets per planetary systemis similar, (cid:104)
NPPS (cid:105) ∼ .
2, meaning that we would require FSWP × f σ i , high × (cid:104) NPPS (cid:105) (cid:39)
87% of stars to host a single planet. Again, thisdoes not leave enough stars to host the planetary systems from thelow mutual-inclination population, for which we require FSWP ×( − f σ i , high ) (cid:39)
46% of stars. Repeating this calculation results inthe same conclusion for 73 of 100 realizations.Therefore, while explaining the
Kepler dichotomy with a pop-ulation of single planet systems cannot be strictly excluded by ourwork, it is disfavoured and strongly constrained. Even in the casewhere we cannot formally exclude it, the fraction of stars with plan-ets is quite high, requiring extremely efficient planet formation.Since our model does not include planets with P < R p < . R ⊕ , including stars with such planets could further in-crease the number of stars required to host multiple planet systems.Thus, explaining the Kepler dichotomy with a population of excesssingles leaves very little, if any, room for stars harbouring no planets.Similarly, we could consider an alternative model in which theexcess of stars observed as a single transit system is due to a pop-ulation of highly excited two-planet systems. In this scenario, thefraction of stars hosting two highly inclined planets would need tobe (cid:39) × f σ i , high (cid:39)
24% of stars inour standard fully clustered model, there are sufficient stars that the
Kepler dichotomy could be explained by population of highly in-clined systems with two (or more) planets. However, in order to haveenough stars for such model, we had to relax the assumption thatthe number of cluster per stars is drawn from a Poisson distribution.Instead, we must allow for the number of zero planet systems to bedecoupled from λ c , the rate of clusters per star, which is inferredfrom the observed multiple planet systems. Using our clusteredperiods model leads to similar conclusions. Of course, these con-clusions must be interpreted in the context of the parametrizationof our model, e.g., the use of Poisson and Rayleigh distributions.Altogether, our results suggest that the Kepler dichotomy is mosteasily explained by two populations with different inclination dis-tributions, as opposed to different multiplicity distributions.Zink, Christiansen, & Hansen (2019) proposed that the
Kepler dichotomy could be explained by a combination of geometric transitprobability, a distribution of mutual inclinations, and a detectionefficiency model that accounts for the order of planet detection bythe
Kepler pipeline. For their star and planet sample, they report thatthe observed ratio of double to single transit systems is 4% largerthan that predicted by their model when ignoring how the detectionefficiency depends on the order of detection. When using their modelfor how the detection efficiency differs for subsequent planets, thisdifference is reduced to 2%. However, Zink, Christiansen, & Hansen(2019) did not attempt to make use of information contained inthe distributions of orbital period ratios, transit depth ratios, or
MNRAS000
MNRAS000 , 1–31 (2019) lustered Planetary Systems transit duration ratios of planets in a single planetary system, as wasdone in this study. We find that these provide valuable informationfor characterizing the distribution of planetary architectures andevidence for the Kepler dichotomy.It is also useful to compare the planet detection efficiencymodels of the two studies. Zink, Christiansen, & Hansen (2019) fitseparate detection efficiency curves for the first planet (technicallythe first threshold crossing event, TCE) to be detected by the
Kepler pipeline and for subsequent planets (or other TCEs). For planets withorbital periods less than 200 d (i.e., the vast majority of detectionsin our sample), they find that the expected MES needed for a 50%probability of detecting the planet increases from (cid:39) . (cid:39) . Kepler data provide evidence for a high-inclination popula-tion of planetary systems to explain the observed
Kepler dichotomy.
Planets near a mean-motion resonance (MMR):
We find that ourclustered models can produce the spikes in the period ratio distri-bution near first-order MMRs similar to those observed by
Kepler by assigning planets from the high inclination population to thelow inclination population when near an MMR. Since we do notexplicitly consider models with an explicit excess near MMRs, wedo not exclude that possibility. Nevertheless, our results do provideinsight into how significant an excess would be needed to createsuch spikes.While f σ i , high describes the fraction of systems initially as-signed to the high-inclination population, a fraction of the planetsin these systems are still assigned a low mutual inclinations if theyare near a (first-order) MMR. Since the distribution of period ratiosfor all planets is smooth (see Figures 7 and 8), the spikes in theobserved period ratio distribution come solely from this popula-tion. In our fully clustered model, we find that the fraction of allplanets near an MMR (as defined in §2.1.6) with another planet is f mmr = . ± .
04. Thus, reproducing the spikes near MMRs inthe observed period ratio distribution (without invoking a change inthe distribution of mutual inclinations) would require a populationof FSWP × f mmr × f σ i , high × (cid:104) NPPS (cid:105) = . + . − . planets per star As expected, we do not get any statistically significant spikes if we donot assign planets near an MMR to have low mutual inclinations. / in = ( P j + 2 / P j + 1 )/( P j + 1 / P j )0.000.050.100.150.20 F r a c t i o n Clustered P+RNon-clusteredKepler
Figure 9.
The distributions of ratios of period ratios, for adjacent observedplanets in three-planet chains. The grey shaded histogram shows the
Kepler distribution, while the solid blue and dotted red lines show those of theclustered periods and sizes and non-clustered models, respectively. None ofour models include any explicit correlation between period ratios, but ourclustered model appears to be more peaked than our non-clustered modeldue to the clustering in periods. Still, our clustered models are not as peakedas the data, suggesting that there is a significant correlation in the periodratios of observed planet pairs within the same system. that are near a first-order MMR. This can be compared to a rate ofFSWP × ( − f σ i , high ) × (cid:104) NPPS (cid:105) (cid:39) . + . − . planets per star that areresponsible for the “background” population of multiple planet sys-tems, excluding spikes in the period ratio distribution. We cautionthat these results are based on our distance function that includes aterm for the overall observed period ratio distribution, but does notknow about the dynamical significance of MMRs. Using a distancefunction that explicitly considers the behavior of the period-ratiodistribution near MMRs would be expected to provide a strongerconstraint. Conversely, at least part of the observed spikes in theperiod ratio distribution could be due to shifting the period ratiosof planet pairs from slightly less to slightly more than the value atresonance.The period ratio distribution also provides additional con-straints on proposals to explain the Kepler dichotomy. Ourparametrization of the mutual inclination distribution involves con-tributions to the period ratio distribution from both high and lowmutual inclination populations. While a four-planet system hasthree pairs of adjacent planets, a two-planet system has only oneplanet pair. Therefore, reducing the mean multiplicity of the high-inclination population from greater than four to two (see §4.4) wouldreduce the amplitude of spikes in the period ratio distribution bymore than a factor of three. If the high-inclination population wereto have a significantly lower mean number of planets per plane-tary system than the low-inclination population, then explaining thenear-resonant peaks in the distribution of period ratios would likelyrequire invoking an actual excess of planetary pairs near resonance.
Are planets evenly spaced?
Weiss et al. (2018a) and Zhu (2019)explored the period ratios of inner–outer pairs for observed three-planet chains. Weiss et al. (2018a) found a strong correlation be-tween P in = P j + / P j and P out = P j + / P j + (where planets aresorted so P j ’s are increasing) for small period ratios and concludedthat planets in 3+ systems tend to be uniformly spaced. However,Zhu (2019) suggested that this correlation is largely due to a combi- MNRAS , 1–31 (2019) He, Ford, and Ragozzine nation of detection biases and imposing an arbitrary upper limit onthe period ratio ( P < P in and P out as well as the their ratio( P out /P in ). Our model does not explicitly enforce any correlationbetween P in and P out either in the generation of systems or in thedistance function. Any observed correlation will be due to a com-bination of observational selection effects and the “natural” peak in P given the process for drawing orbital periods.In Figure 9, we plot the observed distribution of the ratio ofperiod ratios, P out /P in , for our fully clustered and non-clusteredmodels along with that of the Kepler data. The intrinsic distributionsof P out /P in (not shown) are peaked around one in each of ourmodels. However, the observed distribution of P out /P in is muchmore highly peaked than that observed in the non-clustered model(dotted red). In contrast, the clustered models produce a peak in P out /P in around unity, but it is slightly wider than the peak in theobserved distribution of P out /P in (solid blue). Furthermore, the Kepler distribution of the ratios of period ratios (shaded grey) iseven more peaked than the predictions of our clustered models.We conclude that there is a significant correlation in the periodratios of neighbouring planet pairs in actual planetary systems. Thisstrengthens the conclusion of Weiss et al. (2018a).
Future studies can build on our model and codes to incorporateadditional physical processes which are not included in our currentmodels into our forward model. One could consider alternative para-metric distributions for the distributions of intrinsic multiplicity,periods, radii, eccentricities, and mutual inclinations. For example,one possible extension would be to incorporate a more sophisticatedmodel for the distribution of planet radii (e.g., incorporating perioddependence and/or allowing for a local minimum in the frequencyas a function of planet size). Zink, Christiansen, & Hansen (2019)showed that adding a radius-valley resulted in only a modest effecton the inferred parameters for the radius and period distribution.Nevertheless, allowing for a more flexible model for the radius-period distribution would likely lead to improving the goodness-of-fit, particularly for the somewhat poorly fit distributions of transitdepth and transit depth ratios.Another possible extension would be to allow for explicitlymodelling a population of systems with resonant or near-resonantplanet pairs or chains. The results presented here show that the
Ke-pler observations (particularly the period ratio and transit durationratio distributions) can be explained by a model with no excessof planets near MMRs relative to the background period ratio dis-tribution. However, this finding does not preclude the possibilitythat there still might be a small excess of near-resonant planetarysystems also contributing to the observed spikes in the period ratiodistribution (see §4.5). Future studies could explore models with pa-rameters for the amplitude and/or width of such spikes in the periodratio distribution. Such additions would likely benefit from addingnew terms to the distance function that explicitly compare behaviornear MMRs. Improvements to the detection efficiency model thattake into account the presence of transit timing variations (TTVs)may also benefit from these efforts.This paper adopted a distance function that takes into accounteach of the observed distributions of the
Kepler population (as de-scribed in §2.4.1 and shown in Figure 2), along with the overallrate of planets per observed star f Kepler . Our summary statisticsare based on the marginal distributions for each observable. There-fore, it is possible that we might find models that reproduce our chosen summary statistics well, but differ in its prediction for theobserved
Kepler catalogue in other ways (e.g., a correlation be-tween orbital period and planet radius, or a correlation with planetoccurrence and host star metallicity). As shown in §3, our choiceof summary statistics and distance function already provide strongconstraints on many physically interesting model parameters. Fu-ture research adding additional summary statistics and componentdistances could shrink the ABC posterior further and help constrainany additional model parameters. For example, future research couldincorporate the detection of TTVs or transit duration variations(TDVs), so as to provide stronger constraints on the abundance andproperties of non-transiting planets.Our forward model makes use of the
Kepler
DR25 complete-ness and reliability products (Burke & Catanzarite 2017a,b,c; Chris-tiansen 2017; Thompson et al. 2018) to account for the vast majorityof factors that influence the detectability of planets. One aspect inwhich our model of the detection efficiency could be improved re-lates to how the
Kepler pipeline’s detection efficiency is affectedby multiple transiting planets. The pipeline first detects the planetwith the largest multiple event statistics (MES, analogous to com-bined transit signal-to-noise), masks out observations near whenthat planet transits, and re-searches the light curve for additionalplanets. Typically, only ∼ .
4% of the light curve is lost for eachadditional planet (Schmitt, Jenkins, & Fischer 2017; Zink, Chris-tiansen, & Hansen 2019). Thus, masking out observations neartransits of a more easily detected planet is expected to have only amodest effect on the integrated transit signal-to-noise and the proba-bility of observing at least three transits. Indeed, Zink, Christiansen,& Hansen (2019) showed that the multiplicity distribution is onlyslightly changed when accounting for the lower detection efficiencyof subsequent planets.
We have developed a framework for generating populations of plan-etary systems and simulating observed catalogues of exoplanetsunder the conditions of a
Kepler -like mission. We compare threephysically motivated models: the non-clustered model, the clusteredperiods model, and the clustered periods and sizes model, to the ob-served
Kepler
Q1-Q17 DR25 population of uniformly vetted planetcandidates. Our analysis is limited to planets with orbital periodsbetween 3 and 300 d and radii between 0.5 and 10 R ⊕ , which formthe bulk of the DR25 catalogue (2137 planets in our subset). Whilemost of our findings are consistent with previous works (Lissauer etal. 2011b; Fang & Margot 2012; Tremaine & Dong 2012; Fabryckyet al. 2014; Mulders et al. 2018), this work improves on previousstudies by incorporating improved and more homogeneous planetand star catalogues, an improved model for Kepler detection effi-ciency, and by providing a forward model that simultaneously fitsall of the marginal distributions described in §2.4.1.Our models are highly flexible and allow for the relatively fastgeneration ( ∼
10 s for a physical and observed catalogue with the 79935 target stars used in this study) of simulated observed cataloguesthat are similar to the
Kepler exoplanet catalogue in terms of manyobservables. We define a set of observable summary statistics and adistance function for comparing simulated models with the
Kepler data. We identify model parameters that result in simulated observedcatalogues that closely match the
Kepler catalogue in terms of theobserved planet multiplicity, period, period ratio, transit duration,period-normalized transit duration ratio ( ξ ), transit depth, and transitdepth ratio distributions. We train a Gaussian process (GP) model MNRAS000
Kepler catalogue in terms of theobserved planet multiplicity, period, period ratio, transit duration,period-normalized transit duration ratio ( ξ ), transit depth, and transitdepth ratio distributions. We train a Gaussian process (GP) model MNRAS000 , 1–31 (2019) lustered Planetary Systems for each of our physical models in order to build an emulator forrapidly predicting the distance as a function of the model parameters.Using the emulator, we draw from the ABC posterior and providecredible intervals for each of the model parameters as constrainedby the data.We find that a non-clustered model with a Poisson distributionfor the true number of planets per system, single and broken powerlaws for the periods and planetary radii, and a simple stability crite-ria, can reasonably fit most of the key observable properties of the Kepler population, including the overall rate of observed planets,as well as the period, transit depth, and transit duration distribu-tions, although there are clear differences. However, the number ofobserved multiplanet systems, the period ratio distribution, and theplanet radius ratio distributions are poorly modelled by the non-clustered model. In contrast, clustered models (of periods) can pro-duce observed planet multiplicity and period ratio distributions thatare significantly better fits to those observed by
Kepler . The transitdurations and duration ratios are also improved to an appreciableextent. Our fully clustered model, with clustered periods and sizes,performs similarly well and is the best description of the marginaldistributions of the
Kepler data to date. While this model allowsfor more similar sizes of planets within the same system, the fit tothe observed radius ratio distribution is not significantly improvedin terms of the KS or AD distances, suggesting that additional fea-tures in this distribution are still inadequately reproduced by any ofour current models. Most notably, the observed radius ratio distri-bution appears to be quite asymmetric, which is likely a signatureof stripped planetary atmospheres due to photoevaporation (Lopez,Fortney, & Miller 2012; Owen & Wu 2017) or core-powered mass-loss (Gupta & Schlichting 2018). Weiss et al. (2018a) find a verysimilar effect of clustering and asymmetry in the radius ratio dis-tribution and also a slight positive correlation between radius ratioand the difference in effective temperatures between adjacent plan-ets. Thus, we show with forward modelling that while the observedplanet radii ratio distribution of the
Kepler data require more thana simple clustering in intrinsic sizes to explain, the data cannot bereproduced in a non-clustered model by detection biases alone, incontrast to the results of Zhu (2019).Several investigations on the numbers of
Kepler single andmultitransiting systems, i.e. the observed multiplicity distribution,suggest that there is an apparent excess of observed singles and thatthis indicates the existence of more than one underlying populationof planetary systems, or a so-called “
Kepler dichotomy” (e.g., Lis-sauer et al. 2011b; Johansen et al. 2012; Hansen & Murray 2013;Ballard & Johnson 2016). While there is evidence for the excessof stars with a single transiting planet (relative to the predictionsbased on abundance of systems with multiple transiting planets),other studies (e.g., Weiss et al. 2018b) show that the stellar prop-erties of the singles and multis are indistinguishable, suggesting acommon origin. We account for the dichotomy by including twopopulations of planetary systems with separate mutual inclinationdispersions (Rayleigh ( σ i , high ) for a fraction f σ i , high of all systemsand Rayleigh ( σ i , low ) for the remaining systems) for each of our mod-els, and find that this is necessary to fit the multiplicity distribution.In our non-clustered model, the occurrence of multitransiting sys-tems is significantly underproduced, due to a tendency for f σ i , high to be very low ( (cid:46) − f σ i , high (cid:39) σ i , low ∼ . ◦ andthe remaining ∼
40% of systems with broad mutual inclinationsrequired to account for the excess observed singles. The log ξ distri-bution for planets not near any apparent MMRs with other observed planets is also well reproduced with these mutual inclinations and aRayleigh distribution of eccentricities, σ e ∼ .
02. Based on modelsthat allow for clustering in period and radius, it is unlikely that theexcess of single transiting planet systems could be explained solelyby a large population of intrinsically single-planet systems. Thus,our results provide new evidence in favour of the
Kepler dichotomybeing due to a population of planetary systems with a broader dis-tribution of mutual inclinations than characteristic of the observed
Kepler multiple planet systems.Previous studies also show that most observed
Kepler systemsare not near low-order MMRs (e.g., Lissauer et al. 2011b; Veras& Ford 2012; Fabrycky et al. 2014; Steffen & Hwang 2015), andour model reliably reproduces this observation. Nevertheless, thesmall fraction of systems near MMRs are particularly interesting todynamicists and for constraining planet formation models. There-fore, we investigated the ability of our model to reproduce the smallspikes in the period ratio distribution slightly wide of first-orderMMRs due purely to geometrical effects by assigning planet pairsnear MMRs with mutual inclinations also drawn from the populationwith a narrow distribution of mutual inclinations ( σ i , low ∼ . ◦ ).While some but not all of the draws from our clustered modelsinclude spikes in the period ratio distribution near MMRs, thesespikes are statistically robust (see the quantiles (red curves) in Fig-ures 5 and 6). Thus, we show that allowing the distribution of mutualinclinations to depend upon the proximity of a planet pair to MMRprovides an alternative explanation for the observed spikes in theperiod ratio distribution. This study did not consider and thus can-not exclude models in which the observed spikes in the period ratiodistribution are due to an intrinsically higher rate (e.g., due to mi-gration leading to resonant trapping). Previous studies have notedthat most Kepler systems are not in an MMR (e.g., Veras & Ford2012), and so an additional mechanism would be required to ex-plain why the systems are near, but not in MMR (e.g., Lithwick &Wu 2012; Lee et al. 2013; Petrovich, Malhotra, & Tremaine 2013;Delisle & Laskar 2014; Goldreich & Schlichting 2014; Xie 2014;Chatterjee & Ford 2015; Izidoro et al. 2017). Future studies may beable to distinguish between these two models and better constrainthe properties of near-MMR systems by making use of additionalobservational constraints (e.g., transit timing variations and transitduration variations). The results of precision radial velocity surveysmay also help constrain the rate of additional non-transiting planetcompanions.Our framework of simulating ensembles of planetary systemsin a forward model provides a way of directly probing the underly-ing populations of exoplanets and their properties, including planetsthat are not observed or do not transit. Our clustered models alsoprovide considerable utility for informing follow-up efforts of newexoplanet surveys such as the Transiting Exoplanet Survey Satellite(TESS) mission, which is expected to discover thousands of addi-tional planetary candidates (Sullivan et al. 2015; Bouma et al. 2017;Stassun et al. 2017). By matching the observed distributions of pe-riod ratios and transit depth ratios for the
Kepler exoplanets andassuming that the
Kepler population is representative of planetarysystems to be observed by TESS, our clustered periods and sizesmodel can be used to compute conditional probabilities of addi-tional RV detectable planets in systems with already known planetsof measured periods and radii.We have made the core SysSim code ( https://github.com/ExoJulia/ExoplanetsSysSim.jl ), inputs collated from numer-ous datafiles ( https://github.com/ExoJulia/SysSimData ),and the code specific to the clustered models ( https://github.com/ExoJulia/SysSimExClusters , Zenodo DOI 10.5281/zen-
MNRAS , 1–31 (2019) He, Ford, and Ragozzine odo.3386372) available to the public via Github. We encourageother researchers to contribute model extensions via pull requestsand/or additional public git repositories. Additionally, researcherscan use our forward modelling pipeline to perform detailed com-parisons of the results of planet formation simulations with
Kepler observations, so as to improve our understanding of planet forma-tion and the architectures of planetary systems more generally. Forresearchers who prefer not to run the SysSim code themselves, weprovide hundreds of physical and observed catalogues in table for-mat, each containing the simulated true and observed properties of ∼ planetary systems, generated from each of our three mod-els using either the best-fitting parameter values or draws from theABC-posterior distribution for parameter values explained in §3. ACKNOWLEDGEMENTS
We thank the entire
Kepler team for years of work leading to asuccessful mission and data products critical to this study. We ac-knowledge many valuable contributions with members of the
Ke-pler
Science Team’s working groups on multiple body systems,transit timing variations, and completeness working groups. Wethank Keir Ashby, Danley Hsu, Neal Munson, Shane Marcus, andRobert Morehead for contributions to the broader SysSim project.We thank Derek Bingham, Earl Lawrence, Ilya Mandell, OdedAahronson, Ben Bar-Oh, Dan Fabrycky, Jack Lissauer, Gijs Mul-ders, Aviv Ofir, and Jason Rowe for useful conversations. We thankRebekah Dawson for reading preliminary drafts of this paper andproviding detailed suggestions. MYH acknowledges the support ofthe Natural Sciences and Engineering Research Council of Canada(NSERC), funding reference number PGSD3-516712-2018. EBFand DR acknowledge support from NASA Origins of Solar Sys-tems grant
Kepler
Participating Scientist Program Cycle II grant corner.py package (Foreman-Mackey 2016). We acknowl-edge the Institute for CyberScience ( http://ics.psu.edu/ ) atThe Pennsylvania State University, including the CyberLAMP clus-ter supported by NSF grant MRI-1626251, for providing advancedcomputing resources and services that have contributed to the re-search results reported in this paper. This study benefited from the2013 SAMSI workshop on Modern Statistical and ComputationalMethods for Analysis of
Kepler
Data, the 2016/2017 Program onStatistical, Mathematical and Computational Methods for Astron-omy, and their associated working groups. This material was basedupon work partially supported by the National Science Foundationunder Grant DMS-1127914 to the Statistical and Applied Math-ematical Sciences Institute (SAMSI). Any opinions, findings, and conclusions or recommendations expressed in this material are thoseof the author(s) and do not necessarily reflect the views of the Na-tional Science Foundation.
REFERENCES
Anderson, T. W. & Darling, D. A. 1952, Ann. Math. Stat., 23, 193Andrae, R., et al. 2018, A&A, 616, A8Ballard, S. & Johnson, J. A. 2016, ApJ, 816, 66Batalha, N. M., et al. 2013, ApJS, 204, 24Bezanson, J., Edelman, A., Karpinski, S., & Shah, V. B. 2014, arXiv e-prints,arXiv:1411.1607Borucki, W. J., et al. 2010, Science, 327, 977Borucki, W. J., et al. 2011a, ApJ, 728, 117Borucki, W. J., et al. 2011b, ApJ, 736, 19Bouma, L. G., Winn, J. N., Kosiarek, J., McCullough, P. R., 2017,arXiv:1705.08891Brakensiek, J. & Ragozzine, D. 2016, ApJ, 821, 47Burke, C. J., & Catanzarite, J. 2017a, Planet Detection Metrics: Windowand One-Sigma Depth Functions for Data Release 25, Tech. rep. NASAAmes Research Center, Moffett Field, CA.Burke, C. J., & Catanzarite, J. 2017b, Planet Detection Metrics: Per-TargetFlux-Level Transit Injection Tests of TPS for Data Release 25, Tech.rep. NASA Ames Research Center, Moffett Field, CA.Burke, C. J., & Catanzarite, J. 2017c, Planet Detection Metrics: Per-TargetDetection Contours for Data Release 25, Tech. rep. NASA Ames Re-search Center, Moffett Field, CA.Carrera, D., Ford, E. B., Izidoro, A., Jontof-Hutter, D., Raymond, S. N.,Wolfgang, A., 2018, ApJ, 866, 104Chambers, J. E., Wetherill, G. W., & Boss, A. P. 1996, Icarus, 119, 261Chatterjee, S. & Ford, E. B. 2015, ApJ, 803, 33Chen, J. & Kipping, D. 2016, ApJ, 834, 17Christiansen, J. L. 2017, Planet Detection Metrics: Pixel-Level Transit In-jection Tests of Pipeline Detection Efficiency for Data Release 25, Tech.rep. NASA Ames Research Center, Moffett Field, CA.Ciardi, D. R., Fabrycky, D. C., Ford, E. B., Gautier, T. N., Howell, S. B.,Lissauer, J. J., Ragozzine, D., Rowe, J. F., 2013, ApJ, 763, 41Coughlin, J. L., et al. 2016, ApJS, 224, 12Coughlin, J. L. 2017, Planet Detection Metrics: Robovetter Completenessand Effectiveness for Data Release 25, Tech. rep. NASA Ames ResearchCenter, Moffett Field, CA.Cressie, N. & Read, T. R. C. 1984, J. R. Stat. Soc. B, 46, 440Dawson R. I., Lee E. J., Chiang E., 2016, ApJ, 822, 54Delisle, J.-B. & Laskar, J. 2014, A&A, 570, L7Fabrycky, D. C., et al. 2014, ApJ, 790, 146Fabrycky, D. C. & Winn, J. N. 2009, ApJ, 696, 1230Fang, J., & Margot, J.-L. 2012, ApJ, 761, 92Figueira, P., et al. 2012, A&A, 541, A139Ford, E. B., Quinn, S. N., & Veras, D. 2008, ApJ, 678, 1407Ford, E. B. 2014, Proc. Natl. Acad. Sci. USA, 111, 12616Ford, E. B., He, M. Y., Hsu, D. C., & Ragozzine, D. 2018b, PlanetarySystems Simulation & Model of Kepler Mission for Characterizingthe Occurrence Rates of Exoplanets and Planetary Architectures, v1.0,Zenodo, doi:10.5281/zenodo.1205172. https://doi.org/10.5281/zenodo.1205172
Foreman-Mackey, D. 2016, corner.py: Scatterplot matrices in Python, J.Open Source Softw., 1(2), 24, doi:10.21105/joss.00024Fulton, B. J., et al. 2017, AJ, 154, 109Gaia Collaboration, 2018, A&A, 616, A1Ginzburg, S., Schlichting, H. E., & Sari, R. 2016, ApJ, 825, 29Ginzburg, S., Schlichting, H. E., & Sari, R. 2016, MNRAS, 476, 759Gladman, B. 1993, Icarus, 106, 247Goldreich, P. & Schlichting, H. E. AJ, 147, 32Gupta, A. & Schlichting, H. E. 2019, MNRAS, 487, 24Hadden, S. & Lithwick, Y. 2014, ApJ, 787, 80Hansen, B. M. S. & Murray, N. 2013, ApJ, 775, 53Howard, A. W., et al. 2012, ApJS, 201, 15 MNRAS000
Foreman-Mackey, D. 2016, corner.py: Scatterplot matrices in Python, J.Open Source Softw., 1(2), 24, doi:10.21105/joss.00024Fulton, B. J., et al. 2017, AJ, 154, 109Gaia Collaboration, 2018, A&A, 616, A1Ginzburg, S., Schlichting, H. E., & Sari, R. 2016, ApJ, 825, 29Ginzburg, S., Schlichting, H. E., & Sari, R. 2016, MNRAS, 476, 759Gladman, B. 1993, Icarus, 106, 247Goldreich, P. & Schlichting, H. E. AJ, 147, 32Gupta, A. & Schlichting, H. E. 2019, MNRAS, 487, 24Hadden, S. & Lithwick, Y. 2014, ApJ, 787, 80Hansen, B. M. S. & Murray, N. 2013, ApJ, 775, 53Howard, A. W., et al. 2012, ApJS, 201, 15 MNRAS000 , 1–31 (2019) lustered Planetary Systems Hsu, D. C., Ford, E. B., Ragozzine, D., & Morehead, R. C., 2018, AJ, 155,205Hsu, D. C., Ford, E. B., Ragozzine, D., & Ashby, K. 2019, AJ, 158, 109Izidoro, A., Ogihara, M., Raymond, S. N., Morbidelli, A., Pierens, A., Bitsch,B., Cossou, C., Hersant, F., 2017, MNRAS, 470, 1750Johansen, A., Davies, M. B., Church, R. P., Holmelin, V., 2012, ApJ, 758,39Kipping, D. M. 2010, MNRAS, 407, 301Kolmogorov, A. N. 1933, Giornale dell’Istituto Italiano degli Attuari, 4, 83Latham, D. W., et al. 2011, ApJL, 732, L24Lee, M. H, Fabrycky, D., & Lin, D. N. C. 2013, ApJ, 774, 52Lissauer, J. J., et al. 2011a, Nature, 470, 53Lissauer, J. J., et al. 2011b, ApJS, 197, 8Lissauer, J. J., et al. 2012, ApJ, 750, 112Lissauer, J. J., et al. 2014, ApJ, 784, 44Lithwick, Y. & Wu, Y. 2012, ApJL, 756, L11Lopez, E. D. & Fortney, J. J. 2014, ApJ, 792, 1Lopez, E. D., Fortney, J. J., & Miller, N. 2012, ApJ, 761, 59Millholland, S., Wang, S., & Laughlin, G. 2017, ApJL, 849, L33Mills, S. M., et al. 2019, AJ, 157, 5Moorhead, A. V., et al. 2011, ApJS, 197, 1Morehead, R. C. 2016, PhD Dissertation, https://etda.libraries.psu.edu/catalog/mc87pq25t
Mulders, G. D., Pascucci, I., Apai, D., Ciesla, F. J., 2018, AJ, 156, 24Mullally, F., et al. 2015, ApJS, 217, 31Ning, B., Wolfgang, A., & Ghosh, S. 2018, ApJ, 869, 5O’Hagan, A. 2004, Rel. Eng. Sys. Safety, 91, 1290.Osada, R., Funkhouser, T., Chazelle, B., & Dobkin, D. 2002, ACM Trans.Graph., Vol. 21., No. 4, 807, 832,
Owen, J. E. & Wu, Y. 2013, Kepler Planets: A Tale of Evaporation, ApJ,775, 105Owen, J. E. & Wu, Y. 2017, ApJ, 847, 29Petigura, E. A., Marcy, G. W., & Howard, A. W. 2013b, ApJ, 770, 69Petrovich, C., Malhotra, R., & Tremaine, S. 2013, ApJ, 770, 24Pettitt, A. N. 1976, Biometrika, 63, 161Price, E. M. & Rogers, L. A. 2014, ApJ, 794, 92Pu, B. & Wu, Y. 2015, ApJ, 807, 44Ragozzine, D. & Holman, M. J. 2010, arXiv:1006.3727Rajpaul, V., Aigrain, S., Osborne, M. A., Reece, S., Roberts, S., 2015,MNRAS, 452, 2269Rasmussen, C. E. & Williams, C. K. I. 2006, Gaussian Processes for MachineLearning, MIT Press, Cambridge, MARowe, J. F., et al. 2014, ApJ, 784, 45Rowe, J. F., et al. 2015, ApJS, 217, 16Seager, S. & Mallén-Ornelas, G. 2003, ApJ, 585, 1038Shabram, M., Demory, B.-O., Cisewski, J., Ford, E. B., Rogers, L., 2015,ApJ, 820, 93Schmitt, J. R., Jenkins, J. M., & Fischer, D. A. 2017, ApJ, 153, 180Smirnov, N. 1948, Ann. Math. Stat., 19, 279Stassun, K. G., Oelkers, R. J., Pepper, J., Gaudi, B. S., 2017, AJ, 156, 102Steffen, J. H., et al. 2010, ApJ, 725, 1226Steffen, J. H. & Hwang, J. A. 2015, MNRAS, 448, 1956Sullivan, P. W., et al. 2015, ApJ, 809, 77Thompson, S. E., Mullally, F., Coughlin, J., Christiansen, J. L., Henze, C.E., Haas, M. R., Burke, C. J., 2015, ApJ, 812, 46Thompson, S. E., et al. 2018, ApJS, 235, 38Tremaine, S. & Dong, S. 2012, AJ, 143, 94Twicken, J. D., et al. 2016, AJ, 152, 6Van Eylen, V. & Albrecht, S. 2015, ApJ, 808, 126Van Eylen, V., Agentoft, C., Lundkvist, M. S., Kjeldsen, H., Owen, J. E.,Fulton, B. J., Petigura, E., Snellen, I., 2018, MNRAS, 479, 4786Van Eylen, V., et al. 2019, AJ, 157, 61Veras, D. & Ford, E. B. 2012, MNRAS, 420, L23Weiss, L. M., et al. 2018, AJ, 155, 48Weiss, L. M., et al. 2018, AJ, 156, 254Weiss, L. M. & Marcy, G. W. 2014, ApJ, 783, L6Weissbein, A., Steinberg, E., & Sari, R. 2012, arXiv:1203.6072 Winn, J. N. & Fabrycky, D. C. 2015, ARA&A, 53, 409Wolfgang, A. & Laughlin, G. 2012, ApJ, 750, 148Wolfgang, A., Rogers, L. A., & Ford, E. B. 2016, ApJ, 825, 19Wu, Y. & Lithwick, Y. 2013, ApJ, 772, 74Wu, D.-H., Zhang, R. C., Zhou, J.-L., Steffen, J. H., 2019, MNRAS, 484,1538Xie, J.-W. 2014, ApJ, 786, 153Xie, J.-W., et al. 2016, Proc. Natl. Acad. Sci. USA, 113, 11431Youdin, A. N. 2011, ApJ, 742, 38Zhu, W., Petrovich, C., Wu, Y., Dong, S., Xie, J., 2018, ApJ, 860, 101Zhu, W. 2019, arXiv:1907.02074Zink, J. K., Christiansen, J. L., & Hansen, B. M. S. 2019, MNRAS, 483,4479This paper has been typeset from a TEX/L A TEX file prepared by the author.MNRAS , 1–31 (2019) He, Ford, and Ragozzine f i ,high = 0.03 +0.020.02 . . . . l n ( p ) ln( p )=1.26 +0.160.15 . . . . P P = 0.13 +0.100.07 . . . . . R R =1.08 +0.310.20 . . . . . R R =5.30 +0.530.40 . . . . e e =0.003 +0.0030.002 i , h i g h () i , high ( )=62.23 +16.0949.59 .
03 0 .
06 0 .
09 0 . f i ,high . . . . i , l o w () .
00 1 .
25 1 .
50 1 . ln( p ) .
30 0 .
15 0 .
00 0 . P . . . . . R . . . . . R .
005 0 .
010 0 .
015 0 . e
20 40 60 80 i , high ( ) .
15 0 .
30 0 .
45 0 . i , low ( ) i , low ( )=0.31 +0.100.10 Figure 1.
Marginal posterior distributions of the model parameters for the non-clustered model, showing the projections of the points that pass a distancethreshold of D W , KS =
55 as drawn from the GP emulator. The prior mean function was set to a constant value of 75. MNRAS000
55 as drawn from the GP emulator. The prior mean function was set to a constant value of 75. MNRAS000 , 1–31 (2019) lustered Planetary Systems f i ,high = 0.42 +0.100.10 . . . . l n ( c ) ln( c )=0.39 +0.320.31 . . . . l n ( p ) ln( p )=1.50 +0.300.50 . . . . . P P = 0.13 +0.730.48 . . . . . R R =0.75 +0.460.46 . . . . . R R =4.91 +0.960.62 . . . . . e e =0.014 +0.0090.007 i , h i g h () i ,high ( )=49.40 +21.9520.73 . . . . . i , l o w () i ,low ( )=1.13 +0.380.33 .
30 0 .
45 0 . f i ,high . . . . . P . . . . ln( c ) . . . . ln( p ) . . . . . P . . . . . R . . . . . R .
008 0 .
016 0 .
024 0 .
032 0 . e
20 40 60 80 i ,high ( ) . . . . . i ,low ( ) .
12 0 .
16 0 .
20 0 .
24 0 . PP = 0.20 +0.040.03 Figure 2.
Marginal posterior distributions of the model parameters for the clustered periods model, showing the projections of the points that pass a distancethreshold of D W , KS =
35 as drawn from the GP emulator. The prior mean function was set to a constant value of 75.MNRAS , 1–31 (2019) He, Ford, and Ragozzine f i ,high = 0.41 +0.110.13 . . . . l n ( c ) ln( c )=0.28 +0.320.27 . . . . l n ( p ) ln( p )=1.27 +0.380.50 . . . . . P P = 0.01 +0.710.29 . . . . R R =0.85 +0.260.26 . . . . R R =5.49 +0.420.33 . . . . e e =0.010 +0.0080.006 i , h i g h () i ,high ( )=49.47 +24.5225.32 . . . . i , l o w () i ,low ( )=1.16 +0.310.26 .
15 0 .
30 0 .
45 0 . f i ,high . . . . . P . . . . ln( c ) . . . . ln( p ) . . . . . P . . . . R . . . . R .
008 0 .
016 0 .
024 0 . e
20 40 60 80 i ,high ( ) . . . . i ,low ( ) .
12 0 .
16 0 .
20 0 .
24 0 . PP = 0.17 +0.040.03 Figure 3.
Marginal posterior distributions of the model parameters for the clustered periods model, showing the projections of the points that pass a distancethreshold of D W , AD (cid:48) =
140 as drawn from the GP emulator. The prior mean function was set to a constant value of 250. MNRAS000
140 as drawn from the GP emulator. The prior mean function was set to a constant value of 250. MNRAS000 , 1–31 (2019) lustered Planetary Systems f i ,high = 0.40 +0.110.12 . . . . . l n ( c ) ln( c )=0.24 +0.330.40 . . . . . l n ( p ) ln( p )=0.73 +0.600.56 . . . . . P P = 0.07 +0.660.45 . . . . . R R =1.27 +0.260.25 . . . . R R =5.08 +0.710.54 . . . . e e =0.014 +0.0100.008 i , h i g h () i ,high ( )=49.26 +23.4725.31 . . . . . i , l o w () i ,low ( )=1.29 +0.350.32 . . . . . R R = 0.32 +0.070.07 .
15 0 .
30 0 .
45 0 . f i ,high . . . . . P . . . . . ln( c ) . . . . . ln( p ) . . . . . P . . . . . R . . . . R .
01 0 .
02 0 .
03 0 . e
20 40 60 80 i ,high ( ) . . . . . i ,low ( ) .
16 0 .
24 0 .
32 0 .
40 0 . R .
12 0 .
16 0 .
20 0 .
24 0 . PP = 0.20 +0.040.04 Figure 4.
Marginal posterior distributions of the model parameters for the clustered periods and sizes model, showing the projections of the points that pass adistance threshold of D W , AD (cid:48) =
140 as drawn from the GP emulator. The prior mean function was set to a constant value of 250.MNRAS , 1–31 (2019) He, Ford, and Ragozzine f i ,high = 0.32 +0.060.06 . . . . . l n ( c ) ln( c )=0.04 +0.460.45 . . . . l n ( p ) ln( p )=1.05 +0.390.34 . . . . . P P = 0.19 +0.690.57 . . . . . R R =1.28 +0.720.57 . . . . . R R =4.50 +1.030.70 . . . . . e e =0.019 +0.0110.009 i , h i g h () i ,high ( )=48.84 +18.8318.18 . . . . i , l o w () i ,low ( )=1.38 +0.430.38 . . . . . R R = 0.30 +0.070.07 . . . . f i ,high . . . . . P . . . . . ln( c ) . . . . ln( p ) . . . . . P . . . . . R . . . . . R .
01 0 .
02 0 .
03 0 .
04 0 . e
20 40 60 80 i ,high ( ) . . . . i ,low ( ) .
16 0 .
24 0 .
32 0 .
40 0 . R .
12 0 .
16 0 .
20 0 .
24 0 . PP = 0.21 +0.040.04 Figure 5.
Marginal posterior distributions of the model parameters for the alternative MMR inclinations model (clustered periods and sizes model without thelowering of mutual inclinations for planets near an MMR), showing the projections of the points that pass a distance threshold of D W , KS =
35 as drawn fromthe GP emulator. The prior mean function was set to a constant value of 75. MNRAS000
35 as drawn fromthe GP emulator. The prior mean function was set to a constant value of 75. MNRAS000 , 1–31 (2019) lustered Planetary Systems f i ,high = 0.31 +0.100.10 . . . . l n ( c ) ln( c )=0.35 +0.320.41 . . . . . l n ( p ) ln( p )=0.56 +0.500.43 . . . . . P P = 0.05 +0.370.33 . . . . R R =1.32 +0.300.29 . . . . R R =4.81 +0.790.71 . . . . e e =0.014 +0.0100.008 i , h i g h () i ,high ( )=42.32 +25.7120.23 . . . . . i , l o w () i ,low ( )=1.29 +0.340.31 . . . . . R R = 0.30 +0.060.06 .
15 0 .
30 0 .
45 0 . f i ,high . . . . . P . . . . ln( c ) . . . . . ln( p ) . . . . . P . . . . R . . . . R .
01 0 .
02 0 .
03 0 . e
20 40 60 80 i ,high ( ) . . . . . i ,low ( ) .
16 0 .
24 0 .
32 0 .
40 0 . R .
12 0 .
16 0 .
20 0 .
24 0 . PP = 0.18 +0.040.03 Figure 6.
Marginal posterior distributions of the model parameters for the alternative MMR inclinations model (clustered periods and sizes model without thelowering of mutual inclinations for planets near an MMR), showing the projections of the points that pass a distance threshold of D W , AD (cid:48) =
140 as drawnfrom the GP emulator. The prior mean function was set to a constant value of 250.MNRAS , 1–31 (2019) He, Ford, and Ragozzine
25 30 35 40 45 50 55 W (KS)0.00.2 F r a c t i o n Clustered P+RClustered PNon-clustered
200 300 400 500 600 W (AD )0.00.2 F r a c t i o n Clustered P+RClustered PNon-clustered D f F r a c t i o n wD f CRPD F r a c t i o n w CRPD KS for { P }0.00.1 F r a c t i o n w KS AD for { P }0.00.2 F r a c t i o n w AD KS for { }0.00.2 F r a c t i o n w KS AD for { }0.000.25 F r a c t i o n w AD KS for { }0.00.2 F r a c t i o n w KS AD for { }0.00.5 F r a c t i o n w AD KS for { i +1 / i }0.00.2 F r a c t i o n w KS AD for { i +1 / i }0.00.2 F r a c t i o n w AD KS for { t dur }0.00.2 F r a c t i o n w KS AD for { t dur }0.00.2 F r a c t i o n w AD KS for { res }0.00.1 F r a c t i o n w KS AD for { res }0.00.2 F r a c t i o n w AD KS for { non res }0.000.25 F r a c t i o n w KS AD for { non res }0.00.5 F r a c t i o n w AD Figure 7.
Histograms of the total (weighted; top row) and individual distance terms (second row and below) for each of our models. Each histogram includes 100simulated catalogues with parameters drawn from the emulator passing our KS distance thresholds for each model. With the exception of the two upper-mostrows, the left-hand panels show the individual KS distance terms while the right-hand panels show the individual AD distance terms. For each panel, the bottomx-axis labels the true KS or AD distances ( D KS or D AD (cid:48) ), while the top x-axis labels the weighted distance (i.e. the bottom axis scaled by a weight w as listedin Table 2). MNRAS000