The Exoplanet Population Observation Simulator. I  The Inner Edges of Planetary Systems
Bestfit parameters
Bestfit parameters
Parameters
Abstract
The Kepler survey provides a statistical census of planetary systems out to the habitable zone. Because most planets are nontransiting, orbital architectures are best estimated using simulated observations of ensemble populations. Here, we introduce EPOS, the Exoplanet Population Observation Simulator, to estimate the prevalence and orbital architectures of multiplanet systems based on the latest Kepler data release, DR25. We estimate that at least of sunlike stars have nearly coplanar planetary systems with 7 or more exoplanets. The fraction of stars with at least one planet within au could be as high as depending on assumptions about the distribution of single transiting planets. We estimate an occurrence rate of planets in the habitable zone around sunlike stars of . The innermost planets in multiplanet systems are clustered around an orbital period of 10 days (0.1 au), reminiscent of the protoplanetary disk inner edge or could be explained by a planet trap at that location. Only a small fraction of planetary systems have the innermost planet at long orbital periods, with fewer than and having no planet interior to the orbit of Mercury and Venus, respectively. These results reinforce the view that the solar system is not a typical planetary system, but an outlier among the distribution of known exoplanetary systems. We predict that at least half of the habitable zone exoplanets are accompanied by (nontransiting) planets at shorter orbital periods, hence knowledge of a closein exoplanet could be used as a way to optimize the search for Earthsize planets in the Habitable Zone with future direct imaging missions.
Gijs D. Mulders
1 Introduction
The Kepler exoplanet survey has revolutionized our understanding of planetary systems around other stars (e.g. Borucki, 2017). During the first four years of the mission, hereafter Kepler, it discovered thousands of exoplanets and exoplanet candidates, the majority of which are smaller than Neptune and orbit close to their hosts stars. Because transit surveys can detect only a small fraction of the exoplanets, a completeness correction is necessary to understand the true population of exoplanets (e.g. Batalha, 2014).
The planet occurrence rate, defined as the average number of planets per star, is typically found to be of order unity for planets in the orbital period and planet radius range detectabe with Kepler (Youdin, 2011; Fressin et al., 2013). Much progress has been made in recent years in understanding the survey detection efficiency of the Kepler mission (e.g. Christiansen et al., 2015, 2016). While earthsized planets in the habitable zone of sunlike stars are just below the detection limits of Kepler (Thompson et al., 2018), occurrence rates of earthsized planets in the habitable zone, , can be estimated by considering smaller M and K dwarf stars (Dressing & Charbonneau, 2013, 2015) or by extrapolating from shorter orbital periods and/or larger planet radii (Petigura et al., 2013a; ForemanMackey et al., 2014; Burke et al., 2015).
The exoplanet population discovered by the Kepler mission has been characterized in great detail. Key findings include that planet occurrence rates increase with decreasing planet size (Howard et al., 2012) and remains roughly constant for planets smaller than (Petigura et al., 2013b; Morton & Swift, 2014; Mulders et al., 2015b). Recently, a gap in the planet radius distribution has been identified around (Fulton et al., 2017; Van Eylen et al., 2017). The occurrence of subNeptunes increases with distance from the host star up to an orbital period of days, after which it remains roughly constant (Howard et al., 2012; Mulders et al., 2015a), in contrast to planets larger than Neptune whose occurrence increases with orbital period (Dong & Zhu, 2013).
The presence of multiple transiting planets around the same star provides important additional constraints on planetary orbital architectures. The large number of observed systems with multiple transiting planets indicate that planets are preferentially located in systems with small mutual inclinations (Lissauer et al., 2011, 2012; Fang & Margot, 2012). Their orbital spacings follow a peaked distribution without a clear preference for being in orbital resonances (Fabrycky et al., 2014). However, the majority of planets in multiplanet system are not transiting, and modifying planetary occurrence rates from those based on the observed planetary architectures to account for nontransiting or otherwise undetected planets is not straightforward. The true multiplicity of planetary systems can only be constrained from the observed multiplicity and making assumptions about their mutual inclination distribution (Tremaine & Dong 2012, see also Brakensiek & Ragozzine 2016). Ensemble populations of systems with six or more planets with mutual inclinations of a few degrees provide a good match to the observed population of multiplanet systems (Lissauer et al., 2011; Fang & Margot, 2012; Fabrycky et al., 2014; Ballard & Johnson, 2016). However, these simulations underpredict the number of single transiting systems by a factor 2, indicating that an additional populations of planets must be present that is either intrinsically single or has high mutual inclinations.
The occurrence rates and planetary architectures place strong constraints on planet formation models (e.g. Hansen & Murray, 2013; Dawson et al., 2016). However, the quantitative comparison of planet formation and orbital evolution models and the observed exoplanet population has been greatly complicated by observational biases. Therefore, published comparisons resorted to trying to explain the presence of specific planet subpopulations by identifying processes that can give rise to such planets. While these studies are essential, such qualitative comparisons can not make full use of the wealth of information represented by the overall exoplanet population statistics. To address this limitation we introduce EPOS ^{1}^{1}1https://github.com/GijsMulders/epos, the Exoplanet Population Observation Simulator, that provides a Python interface comparing planet population synthesis models to exoplanet survey data, taking into account the detection biases that skew the observed distributions of planet properties. EPOS uses a forward modeling approach to simulate observable exoplanet populations and constrain their properties via Markov Chain Monte Carlo simulation using emcee (ForemanMackey et al., 2013), and has already been employed in two different studies (Kopparapu et al., 2018; Pascucci et al., 2018).
In this paper, we verify this approach using parametric models of planet populations, that we compare to the final data release of the Kepler mission, DR25. From this, we are able to make a statistical evaluation of the properties of exoplanetary systems. Among other results we report on the location of the innermost planet in planetary systems to place our Solar System in the context of exoplanet systems. We examine how the innermost planets in most of the Kepler systems are located much closer in than Mercury and Venus. In an upcoming paper, we will use EPOS to make a direct comparison between planet formation models and exoplanet populations.
2 Code description
The Exoplanet Population Observation Simulator, EPOS, employs a forward modeling approach to constrain exoplanet populations through the following iterative procedure:
 Step 1

Define a distribution of planetary systems from analytic forms for planet occurrence rates and planetary architectures or from a planet population synthesis model.
 Step 2

From this derive a transiting planet population by assigning random orientations to each system and evaluating which planets transit their host stars.
 Step 3

Determine which of the transiting planets would be detected by Kepler, accounting for detection efficiency.
 Step 4

Compare the detectable planet population with exoplanet survey data using a summary statistic.
 Step 5

Repeat steps 14 until the simulated detectable planet population matches the observed planet population to constrain the intrinsic distribution of planetary systems.
In this section we will describe each step in greater detail. Steps 14 take less than a second to run on a singlecore CPU, allowing for an efficient sampling of the parameter space using MCMC methods. Figure 1 summarizes these steps and their quantitative implementation in a flowchart. A description of all mathematical symbols used in this paper can be found in the Appendix, Table 3.
2.1 Step 1: Planet Distributions
Planet populations are generated using a Monte Carlo simulation by random draws from a multidimensional probability distribution function, , which represents the intrinsic distribution of planets and planetary systems. We simulate a planet survey equal in size to the Kepler survey, roughly stars. The parametric descriptions of are based on studies of planet occurrence rate and planet multiplicity from Kepler. Distributions based on planet formation models will be described in an upcoming paper.
Here we describe two parametric descriptions for the planet population that correspond to two different modes in EPOS:
 Occurrence rate mode

