Exoplanet Population Observation Simulator

The Exoplanet Population Observation Simulator. I - The Inner Edges of Planetary Systems

Gijs D. Mulders Lunar and Planetary Laboratory, The University of Arizona, Tucson, AZ 85721, USA Earths in Other Solar Systems Team, NASA Nexus for Exoplanet System Science Ilaria Pascucci Lunar and Planetary Laboratory, The University of Arizona, Tucson, AZ 85721, USA Max Planck Institute for Astronomy, K’́onigstuhl 17, Heidelberg, D-69117, Germany Earths in Other Solar Systems Team, NASA Nexus for Exoplanet System Science Dániel Apai Department of Astronomy, The University of Arizona, Tucson, AZ 85721, USA Lunar and Planetary Laboratory, The University of Arizona, Tucson, AZ 85721, USA Max Planck Institute for Astronomy, K’́onigstuhl 17, Heidelberg, D-69117, Germany Earths in Other Solar Systems Team, NASA Nexus for Exoplanet System Science Fred J. Ciesla Department of the Geophysical Sciences, The University of Chicago, 5734 South Ellis Avenue, Chicago, IL 60637 Earths in Other Solar Systems Team, NASA Nexus for Exoplanet System Science

Best-fit parameters

Gijs D. Mulders Lunar and Planetary Laboratory, The University of Arizona, Tucson, AZ 85721, USA Earths in Other Solar Systems Team, NASA Nexus for Exoplanet System Science Ilaria Pascucci Lunar and Planetary Laboratory, The University of Arizona, Tucson, AZ 85721, USA Max Planck Institute for Astronomy, K’́onigstuhl 17, Heidelberg, D-69117, Germany Earths in Other Solar Systems Team, NASA Nexus for Exoplanet System Science Dániel Apai Department of Astronomy, The University of Arizona, Tucson, AZ 85721, USA Lunar and Planetary Laboratory, The University of Arizona, Tucson, AZ 85721, USA Max Planck Institute for Astronomy, K’́onigstuhl 17, Heidelberg, D-69117, Germany Earths in Other Solar Systems Team, NASA Nexus for Exoplanet System Science Fred J. Ciesla Department of the Geophysical Sciences, The University of Chicago, 5734 South Ellis Avenue, Chicago, IL 60637 Earths in Other Solar Systems Team, NASA Nexus for Exoplanet System Science

Best-fit parameters

Gijs D. Mulders Lunar and Planetary Laboratory, The University of Arizona, Tucson, AZ 85721, USA Earths in Other Solar Systems Team, NASA Nexus for Exoplanet System Science Ilaria Pascucci Lunar and Planetary Laboratory, The University of Arizona, Tucson, AZ 85721, USA Max Planck Institute for Astronomy, K’́onigstuhl 17, Heidelberg, D-69117, Germany Earths in Other Solar Systems Team, NASA Nexus for Exoplanet System Science Dániel Apai Department of Astronomy, The University of Arizona, Tucson, AZ 85721, USA Lunar and Planetary Laboratory, The University of Arizona, Tucson, AZ 85721, USA Max Planck Institute for Astronomy, K’́onigstuhl 17, Heidelberg, D-69117, Germany Earths in Other Solar Systems Team, NASA Nexus for Exoplanet System Science Fred J. Ciesla Department of the Geophysical Sciences, The University of Chicago, 5734 South Ellis Avenue, Chicago, IL 60637 Earths in Other Solar Systems Team, NASA Nexus for Exoplanet System Science

Parameters

Gijs D. Mulders Lunar and Planetary Laboratory, The University of Arizona, Tucson, AZ 85721, USA Earths in Other Solar Systems Team, NASA Nexus for Exoplanet System Science Ilaria Pascucci Lunar and Planetary Laboratory, The University of Arizona, Tucson, AZ 85721, USA Max Planck Institute for Astronomy, K’́onigstuhl 17, Heidelberg, D-69117, Germany Earths in Other Solar Systems Team, NASA Nexus for Exoplanet System Science Dániel Apai Department of Astronomy, The University of Arizona, Tucson, AZ 85721, USA Lunar and Planetary Laboratory, The University of Arizona, Tucson, AZ 85721, USA Max Planck Institute for Astronomy, K’́onigstuhl 17, Heidelberg, D-69117, Germany Earths in Other Solar Systems Team, NASA Nexus for Exoplanet System Science Fred J. Ciesla Department of the Geophysical Sciences, The University of Chicago, 5734 South Ellis Avenue, Chicago, IL 60637 Earths in Other Solar Systems Team, NASA Nexus for Exoplanet System Science
Abstract