Simulates only the distribution of planet radius, , and orbital period, , as . This mode is similar to occurrence rate calculations that estimate the average number of planets per star, , in a certain period and radius range. This mode can also be used to estimate the planet mass distribution for comparison with microlensing data, see Pascucci et al. (2018) for details.
 Multiplanet mode

Simulates multiple planets per system, taking into account their relative spacing and mutual inclination, , to estimate orbital architectures. The properties of each planet in the system are drawn from a distribution , where is the planet’s orbital inclination and is an index for each planet in the system.
2.1.1 Occurrence rate mode
The central assumption we use in this paper is that the planet occurrence rate distribution is a separable function in planet radius and orbital period,
(1) 
The distribution is normalized such that integral over the simulated period and planet radius range equals the number of planets per star, :
(2) 
The planet orbital period distribution is described by a broken power law
(3) 
where the break reflects the observed flattening of subNeptune occurrence rates around an orbital period of ten days (e.g. Howard et al., 2012; Mulders et al., 2015a). In this example we use for the power law index of the increase interior to days, and a flat occurrence rate , at longer periods but we will refine these values later.
The planet radius distribution function is similarly described by a broken power law
(4) 
reflecting the observed increase of small planets with decreasing radius (), and the “plateau” () of constant occurrence rates for planets smaller than (e.g. Howard et al., 2012; Petigura et al., 2013b).
The probability distribution function of planet occurrence thus has 7 free parameters (, , , , , , ). An example for is shown in Figure 2.
2.1.2 Multiplanet mode
In its multiplanet mode EPOS simulates a planetary system for each star. Each planet in the system, denoted by subscript , is characterized by its radius, orbital period, and orbital inclination with respect to the observer, , and there are planets per system. Because the properties of planets in multiplanet systems tend to be correlated, we draw a typical set of planet properties for each system () on which the properties of each planet in the system () are dependent:
(5)  
where we make the assumption that distributions of planet size, period, and inclination are not interdependent. Previous studies have shown that planets in multiplanet systems tend to be more similar in size and more regularly spaced than random pairings from the overall distribution (Millholland et al., 2017; Weiss et al., 2018) and have low mutual inclinations of a few degrees (e.g. Fabrycky et al., 2014). This motivates our choice to not draw the properties of each planet independently from a distribution .
The system properties () are drawn from a distribution
(6) 
where once again we made the assumption that the distributions of planet size, period, and inclination are separable functions of each variable. The distribution of is normalized such that the integral over the simulated parameter range is equal to the fraction of stars with planetary systems, :
(7) 
The parameterization of the system and planet properties are described below. The system architecture is illustrated in Figure 3
Orbital Inclination
A planet’s orbital inclination with the respect to the line of sight, , is dependent on the inclination of the system, , the mutual inclination of the planet with respect to that system, , and the longitude of the ascending node, .
We assume that the planet mutual inclinations, , follow a Rayleigh distribution:
(8) 
where is the mode of the mutual inclination distribution of planetary orbits, which is typically (Fabrycky et al., 2014). The distribution of system inclinations with respect to the observer, , is proportionate to .
(9) 
where is an edgeon orbit. The distribution of planet inclinations with respect to the observer is then
(10) 
where is the longitude of the ascending node with respect to the observer.
Orbital Period
We assume that the planets are regularly spaced in the logarithm of the orbital period, with the location of the first planet () defining the location of additional planets () in the system. The period ratio of adjacent planets is denoted by where is always the planet with the larger orbital period.
The location of the first planet in the system, (), is parameterized by the broken power–law distribution described by Equation 3, but with a steeper decline in planet occurrence with orbital period () that we will constrain in the fitting process.
Observed orbital spacings follow a broad range in period ratios, and we follow Malhotra (2015) in parameterizing the distribution of dimensionless spacings, , as a lognormal distribution. The orbital period distribution of the kth planet in the system is then given by
(11) 
with respect to the orbital period of the previous planet, . and are free parameters that characterize the median and width of the distribution, with typical values of and .
Planet Radius
We assume all planets in the system are of equal size
(12) 
where we assume the radius distribution follows Equation 4. We choose equalsized planets over randomly assigned sizes because planets in multiplanet systems tend to be of similar size (Lissauer et al., 2011; Millholland et al., 2017; Weiss et al., 2018). While there are observed trends of increasing planet size with orbital period, these trends are strongest at short orbital periods (e.g. Carrera et al., 2018) and for large planet sizes (Ciardi et al., 2013; Huang et al., 2016) that we will exclude from our observational comparison, see §2.4.
Thus the number of free parameters describing the population of planetary systems is 11 (, , , , , , , , , , ).
2.2 Step 2: Transiting Planet Populations
The synthetic planet population is simulated using a Monte Carlo approach by sampling the distribution function of planet properties (Eq. 1 or 5) times. We simulate planetary systems in the period range days and (Figure 4). For each planet, 2 random numbers are drawn to determine its properties. The first random number is used to determine the planet orbital period from the cumulative distribution function of Eq. 3. The second random number is used to determine the planet radius from the cumulative distribution function of Eq. 4.
The total number of draws in the synthetic survey is where is the number of stars in the survey (), multiplied by the average number of planets per star, . In occurrence rate mode, with , the simulated sample consists of planets. In multiplanet mode, with and , the simulated sample consists of the same number of planets distributed across systems. Each star in the survey is assigned a unique identifier, ID, to keep track of the observable planet multiplicity.
In occurrence rate mode, the subset of transiting planets can be simulated by considering the geometric transit probability
(13) 
The semimajor axis is calculated using Kepler’s third law. The stellar mass and radius are the average values of the surveyed stars (see next section). A planet is transiting if
(14) 
where is a continuous random variable between and . The transiting planet sample, , typically contains planets. An example is shown in Figure 4.
2.2.1 Multiplanet transit probability
The Monte Carlo simulation of transiting planetary systems goes as follows. First, the system inclination, is drawn according to Equation 9. The draw from a distribution is performed by generating a random number between 0 and 1 and interpolate the parameter value from the normalized cumulative distribution. Next, the mutual inclination of the orbit of each planet in the system, , is drawn from Equation 8. Then the longitude of the ascending node with respect to the observer, , is randomly drawn for each planet from a uniform distribution, . The inclination, of each planets orbit with respect to the lineofsight can then be calculated from equation 10. The transiting planet population, , is defined by planets that traverse the stellar limb, given by
(15) 
where is the geometric transit probability from Eq. 13. We note that for a single planet (), this expression is equivalent to the geometric transit probability of Eq. 13. This expression is only valid for small mutual inclinations.
To account for the Kepler dichotomy, the apparent excess of single transiting systems (e.g. Ballard & Johnson, 2016), we assume that a fraction of planetary systems has an isotropic distribution of orbits described by Eq. 13 instead of Eq. 15. Typically, , but we will treat it as a free parameter to be constrained from the data.
2.3 Step 3: Survey Detection Efficiency
We simulate a detectable planet sample, , from the simulated sample of transiting planets, , by taking into account the survey detection efficiency, , and the vetting completeness, , described below and displayed in Figure 5. A planet is detectable if:
(16) 
where is a continuous random variable between 0 and 1. is the detection efficiency based on the combined signaltonoise ratio of all planet transits. is the detection efficiency of the Kepler Robovetter for a planet candidate sample with high reliability.
We calculate detection efficiency contours for each individual star observed by Kepler using KeplerPORTs^{2}^{2}2https://github.com/nasa/KeplerPORTs (Burke & Catanzarite, 2017). We included all stars that were fully searched by the Kepler pipeline, for which stellar properties were available in the Mathur et al. (2017) catalog, and for which all detection metrics were available on the exoplanet archive^{3}^{3}3https://exoplanetarchive.ipac.caltech.edu/docs/Kepler_completeness_reliability.html.
The detection efficiencies were evaluated on a grid with 20 logarithmicallyspaced orbital period bins between and days assuming planets on circular orbits, and 21 logarithmically spaced planet radius bins between and . After removing giant and subgiant stars according to the effective temperature dependent surface gravity criterion in Huber et al. (2016), the sample consists of stars, with a median mass of and median radius of . We calculate the survey detection efficiency, , by averaging the individual contributions of each star, displayed in the upper panel of Figure 5.
Not all transiting planets that are detectable based on the detection efficiency are vetted as planet candidates by the Kepler Robovetter (Coughlin et al., 2016). Including only reliable planet candidates, here defined with a disposition score larger than , further reduces the vetting completeness. The vetting efficiencies were calculated following Thompson et al. (2018) based on the Kepler simulated data products^{4}^{4}4https://exoplanetarchive.ipac.caltech.edu/docs/KeplerSimulated.html. The vetting efficiency was evaluated on the same radius and period grid as the detection efficiency and using the same criterion to select mainsequence stars. We use a powerlaw in planet radius and a broken powerlaw in orbital period to obtain a smooth function for the vetting completeness
(17) 
with , , , , (see Appendix B for details). The bestfit vetting efficiency is shown in Figure 5. The vetting efficiency varies between close to 100% for Neptunesized planets at short orbital periods to 25% near the habitable zone. The break in the powerlaw is needed to match the reduced vetting efficiency for planets at orbital periods larger than days as described in Thompson et al. (2018).
The detectable planet sample, , typically contains about planets (Figure 4). In multiplanet mode, EPOS also keeps track of which planets are observable as part of multiplanet systems (Figure 6). From the observed multiplanet systems we generate a set of summary statistics: , the number of stars with detectable planets; , the orbital period ratio between adjacent planets (Fig. 7); , the orbital period of the innermost planet in the system (Fig. 8).
It is worth noting that the properties of the observable multiplanet systems are significantly different from the distributions from which they are generated. When the intrinsic population contains only planetary systems with planets, the majority of systems have transiting detectable planets. In addition, the observable period and period ratio distributions are skewed by detection biases.
The geometric transit probability favors the detection of planets at shorter orbital periods. The distribution of the location of the innermost observed planet in the systems peaks at days compared to days in the intrinsic distribution (Fig. 8). We also find that planetary systems at very short orbital periods ( day) are overrepresented in the detectable distribution by an order of magnitude, while systems with orbital periods similar to the terrestrial planets ( days) are underrepresented by an order of magnitude, due to the same detection biases (Fig. 8).
Pairs of planets with smaller orbital period ratios are more likely to be both transiting, shifting the peak of the observable orbital period ratio distribution to smaller period ratios (Fig. 7). The hatched area in the bottom panel of Figure 7 shows simulated planet pairs that are observable as adjacent but have a nontransiting planet between them. These planet pairs dominate the period ratio distribution at large orbital period ratios ().
2.4 Step 4: Observational Comparison
In this step we compare the simulated planet populations to the observed exoplanet properties to evaluate how well the simulated planet population reproduces the collection of observed planetary systems. We generate a set of summary statistics for both the simulated data and the Kepler survey. We evaluate the survey as a whole, and do not consider dependencies on stellar properties such as stellar mass (Mulders et al., 2015a, b) or metallicity (Mulders et al., 2016; Petigura et al., 2018).
In occurrence rate mode, the summary statistics are the planet radius distribution, , the orbital period distribution , and the total number of planets, . In multiplanet mode, additional summary statistics are calculated for the number of stars with planets, , the period ratio distribution between adjacent planets, , and the period of the innermost planet in the system, . We evaluate this summary statistic in the range of and days. These ranges exclude two regions where the assumptions of separability of parameters clearly breaks down. The maximum planet size of is chosen to exclude giant planets, which have a different distribution of orbital periods than subNeptunes (Dong & Zhu, 2013; Santerne et al., 2016), are less often part of multiplanet systems (Steffen et al., 2012), or have very dissimilar sizes than other planets in the system (Huang et al., 2016). The minimum orbital period of 2 days is chosen to exclude the photoevaporation desert (e.g. Lundkvist et al., 2016) where the planetradius distribution deviates significantly from that at larger orbital periods. The other bounds are chosen because there are very few planet detection outside this range.
We then generate a summary statistic from the Kepler DR25 catalog. The planet candidate list is taken from Thompson et al. (2018). We include only main sequence planet hosts by removing giant and subgiant stars according to the effective temperature dependent surface gravity criterion in (Huber et al., 2016). We also use a disposition score cut of 0.9 to select a more reliable sample of planet candidates (see Thompson et al. 2018 for details). We account for the lower completeness of this highreliability planet sample by explicitely taking into account the vetting completeness in the calculation of the survey detection efficiency.
The final list containing 3,041 planet candidates is shown in Figure 9. The list contains 1,840 observed single systems and 324 double, 113 triple, 38 quadruple, 10 quintuples, and 2 sextuple systems within the region where the summary statistic is evaluated.
2.4.1 Occurrence rate mode
Figure 10 shows how the summary statistics of planet radius and orbital period are generated from the detectable planet population (blue) and from the Kepler exoplanet population (orange). We compare the planet radius distributions and orbital period distributions separately. While this approach ignores any covariances between planet radius and orbital period that are present in the Kepler data, it is consistent with the assumption made in Eq. 1 that these functions are separable. We quantify the distance between the two distributions using the twosample KolmogorovSmirnoff (KS) test and calculated the associated probabilities, and , that the observed and simulated distributions are drawn from the same data. We will minimize the differences between these distributions in the fitting step.
2.4.2 Multiplanet statistics
We calculate three additional summary statistics for multiplanet systems: the frequency of multiplanet systems (, Fig. 11); the period ratio of adjacent planet pairs (, Fig. 12); and the distribution of the locations of the innermost planet in each system (, Fig. 13), all evaluated within and . As discussed in the previous section, the observed planet populations are subject to detection biases, and a proper comparison requires comparing the observed distribution to the simulated distribution, which is what EPOS does. All summary statistics are calculated in the same way for observed planets (orange lines) and simulated planets (blue histograms) from the set of orbital periods and host star IDs. We use 2sample KS tests to calculate the probabilities that the distribution of inner planet orbital periods () and orbital period ratios () are drawn from the same distribution as the observations. We use a Pearson test to calculate the probability, , that the multiplanet frequencies of the simulated sample are drawn from the same distribution as the observations.
The frequency distribution of planets in multiplanet systems is shown in Figure 11 for a model with planets per system, a mode for the mutual inclination of , and . The hatched region indicates simulated planets that are on isotropic orbits to match the excess of singletransiting systems. We also show a model with only coplanar orbits ( and , green) that overpredicts the frequency of multiplanet systems, particularly at high numbers of planets per system. A model with planets on isotropic orbits (, red) underpredict the frequency of all multiplanet systems.
The period ratio of adjacent planet pairs is shown in Figure 12. The blue histogram shows the observable period ratios of the example model where planets are regularly spaced according to the period ratio distribution of Eq. 11. The shape of the period ratio distribution qualitatively reproduces the observed distribution (orange), with a peak near a period ratio of and a tail towards large orbital period ratios. For comparison, the red line shows a simulation where planetary systems are not regularly spaced, which is constructed by randomly drawing planets from the periodradius distribution of Eq. 1 (Figure 2). These simulated systems have an observable period ratio distribution that is much wider than observed, indicating that planetary systems are regularly spaced (See also Weiss et al. 2018).
The distribution of the location of the innermost detected planet in each multiplanet systems is shown in Figure 13. We only focus on the innermost planet of a detected multiplanet system, as for observed singleplanet systems it can not be derived from the observable if they are intrinsically single or part of a multiplanet system with nontransiting planets. The distribution of innermost detected planets constrains the fraction of planetary systems without planets at short ( days) orbital periods.
2.5 Step 5: Fitting Procedure
With the framework to compare the parameterized distributions of exoplanets to observables we proceed to constrain the parameters to provide the best match to the observed planetary systems. The runtime of steps 14 is less than a tenth of a second in occurrence rate mode and less than a second in multiplanet mode, allowing for an comprehensive sampling of parameter space. We use emcee (ForemanMackey et al., 2013), an opensource Python implementation of the algorithm by Goodman & Weare (2010) to sample the parameter space and estimate the posterior distribution of parameters.
2.5.1 Occurrence rate mode
In occurrence rate mode, we explore the 7dimensional parameter space consisting of the number of planets per star, , the orbital period distribution, , , , and the planet radius distribution, , , . The summary statistics to evaluate the model given the observations are the number of detected planets, ; the planet orbital period distribution, ; and the planet radius distribution . We use Fisher’s method (Fisher, 1925) to combine the probabilities from the summary statistic into a single parameter:
(18) 
which we use as the likelihood of the model given the data in emcee, e.g. .
Parameter  Example  Bestfit 

(days)  
()  
We run emcee with 200 walkers for 5,000 iterations, allowing for a 1000step burnin to reach convergence. For the initial positions of the walkers we use the parameters of the example model, see Table 1. The MCMC chain and parameter covariances are shown in the Appendix (Figures 22 and 23). The bestfit values and their confidence intervals are calculated as the , and and percentiles, respectively. The simulated model for the bestfit parameters and for 30 samples of the posterior is shown in Figure 14.
The posterior distribution of planet radius and orbital period are shown in Figure 15. As a sanity check, we calculate occurrence rates as function of planet radius and orbital period using the inverse detection efficiency method. The occurrence per bin is calculate as
(19) 
where is the survey completeness evaluated at the radius and orbital period of each planet in the bin (Figure 16). The posterior distributions provide a decent match to the binned occurrence rates, with three notable deviations. At days and occurrence rates are lower than the posterior. We attribute this to these bins including regions where Kepler has not detected planet candidates (see Figure 16) and occurrence rates are therefore underestimated. The broken power–law in planet radius does not describe the population of giant planets at , and we therefore in the following restrict our observational comparison to .
Parameter  Example  Bestfit 

(days)  
()  
()  
2.6 Multiplanet systems
In multiplanet mode, we explore the 9dimensional parameter space consisting of: the fraction of stars with planetary systems, ; the mode of the mutual inclination distribution, ; the orbital period distribution of the inner planet, , , ; the period ratio distribution, and ; and the fraction of isotropic systems, . The summary statistics to evaluate the model given the observations are the frequency of multiplanet systems, ; the orbital period distribution, ; the planet orbital period ratio distribution, ; the distribution of the innermost planet in the system,. We use Fisher’s method to combine the probabilities from the summary statistic into a single parameter:
(20) 
which we use as the likelihood of the model given the data in emcee, e.g., .
The planet radius distribution is fixed to minimize the amount of free parameters, and we use the bestfit values constrained in the previous section and listed in Table 2. The number of planets per system, , is fixed to because it is not well constrained in the fitting: Systems with fewer than 6 planets do not reproduce the observed multiplicity distribution as there are two detected sextuple systems (). Systems with more than 7 planets have the same observation signature as the 8th, 9th, etc., planet in the system will typically remain undetected, as the combined likelihood of detecting all planets in such systems is too small to allow detections given the size of the Kepler sample.
We sample the parameter space with emcee using 100 walkers for 2,000 iterations, allowing for a 500step burnin. The MCMC chain and parameter covariances are shown in the appendix (Figure 24). The simulated observables of the bestfit models and of 30 samples from the posterior are shown in Figure 17. We will discuss the results in the next section.
3 Results
By using EPOS in parametric mode find that the planet occurrence rate of the Kepler mission are well described by a broken power–law in orbital period and planet radius in the region days and . We estimate the planet occurrence rate, the average number of planets per star, to be in this regime and in the simulated range ( days and ). This number is higher than the occurrence rate calculated from inverse detection efficiencies (), which underestimate the occurrence rates for small planets at long orbital periods. These results are largely consistent with previous occurrence rate studies collected in the SAG13^{5}^{5}5https://exoplanets.nasa.gov/system/internal_resources/details/original/680_SAG13_closeout_8.3.17.pdf literature study on planet occurrence rates (see Kopparapu et al. 2018 for details).
The observed population of exoplanets is well described by a population of regularlyspaced planetary systems with 7 or more planets that orbits of stars. We find that the mutual inclinations are consistent with the Kepler dichotomy: Half the planet population is in planetary systems that have nearly coplanar orbits, here described by a Rayleigh distribution with , consistent with previous estimates that find ranges between (Lissauer et al., 2011; Fang & Margot, 2012; Fabrycky et al., 2014; Ballard & Johnson, 2016). The other half of planets appear as singleplanet transiting systems, which we model as multiplanet systems with isotropically distributed orbital inclinations. The typical orbital period ratio between adjacent planets is but with a wide distribution, consistent with previous analysis of spacings between planets (Fang & Margot, 2013).
We discover that the inner edge of planetary systems are clustered around an orbital period of 10 days (Fig. 18), reminiscent of the protoplanetary disk inner edge which is located at au for premainsequence sunlike stars (MillanGabet et al., 2007). While the break in planet occurrence rate around 10 days has been previously connected to the disk inner edge (Mulders et al., 2015a; Lee & Chiang, 2017), this is the first time a peak in the occurrence rate distribution of subNeptunes has been identified. The posterior distribution of the innermost planet peaks at an orbital period of days and decays towards shorter and longer orbital period with a power–law index of and , respectively. The decay towards long orbital periods is surprising, since planet occurrence rates are constant or slightly increasing in this range (green line). Planets exterior to the break are therefore mostly the 2nd, 3rd, etc., planet in the system. We also show that systems with inner planets at orbital periods of days are intrinsically rare, and we discuss the implications for planet formation theories and the origins of the solar system below.
4 Discussion
4.1 How rare is the solar system?
Our modeling analysis indicates that multiplanet systems without planets interior to days are intrinsically rare. We estimate, based on parametric distributions of planet parameters, that of planetary systems have no planet interior to the orbit of Mercury and have no planet interior to Venus. This implies that solar system may simply be in the tail of the distribution of exoplanet systems.
On the other hand, this comparison with the solar system is made under the assumption that planet orbital architectures (inclinations, spacings, inner planet location) are independent of radius and orbital period. This assumption will need to be verified with additional data. The orbital architectures of the Kepler planetary systems are most constrained by planets that are larger and closer in, typically days and (see Figure 14). A solar system analogue, if detected, would most likely appear as a single transiting systems in the Kepler data due the low probability that multiple terrestrial planets transit (e.g. Brakensiek & Ragozzine, 2016). Hence, we can not rule out that the solar system may be part of a population of planetary systems with different orbital characteristics than those detected by the Kepler mission.
We make a direct estimate of how many planetary systems without planets interior to Mercury or Venus could be present in the Kepler data, without extrapolating the distribution of orbital periods from closerin planets. We generate a population of planetary systems with similar orbital properties as the solar system where we place the innermost planet at the orbital period of mercury (or Venus) but otherwise keep the same parameters as in the bestfit model. Such a simulated planet population with and matches planet occurrence rates exterior to days for days. The simulated observations predict multiplanet systems with an innermost planet exterior to days, while only are observed. This indicates that if a population of planetary systems without planets interior to Mercury existed it would be detectable in the Kepler data. We do note that most detectable planets in this range are larger than one earth radius, so the data do not directly constrain true solar system analogues. Using a mixture of coplanar systems and highly inclined systems (which appear as intrinsically single), we can rule out that more than of systems () have days and .
We repeat this exercise for systems where the innermost planet shares the orbital period of Venus ( days). None of these simulated systems would be detectable as multiplanet systems, and no planetary systems with an orbital period larger than days are detected with Kepler. Hence we can not rule out that the Kepler exoplanet population contains a significant population of multiplanet systems without planets interior to Venus, though we do not find evidence that such a population exists.
Overall, the solar system seems to be a typical planetary system in most diagnostics of orbital architectures (planetary radii, periods, period ratios, inclinations, and planets per system, Figure 19). The only clear difference is in the period of the innermost planet where the solar system is an outlier. While the lack of superearths/miniNeptunes in the solar system is also notable (Martin & Livio, 2015), earthsized planets are not intrinsically rare, and the size of the terrestrial planets does not make the solar system an outlier in the exoplanet distribution.
4.2 Systems with habitable zone planets
We also estimate , the number of earthsized planets with earthlike orbital periods (“habitable zone” planets), here defined as and based on the Kopparapu et al. (2013) conservative habitable zone for a sunlike star that is representative of the most common star in the Kepler sample. By integrating the posterior distribution over the radius and orbital period range we find or . These results are consistent with the estimate of for GK dwarfs by Burke et al. (2015) though we note that there is a large dispersion in the literature^{6}^{6}6https://exoplanets.nasa.gov/system/internal_resources/details/original/680_SAG13_closeout_8.3.17.pdf on habitable zone planet occurrence rates based on the adopted completeness correction, the method of extrapolation into the habitable zone, and the planet sample selection. The confidence intervals on include counting statistics and systematic uncertainties in extrapolating the planet radius and orbital period distribution out to the habitable zone. They do not include systematic uncertainties on the adopted stellar parameters, in particular the stellar radii which directly impacts the planet radii. For example, large uncertainties in the stellar radius of unresolved stellar binaries lead to overestimating the occurrence of small planets (Ciardi et al., 2015; Silburt et al., 2015). Better observational constraints on the stellar properties are expected based on parallax measurments from the ESO Gaia mission. As a consistency check, we have repeated all calculations in this paper with the improved stellar radii and giant star classification from (Berger et al., 2018). The bestfit parameteric distributions are consistent with those in table 1 within errors, and the habitable zone planet occurrence rate of is not significantly different.
The habitable zone planet occurrence rates we estimate is consistent with that of Burke et al. (2015) based on the Q1Q16 catalog, but about a factor four higher than of the earliest estimate based on Kepler data from Petigura et al. (2013a), which we attribute to an improved understanding of detection efficiency of Kepler. This rate is comparable to that of M dwarfs estimated with inverse detection efficiency methods (Dressing & Charbonneau, 2015), though we leave a comparison between M dwarfs and FGK stars with consistent methodology for a future paper.
The majority of simulated planets at long orbital periods ( days) are not intrinsically single, but are part of multiplanet systems where the innermost planets are typically at an orbital period of days. While the multiplanet statistics do not directly constrain planetary systems of earthsized planets in the habitable zone of sunlike stars, the presence of planets on orbital periods of 10 days does provide a way to identify targets for future direct imaging missions, such as the HabEx and LUVOIR mission concepts.
Kipping & Lam (2017) have shown that the presence of transiting planets at short orbital periods increases the chance of finding transiting planets at larger orbital periods. Our analysis indicates that, because most planets in multiplanet systems are not transiting, the probability of finding any planet (transiting or nontransiting) is even higher. For example, the detection of an earthmass planets at an orbital period of ten days with radial velocity measurements indicates a higher probability of finding a planet in the habitable zone. Depending on whether singletransiting planets are intrinsically single or highlyinclined multiples, of systems with closein planets may also have planets in the habitable zone.
4.3 Protoplanetary disk inner edges
The peak in the location of the innermost planet in the system points to a preferred location of planet formation/migration in protoplanetary disks. This location may reflect either the inner edge of the region where planets form (Chiang & Laughlin, 2013; Mulders et al., 2015a; Lee & Chiang, 2017) or that of a planet trap where planet migration stalls (Terquem & Papaloizou, 2007). While the break in planet occurrence at days has been attributed to the disk inner edge before (Mulders et al., 2015a; Lee & Chiang, 2017), the peak in the location of the innermost planet presents solid evidence that there is a preferred location for planet formation at 0.1 au, rather than a continuum of location between 0.11 au. In particular, this result is inconsistent with a wide regions between 0.1 and 1 au acting as a trap for individual planets (e.g. Dittkrist et al., 2014), and is more reminiscent of inward migration of multiple planets where the first planet is trapped at the inner disk edge and halts the migration of other planets (Cossou et al., 2014; Coleman & Nelson, 2016; Izidoro et al., 2017).
4.4 Do most stars have planets?
Our analysis confirms that the average number of planets per (sunlike) star is larger than one (§2.1.1). However, this does not imply that most stars have planets, since planets can be unevenly distributed among stars. In (§2.1.1) we show that of stars can host multiplanet systems. An important assumption we made is that the apparent excess of single transits is because of multiplanet systems have high mutual inclinations. The inclination distribution of exoplanets can not be uniquely constrained from transit survey data (Tremaine & Dong, 2012). Different assumptions on the mutual inclination distribution lead to different estimates of the fraction of stars with planets, illustrated in Figure 20.
On the conservative end, the highlyinclined planets in otherwise coplanar multiplanet systems could contribute to the observed number of single transits (Figure 20, right panel). In this case, the dichotomy is an artifact of the assumption that mutual inclinations follow a Rayleigh distribution. A distribution with a larger tail towards higher inclinations could possibly fit the multiplanet statistics with a single population. In this case, the fraction of stars with planets drops to . Alternatively, a population of intrinsically single planets could contribute to the observed single transiting systems (Figure 20, right panel). In this case, the estimated fraction of stars with planets increases to unity.
A way of discriminating between both scenarios may be the use of stellar obliquity. Transiting planets with high mutual inclinations or isotropic distributions have large obliquities with respect to the line of sight. There are indications that the stellar obliquity is larger for singletransit systems (Morton & Winn, 2014), indicating that they may have highly inclined orbits. Several mechanisms have been proposed to disrupt initially coplanar systems, reducing their multiplicity and/or increasing their mutual inclinations, and giving rise to the population of single transiting systems. These mechanisms include dynamical instabilities within systems (Volk & Gladman, 2015; Pu & Wu, 2015), external perturbations from giant planets or stars (Johansen et al., 2012; Lai & Pu, 2017; Hansen, 2017; Izidoro et al., 2017; Cai et al., 2018) and hoststar misalignment (Spalding & Batygin, 2016). Alternatively, planets may form as intrinsically single systems in half the cases (Moriarty & Ballard, 2016). Additional constraints on the single planet population are needed to constrain these scenarios. However, we stress that roughly half of the observed single transiting systems are part of multiplanet systems and have nontransiting or nondetected planets in the system, and hence the sample of “singles” is diluted with multiplanets, weakening any potential trends. Transit timing variations can also be used to detect additional nontransiting planets in the system (Nesvorný et al., 2012).
The Kepler mission only detects planets out to au and down to . Hence, any estimate of the fraction of stars with planets is likely to be a lower limit. We do not see a clear indication that planet occurrence rate decreases towards the detection limits of Kepler. While planet occurrence rates calculated using the inverse detection efficiency seem to decrease for earthsize and smaller planets (Fig 15), this is also the region where the Kepler exoplanet list becomes incomplete. The posterior planet size distribution does not show this trend and is almost flat in log , indicating that planets below the detection efficiency of Kepler could be extremely common. Since planetary systems tend to have planets of similar size (Millholland et al., 2017; Weiss et al., 2018), planetary systems with only planets smaller than the detection limit would be missed. Extending the planet size distribution down to the size of Ceres () would increase the fraction of stars with planets to .
We find that the orbital period distribution of subNeptunes is nearly flat in logarithm of orbital period. However, planets at long orbital periods are mostly members of systems with planets also at shorter orbital periods: there is no evidence for a large population of planetary systems with only longperiod planets. Extrapolating the distribution of planetary systems to orbital periods larger than a year would add planets to existing systems, increasing the number of planets per star but not the fraction of stars with planets. Of course, such an extrapolation does not take into account the presence of additional populations of planets, which could be important, for example, if planets form more frequently at the snow line (exterior to au). Giant planets orbit of sunlike stars (Cumming et al., 2008), most of them beyond 1 au. Microlensing surveys hint at the existence of a population of Neptunemass planets around a significant fraction of M dwarfs (Cassan et al., 2012; Suzuki et al., 2016). Since these populations remain mostly undetected in the Kepler survey, it is not clear if they belong to the same stars, though there are some indications that they do (Huang et al., 2016).
5 Summary
We present the Exoplanet Population Observation Simulator, EPOS, a Python code to constrain the properties of exoplanet populations from biased survey data. We showcase how planet occurrence rates and orbital architectures can be constrained from the latest Kepler data release, DR25. We find that:

The Kepler exoplanet population between orbital periods of days and planet radii of are well described by broken power–laws. The planet occurrence rate in this regime is consistent with previous works. The estimated planet occurrence rate in the habitable zone is .

The observed multiplanet frequencies are consistent with ensemble populations that are a mix of nearly coplanar systems and a population of planets whose orbital architectures are unconstrained, consistent with the previously reported Kepler dichotomy. of exoplanets are in systems with 6 or more planets, orbiting of sunlike stars. The remaining of exoplanets could be intrinsically single planets or be part of multiplanet systems with high mutual inclinations, raising the fraction of stars with planets to somewhere between and .

The mutual inclinations of planetary orbits can be described by a Rayleigh distribution with a mode of . The spacing between adjacent planets follows a wide distribution with a peak at an orbital period ratio , though detection biases shift the peak of the observed distribution to shorter orbital period ratios.

The distribution of the innermost planet in the system peaks at an orbital period of days. Planetary systems without planets interior to Mercury’s and Venus’ orbit are rare at and , respectively.
The EPOS code presented in this paper provides a first step in a larger effort to constrain planetary orbital architectures from exoplanet populations. The next step will be to use more complex models of planetary system architectures based on planet formation models. Future directions include incorporating additional exoplanet survey data from radial velocity, microlensing, and direct imaging, as well as from transit surveys such as K2 and TESS.
References
 Ballard & Johnson (2016) Ballard, S., & Johnson, J. A. 2016, The Astrophysical Journal, 816, 66
 Batalha (2014) Batalha, N. 2014, in Proceedings of the National Academy of Sciences, 12647–12654
 Berger et al. (2018) Berger, T. A., Huber, D., Gaidos, E., & van Saders, J. L. 2018, arXiv.org, arXiv:1805.00231
 Borucki (2017) Borucki, W. J. 2017, in Proceedings of the American Philosophical Society, Astrobiology and Space Research Directorate Ames Research Center, NASA, 38–65
 Brakensiek & Ragozzine (2016) Brakensiek, J., & Ragozzine, D. 2016, The Astrophysical Journal, 821, 47
 Burke & Catanzarite (2017) Burke, C. J., & Catanzarite, J. 2017, Kepler Science Document
 Burke et al. (2015) Burke, C. J., Christiansen, J. L., Mullally, F., et al. 2015, The Astrophysical Journal, 809, 8
 Cai et al. (2018) Cai, M. X., Portegies Zwart, S., & van Elteren, A. 2018, Monthly Notices of the Royal Astronomical Society, 474, 5114
 Carrera et al. (2018) Carrera, D., Ford, E. B., Izidoro, A., et al. 2018, arXiv.org, arXiv:1804.05069
 Cassan et al. (2012) Cassan, A., Kubas, D., Beaulieu, J. P., et al. 2012, Nature, 481, 167
 Chiang & Laughlin (2013) Chiang, E., & Laughlin, G. 2013, Monthly Notices of the Royal Astronomical Society, 431, 3444
 Christiansen (2017) Christiansen, J. L. 2017, Kepler Science Document
 Christiansen et al. (2015) Christiansen, J. L., Clarke, B. D., Burke, C. J., et al. 2015, The Astrophysical Journal, 810, 95
 Christiansen et al. (2016) —. 2016, The Astrophysical Journal, 828, 99
 Ciardi et al. (2015) Ciardi, D. R., Beichman, C. A., Horch, E. P., & Howell, S. B. 2015, The Astrophysical Journal, 805, 16
 Ciardi et al. (2013) Ciardi, D. R., Fabrycky, D. C., Ford, E. B., et al. 2013, The Astrophysical Journal, 763, 41
 Coleman & Nelson (2016) Coleman, G. A. L., & Nelson, R. P. 2016, Monthly Notices of the Royal Astronomical Society, 457, 2480
 Cossou et al. (2014) Cossou, C., Raymond, S. N., Hersant, F., & Pierens, A. 2014, Astronomy & Astrophysics, 569, A56
 Coughlin (2017) Coughlin, J. L. 2017, KSCI19114001
 Coughlin et al. (2016) Coughlin, J. L., Rowe, J. F., Caldwell, D. A., et al. 2016, The Astrophysical Journal Supplement Series, 224, 12
 Cumming et al. (2008) Cumming, A., Butler, R. P., Marcy, G. W., et al. 2008, Publications of the Astronomical Society of Pacific, 120, 531
 Dawson et al. (2016) Dawson, R. I., Lee, E. J., & Chiang, E. 2016, The Astrophysical Journal, 822, 54
 Dittkrist et al. (2014) Dittkrist, K. M., Mordasini, C., Klahr, H., Alibert, Y., & Henning, T. 2014, Astronomy & Astrophysics, 567, A121
 Dong & Zhu (2013) Dong, S., & Zhu, Z. 2013, The Astrophysical Journal, 778, 53
 Dressing & Charbonneau (2013) Dressing, C., & Charbonneau, D. 2013, The Astrophysical Journal, 767, 95
 Dressing & Charbonneau (2015) —. 2015, The Astrophysical Journal, 807, 45
 Fabrycky et al. (2014) Fabrycky, D. C., Lissauer, J. J., Ragozzine, D., et al. 2014, The Astrophysical Journal, 790, 146
 Fang & Margot (2012) Fang, J., & Margot, J.L. 2012, The Astrophysical Journal, 761, 92
 Fang & Margot (2013) —. 2013, The Astrophysical Journal, 767, 115
 Fisher (1925) Fisher, R. 1925, Statistical methods for research workers (Edinburgh Oliver & Boyd)
 ForemanMackey (2016) ForemanMackey, D. 2016, The Journal of Open Source Software, 24, doi:10.21105/joss.00024
 ForemanMackey et al. (2013) ForemanMackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, Publications of the Astronomical Society of Pacific, 125, 306
 ForemanMackey et al. (2014) ForemanMackey, D., Hogg, D. W., & Morton, T. D. 2014, The Astrophysical Journal, 795, 64
 Fressin et al. (2013) Fressin, F., Torres, G., Charbonneau, D., et al. 2013, The Astrophysical Journal, 766, 81
 Fulton et al. (2017) Fulton, B. J., Petigura, E. A., Howard, A. W., et al. 2017, The Astronomical Journal, 154, 109
 Goodman & Weare (2010) Goodman, J., & Weare, J. 2010, Communications in Applied Mathematics and Computational Science, 5, 65
 Hansen (2017) Hansen, B. M. S. 2017, Monthly Notices of the Royal Astronomical Society, 467, 1531
 Hansen & Murray (2013) Hansen, B. M. S., & Murray, N. 2013, The Astrophysical Journal, 775, 53
 Howard et al. (2012) Howard, A. W., Marcy, G. W., Bryson, S. T., et al. 2012, The Astrophysical Journal Supplement, 201, 15
 Huang et al. (2016) Huang, C., Wu, Y., & Triaud, A. H. M. J. 2016, The Astrophysical Journal, 825, 98
 Huber et al. (2016) Huber, D., Bryson, S. T., Haas, M. R., et al. 2016, The Astrophysical Journal Supplement Series, 224, 2
 Hunter (2007) Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90
 Izidoro et al. (2017) Izidoro, A., Ogihara, M., Raymond, S. N., et al. 2017, Monthly Notices of the Royal Astronomical Society, 470, 1750
 Johansen et al. (2012) Johansen, A., Davies, M. B., Church, R. P., & Holmelin, V. 2012, The Astrophysical Journal, 758, 39
 Jones et al. (2001–) Jones, E., Oliphant, T., Peterson, P., et al. 2001–, SciPy: Open source scientific tools for Python
 Kipping & Lam (2017) Kipping, D. M., & Lam, C. 2017, Monthly Notices of the Royal Astronomical Society, 465, 3495
 Kopparapu et al. (2013) Kopparapu, R. K., Ramirez, R., Kasting, J. F., et al. 2013, Astrophysical Journal, 765, 131
 Kopparapu et al. (2018) Kopparapu, R. K., Hebrard, E., Belikov, R., et al. 2018, arXiv.org, arXiv:1802.09602
 Lai & Pu (2017) Lai, D., & Pu, B. 2017, The Astronomical Journal, 153, 42
 Lee & Chiang (2017) Lee, E. J., & Chiang, E. 2017, The Astrophysical Journal, 842, 40
 Lissauer et al. (2011) Lissauer, J. J., Ragozzine, D., Fabrycky, D. C., et al. 2011, The Astrophysical Journal Supplement, 197, 8
 Lissauer et al. (2012) Lissauer, J. J., Marcy, G. W., Rowe, J. F., et al. 2012, The Astrophysical Journal, 750, 112
 Lundkvist et al. (2016) Lundkvist, M. S., Kjeldsen, H., Albrecht, S., et al. 2016, Nature Communications, 7, 11201
 Malhotra (2015) Malhotra, R. 2015, The Astrophysical Journal, 808, 71
 Martin & Livio (2015) Martin, R. G., & Livio, M. 2015, The Astrophysical Journal, 810, 105
 Mathur et al. (2017) Mathur, S., Huber, D., Batalha, N., et al. 2017, The Astrophysical Journal Supplement Series, 229, 30
 MillanGabet et al. (2007) MillanGabet, R., Malbet, F., Akeson, R., et al. 2007, Protostars and Planets V, 539
 Millholland et al. (2017) Millholland, S., Wang, S., & Laughlin, G. 2017, The Astrophysical Journal Letters, 849, L33
 Moriarty & Ballard (2016) Moriarty, J., & Ballard, S. 2016, The Astrophysical Journal, 832, 34
 Morton & Swift (2014) Morton, T. D., & Swift, J. 2014, The Astrophysical Journal, 791, 10
 Morton & Winn (2014) Morton, T. D., & Winn, J. N. 2014, The Astrophysical Journal, 796, 47
 Mulders (2018) Mulders, G. D. 2018, EPOS: The Exoplanet Population Observation Simulator, doi:10.5281/zenodo.1247569
 Mulders et al. (2015a) Mulders, G. D., Pascucci, I., & Apai, D. 2015a, The Astrophysical Journal, 798, 112
 Mulders et al. (2015b) —. 2015b, The Astrophysical Journal, 814, 130
 Mulders et al. (2016) Mulders, G. D., Pascucci, I., Apai, D., Frasca, A., & MolendaZakowicz, J. 2016, The Astronomical Journal, 152, 187
 Nesvorný et al. (2012) Nesvorný, D., Kipping, D. M., Buchhave, L. A., et al. 2012, Science, 336, 1133
 Pascucci et al. (2018) Pascucci, I., Mulders, G. D., Gould, A., & Fernandes, R. 2018, arXiv.org, arXiv:1803.00777
 Petigura et al. (2013a) Petigura, E. A., Howard, A. W., & Marcy, G. W. 2013a, in Proceedings of the National Academy of Sciences, 19273–19278
 Petigura et al. (2013b) Petigura, E. A., Marcy, G. W., & Howard, A. W. 2013b, The Astrophysical Journal, 770, 69
 Petigura et al. (2018) Petigura, E. A., Marcy, G. W., Winn, J. N., et al. 2018, The Astronomical Journal, 155, 89
 Pu & Wu (2015) Pu, B., & Wu, Y. 2015, The Astrophysical Journal, 807, 44
 Santerne et al. (2016) Santerne, A., Moutou, C., Tsantaki, M., et al. 2016, Astronomy & Astrophysics, 587, A64
 Silburt et al. (2015) Silburt, A., Gaidos, E., & Wu, Y. 2015, The Astrophysical Journal, 799, 180
 Spalding & Batygin (2016) Spalding, C., & Batygin, K. 2016, The Astrophysical Journal, 830, 5
 Steffen et al. (2012) Steffen, J. H., Ragozzine, D., Fabrycky, D. C., et al. 2012, in Proceedings of the National Academy of Sciences, 7982–7987
 Suzuki et al. (2016) Suzuki, D., Bennett, D. P., Sumi, T., et al. 2016, The Astrophysical Journal, 833, 145
 Terquem & Papaloizou (2007) Terquem, C., & Papaloizou, J. 2007, The Astrophysical Journal, 654, 1110
 The Astropy Collaboration et al. (2018) The Astropy Collaboration, PriceWhelan, A. M., Sipőcz, B. M., et al. 2018, ArXiv eprints, arXiv:1801.02634
 Thompson et al. (2018) Thompson, S. E., Coughlin, J. L., Hoffman, K., et al. 2018, The Astrophysical Journal Supplement Series, 235, 38
 Tremaine & Dong (2012) Tremaine, S., & Dong, S. 2012, The Astronomical Journal, 143, 94
 van der Walt et al. (2011) van der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, Computing in Science & Engineering, 13, 22
 Van Eylen et al. (2017) Van Eylen, V., Agentoft, C., Lundkvist, M. S., et al. 2017, arXiv.org, arXiv:1710.05398
 Volk & Gladman (2015) Volk, K., & Gladman, B. 2015, The Astrophysical Journal Letters, 806, L26
 Weiss et al. (2018) Weiss, L. M., Marcy, G. W., Petigura, E. A., et al. 2018, The Astronomical Journal, 155, 48
 Youdin (2011) Youdin, A. N. 2011, The Astrophysical Journal, 742, 38