The Kepler survey provides a statistical census of planetary systems out to the habitable zone. Because most planets are non-transiting, 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 multi-planet systems based on the latest Kepler data release, DR25. We estimate that at least of sun-like 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 sun-like stars of . The innermost planets in multi-planet 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 (non-transiting) planets at shorter orbital periods, hence knowledge of a close-in exoplanet could be used as a way to optimize the search for Earth-size planets in the Habitable Zone with future direct imaging missions.

planetary systems — planets and satellites: formation — protoplanetary disks — methods: statistical — surveys
\correspondingauthor

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 earth-sized planets in the habitable zone of sun-like stars are just below the detection limits of Kepler (Thompson et al., 2018), occurrence rates of earth-sized 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; Foreman-Mackey 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 sub-Neptunes 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 multi-planet system are not transiting, and modifying planetary occurrence rates from those based on the observed planetary architectures to account for non-transiting 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 multi-planet systems (Lissauer et al., 2011; Fang & Margot, 2012; Fabrycky et al., 2014; Ballard & Johnson, 2016). However, these simulations under-predict 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 sub-populations 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 111https://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 (Foreman-Mackey et al., 2013), and has already been employed in two different studies (Kopparapu et al., 2018; Pascucci et al., 2018).

Figure 1: Flowchart of the Exoplanet Population Observation Simulator. A description of all mathematical symbols can be found in table 3.

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 1-4 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 1-4 take less than a second to run on a single-core 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 multi-dimensional 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.

Multi-planet 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.

Figure 2: Example planet probability distribution function, . The color indicates the planet occurrence rate in percent per unit area of . The side panels show the marginalized distributions, in units of and in units of .

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 sub-Neptune 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 Multi-planet mode

In its multi-planet 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 multi-planet 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 multi-planet 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

Figure 3: Illustration of the planetary system architecture. denotes the orbital period of the innermost planet, while denotes the period ratio between adjacent planets. The system is inclined with respect to the observer by an angle , and each planet has a mutual inclination with respect to this system inclination. Note that each planet has a different mutual inclination and therefore a different inclination with respect to the line of sight. All orbits are assumed to be circular and the ellipiticity of the projected orbits are due to their respective inclinations. The probability distributions for , , , and are described in the text.
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 edge-on 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 log-normal distribution. The orbital period distribution of the k-th 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 equal-sized planets over randomly assigned sizes because planets in multi-planet 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 (, , , , , , , , , , ).

Figure 4: Simulated sample of transiting (pink) and detectable (blue) planets generated in occurrence rate mode, see Figure 2.

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 multi-planet 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 semi-major 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 Multi-planet 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 line-of-sight 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.

Figure 5: Planet detection efficiency for dwarf stars from the Kepler mission for the DR25 data release as function of orbital period and planet radius. Left Panel: The detection efficiency, , the probability that a planet candidate is detected based on the signal-to-noise of the transit. Middle Panel: The vetting efficiency, , the probability that a planet is classified as a reliable planet candidate (a Robovetter disposition score of ). Right Panel: The survey completeness, , which includes both the detection and vetting efficiency as well as the the geometric transit probability,

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 signal-to-noise 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 KeplerPORTs222https://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 archive333https://exoplanetarchive.ipac.caltech.edu/docs/Kepler_completeness_reliability.html.

The detection efficiencies were evaluated on a grid with 20 logarithmically-spaced orbital period bins between and days assuming planets on circular orbits, and 21 logarithmically spaced planet radius bins between and . After removing giant and sub-giant 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 products444https://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 main-sequence stars. We use a power-law in planet radius and a broken power-law in orbital period to obtain a smooth function for the vetting completeness

(17)

with , , , , (see Appendix B for details). The best-fit vetting efficiency is shown in Figure 5. The vetting efficiency varies between close to 100% for Neptune-sized planets at short orbital periods to 25% near the habitable zone. The break in the power-law is needed to match the reduced vetting efficiency for planets at orbital periods larger than days as described in Thompson et al. (2018).

Figure 6: Simulated sample of detectable planetary systems. Planets with no additional planets detected in the system are color-coded in gray. Colors indicate the number of observed planets per system.

The detectable planet sample, , typically contains about planets (Figure 4). In multi-planet mode, EPOS also keeps track of which planets are observable as part of multi-planet systems (Figure 6). From the observed multi-planet 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).

Figure 7: Period ratio distribution of adjacent planets in multi-planet systems, . The intrinsic distribution is shown with the red line while the solid blue histograms show the distribution of detections in the simulated survey. The observable distribution is skewed towards shorter orbital period ratios by detection biases. The hatched region indicates the observable planet pairs with at least one non-detected planet between them.
Figure 8: Orbital period distribution of the innermost planet in each system. The intrinsic distribution is shown with the red line while the solid blue histograms show the distribution of detections in the simulated survey. The observable distribution is skewed towards shorter orbital periods by detection biases.

It is worth noting that the properties of the observable multi-planet 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 non-transiting 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 multi-planet 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 sub-Neptunes (Dong & Zhu, 2013; Santerne et al., 2016), are less often part of multi-planet 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 photo-evaporation desert (e.g. Lundkvist et al., 2016) where the planet-radius distribution deviates significantly from that at larger orbital periods. The other bounds are chosen because there are very few planet detection outside this range.

Figure 9: Observed sample of planetary systems. Planets with no additional planets detected in the system are color-coded in gray. Colors indicate the number of observed planets per system.

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 sub-giant 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 high-reliability 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.