Appendix A Table of Mathematical Symbols
All mathematical symbols used in this paper are summarized in table 3.
Parameter  Unit  Description 
Planet radius  
day  Planet orbital period  
Orbital period ratio of adjacent planets  
Orbital inclination  
Planet mutual inclination  
Mode of mutual inclination distribution  
rad  Longitude of ascending node  
Planet probability density function  
Planetary system probability density function  
Continuous random variable  
Survey completeness  
Geometric transit probability  
Detection efficiency  
Vetting efficiency  
Fraction of isotropic systems  
Average number of planets per star  
Fraction of stars with planets  
Planets per system (multiplicity)  
ID  Host star identifier  
{}  Ensemble of parameters  
n  Sample size  
N  Number of detected planets  
Subscript for planetary systems  
Subscript for planets  
Subscript for transiting planets  
Subscript for detectable planets  
in  Subscript for innermost planet  
Index for kth planet in system  
Power law indices  
Normalization factor  
Dimensionless spacing parameters 
Appendix B Vetting Completeness
The Kepler Robovetter (Coughlin et al., 2016; Coughlin, 2017) was designed to discriminate between real transiting events and instrumental or astrophysical signals that mimic that of a transiting planet, a process called “vetting”. For planets with few transits and a low signaltonoise ratio, vetting is a tradeoff between detecting planets (completeness) or removing false positives (reliability). A planet candidate sample with high completeness has a low reliability, and vice versa. The Robovetter was calibrated by evaluating the ability to identify the signals of a number of simulated transiting events to better quantify the efficiency and confidence with which planets can be detected.
Each planet candidate is assigned a disposition score between and , which can be used to create a more reliable sample (fewer false positives classified as planet candidates) at the the costs of having a less complete sample (more planets classified as false positives), see Thompson et al. (2018) for details. For EPOS we opt to use a planet sample with high reliability by using a disposition score cut of , reducing the sample completeness which, however, can be quantified and corrected for as we discuss below. The score cut of was chosen to eliminate the peak of false positives in the habitable zone that coincides with the orbital period of the spacecraft.
An empirical procedure for estimating the vetting completeness is described in Thompson et al. (2018), based on injection of simulated transits into the Kepler data (Christiansen, 2017) and analyzed by the Robovetter (Coughlin, 2017). We downloaded the Robovetter results for simulated planet injections (group 1) from the Kepler simulated data page^{7}^{7}7https://exoplanetarchive.ipac.caltech.edu/docs/KeplerSimulated.html. We removed all injections not on mainsequence stars according to the prescription in Huber et al. (2016). We only considered injections inside the orbital period and radius range where we use EPOS ( days and ). There are injections, of which () are classified as planet candidates and () have a disposition score .
The top left panel of Figure 21 shows the vetting completeness, the fraction of injections classified as planet candidates with a score larger than , binned to the same grid as the survey detection efficiency. Planet injections are not equally distributed across this parameter range, with the majority of injections around subNeptunes at long orbital periods, as shown in the top right panel. The vetting completeness is at orbital periods less than roughly days and decreases below at the longest orbital periods. The completeness also drops below for giant planets, though this is most likely an artifact of those planets only being injected around the most noisy stars, see Thompson et al. (2018), and does not reflect a real decrease in detectability of giant planets in the Kepler pipeline.
We parametrize the vetting completeness to obtain a smooth distribution of radius and orbital period that covers all of the regions we are simulating with EPOS. We use a functional form for the vetting completenss, , that is a double broken powerlaw:
(B1) 
with
(B2) 
and
(B3) 
We find the bestfit parameters using scipy.optimize.curve_fit, where we minimize the difference between the measured and parametrized vetting completeness over all grid cells. The bestfit solution is shown in the bottom left panel of Figure 21 and the bottom right panel shows the residuals per grid cell. The bestfit parameters are , days, , , , , and . The break in the powerlaw with orbital period corresponds to the decreased detection efficiency at long orbital periods as documented in Thompson et al. (2018); Coughlin (2017). The break in the powerlaw of planet radius corresponds to a decreased detection efficiency for giant planets. However, as mentioned above, this is likely an artifact of the way the injection tests were defined, and there is no reason to assume that the low vetting efficiency for giant planets is a real feature (Thompson et al., 2018). Therefore, we do no include the powerlaw index radius break in EPOS. Instead, we assume the planet radius dependence of the detection is described by the powerlaw in planet radius that was fitted to the small planets. This yields a vetting efficiency of
(B4) 
with , which is shown in Figure 5.