Figure 10: Comparison of simulated planets for the example model (blue) with detected planets (orange). The comparison region (black box) excludes hot Neptunes ( days) and giant planets ().

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 two-sample Kolmogorov-Smirnoff (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.

Figure 11: Simulated versus observed frequency of multi-planet systems. The blue histogram shows the example model with an average mutual inclination of . Multi-planet statistics from Kepler derived in the same radius and period range are shown in orange. The hatched region indicates the excess of single transiting planets, here of systems (). Crosses indicate a population of planetary systems on co-planar orbits (green, ) and on isotropic orbits (red, ).

2.4.2 Multi-planet statistics

We calculate three additional summary statistics for multi-planet systems: the frequency of multi-planet 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 2-sample 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 multi-planet frequencies of the simulated sample are drawn from the same distribution as the observations.

The frequency distribution of planets in multi-planet 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 single-transiting systems. We also show a model with only co-planar orbits ( and , green) that over-predicts the frequency of multi-planet systems, particularly at high numbers of planets per system. A model with planets on isotropic orbits (, red) under-predict the frequency of all multi-planet systems.

Figure 12: Simulated (blue) versus observed (orange) period ratio between adjacent planets in multi-planet systems. The red line shows, for comparison, a simulation where planet orbits are randomly drawn from the period-radius distribution (Figure 2) instead of regularly spaced.
Figure 13: Simulated (blue) versus observed (orange) location of the innermost planet in each system.

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 period-radius 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 multi-planet systems is shown in Figure 13. We only focus on the innermost planet of a detected multi-planet system, as for observed single-planet systems it can not be derived from the observable if they are intrinsically single or part of a multi-planet system with non-transiting 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 1-4 is less than a tenth of a second in occurrence rate mode and less than a second in multi-planet mode, allowing for an comprehensive sampling of parameter space. We use emcee (Foreman-Mackey et al., 2013), an open-source 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 7-dimensional 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. .

Figure 14: Period-radius distribution of detected planets in occurrence rate mode. The green population shows the population generated from the best-fit parameters. The black box and dashed lines indicate the range of orbital periods and planet sizes that is included in the observational comparison. The side histograms show the marginalized simulated distributions compared with the observed populations in orange. The blue lines show 30 samples from the posterior.
Parameter Example Best-fit
(days)
()
Table 1: Fit parameters for the example model and the best-fit solutions with confidence intervals.

We run emcee with 200 walkers for 5,000 iterations, allowing for a 1000-step burn-in 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 best-fit values and their confidence intervals are calculated as the , and and percentiles, respectively. The simulated model for the best-fit parameters and for 30 samples of the posterior is shown in Figure 14.

Figure 15: Posterior orbital period distribution (top) and planet radius distribution (bottom). The red bars show the occurrence rates estimated using the inverse detection efficiencies for comparison. Note that the occurrence rates underestimate the true distribution in bins that include regions where Kepler has not detected any planet candidates, in particular and days (see Figure 16).
Figure 16: Kepler DR25 candidate list, color coded by survey completeness. The sample includes only dwarf stars () and planet candidates with a disposition score larger than 0.9. The planet occurrence rate estimated from Eq. 19 is .

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 Best-fit
(days)
()
()
Table 2: Fit parameters for the example model in multi-planet mode and the best-fit solutions with confidence intervals. Parameters , , and were fixed to their best-fit solutions from occurrence rate mode.
Figure 17: Posterior predictive plots of simulated planetary systems. 30 samples from the posterior are shown in blue and the observed systems are shown in orange. The top left panel show the best-fit period-radius distribution in green, see Figure 14. The black box and dashed lines indicate the range of orbital periods and planet sizes that is included in the observational comparison. The top right panel shows the frequency distribution of multi-planet systems with planets. The bottom left panel show the period ratio distribution of adjacent planets. The bottom right panel shows the orbital period distribution of the innermost planet in each multi-planet system.

2.6 Multi-planet systems

In multi-planet mode, we explore the 9-dimensional 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 multi-planet 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 best-fit 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 500-step burn-in. The MCMC chain and parameter covariances are shown in the appendix (Figure 24). The simulated observables of the best-fit 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 SAG13555https://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 regularly-spaced 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 co-planar 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 single-planet transiting systems, which we model as multi-planet 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 pre-main-sequence sun-like stars (Millan-Gabet 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 sub-Neptunes 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.

Figure 18: Marginalized orbital period distribution of the innermost planet in the system. The green line indicates the distribution of all planets in the system.

4 Discussion

4.1 How rare is the solar system?

Our modeling analysis indicates that multi-planet 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 closer-in 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 best-fit model. Such a simulated planet population with and matches planet occurrence rates exterior to days for days. The simulated observations predict multi-planet 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 co-planar 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 multi-planet 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 multi-planet systems without planets interior to Venus, though we do not find evidence that such a population exists.

Figure 19: Best-fit distribution of planetary system properties from Kepler (blue) compared to the solar system terrestrial planets (red letters). The inclinations of the terrestrial planets are with respect to the invariable plane. The orbital periods, radii, inclinations, and period ratios of the terrestrial planets lie near the peak of the distribution of kepler systems. The only notable exception is the orbital period of the innermost planet, where Mercury (93th percentile) and Venus (97th percentile) lie in the tail of the distribution.

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 super-earths/mini-Neptunes in the solar system is also notable (Martin & Livio, 2015), earth-sized 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 earth-sized planets with earth-like orbital periods (“habitable zone” planets), here defined as and based on the Kopparapu et al. (2013) conservative habitable zone for a sun-like 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 literature666https://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 over-estimating 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 best-fit 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 Q1-Q16 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 multi-planet systems where the innermost planets are typically at an orbital period of days. While the multi-planet statistics do not directly constrain planetary systems of earth-sized planets in the habitable zone of sun-like 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 multi-planet systems are not transiting, the probability of finding any planet (transiting or non-transiting) is even higher. For example, the detection of an earth-mass 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 single-transiting planets are intrinsically single or highly-inclined multiples, of systems with close-in 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.1-1 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).

Figure 20: Illustration of how the estimated fraction of stars with planetary systems depends on the distribution of planets across stars. Nearly-coplanar systems orbit of stars (green) while the distribution of observed single transiting planets (purple) is less well constrained.

4.4 Do most stars have planets?

Our analysis confirms that the average number of planets per (sun-like) 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 multi-planet systems. An important assumption we made is that the apparent excess of single transits is because of multi-planet 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 highly-inclined planets in otherwise coplanar multi-planet 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 multi-planet 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 single-transit systems (Morton & Winn, 2014), indicating that they may have highly inclined orbits. Several mechanisms have been proposed to disrupt initially co-planar 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 host-star 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 multi-planet systems and have non-transiting or non-detected planets in the system, and hence the sample of “singles” is diluted with multi-planets, weakening any potential trends. Transit timing variations can also be used to detect additional non-transiting 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 earth-size 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 sub-Neptunes 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 long-period 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 sun-like stars (Cumming et al., 2008), most of them beyond 1 au. Microlensing surveys hint at the existence of a population of Neptune-mass 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 multi-planet frequencies are consistent with ensemble populations that are a mix of nearly co-planar 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 sun-like stars. The remaining of exoplanets could be intrinsically single planets or be part of multi-planet 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.

We thank an anonymous referee for a constructive review of the manuscript. We also thank Shannon Dulz and Peter Plavchan for feedback on an early manuscript draft, Ed Bedrick for advice on statistical methods, and Daniel Huber for advice regarding stellar properties. We acknowledge helpful conversations with Carsten Dominik, Eric Ford, Daniel Carrera, Renu Malhotra, and Rachel Fernandes. This material is based upon work supported by the National Aeronautics and Space Administration under Agreement No. NNX15AD94G for the program “Earths in Other Solar Systems”. The results reported herein benefited from collaborations and/or information exchange within NASA’s Nexus for Exoplanet System Science (NExSS) research coordination network sponsored by NASA’s Science Mission Directorate. \software NumPy (van der Walt et al., 2011) SciPy (Jones et al., 2001–) Matplotlib (Hunter, 2007) Astropy (The Astropy Collaboration et al., 2018) EPOS (Mulders, 2018) emcee (Foreman-Mackey et al., 2013) corner (Foreman-Mackey, 2016) KeplerPORTs (Burke & Catanzarite, 2017)

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, KSCI-19114-001
  • 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)
  • Foreman-Mackey (2016) Foreman-Mackey, D. 2016, The Journal of Open Source Software, 24, doi:10.21105/joss.00024
  • Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, Publications of the Astronomical Society of Pacific, 125, 306
  • Foreman-Mackey et al. (2014) Foreman-Mackey, 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
  • Millan-Gabet et al. (2007) Millan-Gabet, 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., & Molenda-Zakowicz, 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, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, ArXiv e-prints, 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 k-th planet in system
Power law indices
Normalization factor
Dimensionless spacing parameters
Table 3: List of mathematical symbols used in this paper.
Figure 21: The top left panel show the vetting efficiency per grid cell for a Robovetter score larger than . The side panels show marginalized detection efficiency (solid purple) compared to the best-fit double broken power-law model (green dashed). The top right panel shows the locations of the planet injections from Christiansen (2017). The bottom left panel shows the best-fit detection efficiency evaluated per grid cell. The bottom right panel shows the difference between the measured and fitted vetting completeness per grid cell.

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 signal-to-noise ratio, vetting is a trade-off 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 page777https://exoplanetarchive.ipac.caltech.edu/docs/KeplerSimulated.html. We removed all injections not on main-sequence 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 sub-Neptunes 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 power-law:

(B1)

with

(B2)

and

(B3)

We find the best-fit parameters using scipy.optimize.curve_fit, where we minimize the difference between the measured and parametrized vetting completeness over all grid cells. The best-fit solution is shown in the bottom left panel of Figure 21 and the bottom right panel shows the residuals per grid cell. The best-fit parameters are , days, , , , , and . The break in the power-law 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 power-law 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 power-law index radius break in EPOS. Instead, we assume the planet radius dependence of the detection is described by the power-law in planet radius that was fitted to the small planets. This yields a vetting efficiency of

(B4)

with , which is shown in Figure 5.

Appendix C Posterior Distributions

The MCMC chain and corner plots for occurrence rate mode are shown in Fig. 22 and Fig. 23. The corner plots for multi-planet mode are shown in Fig. 24

Figure 22: Position in 7-dimensional parameter space of 200 walkers in the MCMC chain. The dashed line indicates the end of the burn-in phase.
Figure 23: Corner plots for the parametric distribution, generated using open-source Python package corner. Blue lines indicate the initial guess that was based on previous studies.
Figure 24: Corner plots for the multi-planet parametric distribution, generated using open-source Python package corner. Blue lines indicate the initial guess.
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
""
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
   
Add comment
Cancel
Loading ...
200009
This is a comment super asjknd jkasnjk adsnkj
Upvote
Downvote
""
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters
Submit
Cancel

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test
Test description