The Exoplanet Population Observation Simulator. II -- Population Synthesis in the Era of Kepler

TL;DR

This paper introduces epos, a simulator that compares planet formation models with Kepler data, providing benchmarks and insights into physical processes shaping planetary systems.

Contribution

It demonstrates how epos constrains planet formation models by matching Bern synthesis models to Kepler observations and offers diagnostics for future modeling efforts.

Findings

A model with 100 seed cores per disk fits most diagnostics.

Moving the planet trap outward improves model accuracy.

Increasing mini-Neptune formation efficiency enhances model fit.

Abstract

The collection of planetary system properties derived from large surveys such as Kepler provides critical constraints on planet formation and evolution. These constraints can only be applied to planet formation models, however, if the observational biases and selection effects are properly accounted for. Here we show how epos, the Exoplanet Population Observation Simulator, can be used to constrain planet formation models by comparing the Bern planet population synthesis models to the Kepler exoplanetary systems. We compile a series of diagnostics, based on occurrence rates of different classes of planets and the architectures of multi-planet systems, that can be used as benchmarks for future and current modeling efforts. Overall, we find that a model with 100 seed planetary cores per protoplanetary disk provides a reasonable match to most diagnostics. Based on these diagnostics we identify physical properties and processes that would result in the Bern model more closely matching the known planetary systems. These are: moving the planet trap at the inner disk edge outward; increasing the formation efficiency of mini-Neptunes; and reducing the fraction of stars that form observable planets. We conclude with an outlook on the composition of planets in the habitable zone, and highlight that the majority of simulated planets smaller than 1.7 Earth radii have substantial hydrogen atmospheres. The software used in this paper is available online for public scrutiny at https://github.com/GijsMulders/epos

Full text

The software used in this paper is available online for public scrutiny at \faGithub https://github.com/GijsMulders/epos.

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

††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, 2018b) KeplerPORTs (Burke & Catanzarite, 2017)

1 Introduction

Exoplanet surveys have revealed a diverse population of planets with a wide range of sizes, orbital period, and system architectures. Their origins, however, are not well constrained observationally. While follow-up observations can provide important clues to the chemical composition of planets and their atmospheres, such observations are limited to a small and biased set of predominantly large and hot exoplanets. Thus, studying the formation of planets and their emergent ensemble properties are crucial for understanding the origins of the entire exoplanet population (e.g. 2019BAAS...51c.475A).

The diversity in exoplanet properties likely results in part from a range of initial protoplanetary disk properties and evolutionary pathways. Further, the late stages of planet formation are inherently stochastic, with simulations with similar initial conditions producing a wide range of planets and planetary systems (e.g. Raymond et al., 2009; Fischer & Ciesla, 2014). Thus linking the properties of the observed exoplanet population to a formation history can be quite challenging and requires a robust statistical analysis comparing model outcomes to the available data.

Population synthesis models have been used to understand how the interplay between initial disk properties and various physical processes shapes the planetary systems that we ultimately observe (see Mordasini 2018 for a recent review). The processes accounted for in these models include, but are not limited to, the evolution and dispersal of protoplanetary disks, the growth of planetary cores, the accretion and loss of gaseous envelopes, the migration of planets, and gravitational interactions between protoplanets. These models then make quantitative predictions about the physical properties and chemical compositions of planets, and could ultimately be used to asses the relative importance of different planet formation channels.

For example, the majority of exoplanets found by Kepler have substantial hydrogen atmospheres (e.g. Wu & Lithwick, 2013) that must have accreted in a gas-rich environment — in contrast with the gas-poor assembly of the solar system terrestrial planets. The planet population identified by Kepler extend close to the habitable zones of Sun-like stars and can be used to estimate $\eta_{\oplus}$, the fraction of stars with Earth-sized planets in the habitable zone. Studies that extrapolate the planet occurrence rates from regions where those statistics are reliable typically find $\eta_{\oplus}\sim 30\%$ (Burke et al., 2015; Mulders et al., 2018). Studies that estimate $\eta_{\oplus}$ directly from planet candidates in the habitable zone find that a wide range of rates is still possible due to the low number of reliable detections (Hsu et al., 2019).

If exo-Earths do indeed follow different formation pathways than our Solar System, then it is likely that their bulk compositions would vary significantly from those found around the Sun. Indeed, the need for different formation pathways is seen in looking at the planets around M dwarfs; Traditional Solar System-based models predict that large numbers of low-mass planets would be found around such stars (Raymond et al., 2007; Ciesla et al., 2015) while surveys indicate that such systems contain much more massive planets, more so than around the higher mass stars (Dressing & Charbonneau, 2015; Mulders et al., 2015a). The TRAPPIST-1 system (Gillon et al., 2017) is a clear example as it hosts 7 Earth-mass planets inside 0.07 au despite being  10% the mass of the Sun. These planets are thought to have water-rich compositions (Unterborn et al., 2018), possibly indicating that inward migration was a critical factor in their formation (e.g. Ormel et al., 2017).

While the importance of different formation pathways for exoplanets is evident, the relative importance of various processes in defining these pathways of the observed exoplanet population is less clear. Planet population synthesis models make quantitative predictions about the occurrence rates of planets and planetary systems as different processes are introduced to the simulations. Hence, statistical comparisons of model results with observed exoplanet populations can be used to calibrate the relative importance of these formation channels.

Early planet population synthesis models were mainly focussed on the giant planet populations that were known at the time (Ida & Lin, 2004), and often simulated the formation and growth of one planet per disk (e.g. Mordasini et al., 2009a). These models predicted a dearth of close-in super-earths, but such planets were later found to be abundant in more sensitive exoplanet surveys (Howard et al., 2010, 2012). With the inclusion of multi-planet interactions and type-I migration, the Bern planet population synthesis models (Mordasini 2018; Emsenhuber et al. in prep) are now capable of producing planets across the entire range in planet size and orbital period where they are observed with Kepler (Alibert et al., 2013; Jin & Mordasini, 2018). A crucial next step is to evaluate whether the predicted planet populations also match the statistical properties of observed exoplanets.

In this paper, we use epos, the Exoplanet Population Observation Simulator (Mulders et al., 2018), to evaluate the ability of the latest Bern planet population synthesis model (Mordasini, 2018) to produce planets and planetary systems that are consistent with the Kepler exoplanet population. We assess the relative success of the Bern model in reproducing current exoplanet statistics, and identify areas where the predicted planet populations fall short of doing so. We conclude by identifying parameters and processes that should be investigated further that likely would result in the models producing planets and system architectures more consistent with the known exoplanetary systems.

This study serves as a template for future comparisons between planet formation models and exoplanet surveys. As such, the software to produce the results and figures in this paper is publicly available as a Python repository at https://github.com/GijsMulders/epos or with pip install epospy.

2 The Exoplanet Population Observation Simulator

The Exoplanet Population Observation Simulator, epos 111https://github.com/GijsMulders/epos, is a software package designed to simulate survey observations of synthetic exoplanet populations. By taking into account detection biases and survey completeness, models of the planet population can be constrained from or compared to an observed population of exoplanets (Figure. 1)

In Mulders et al. 2018 (hereafter Paper I) we demonstrated the basic approach of epos v1.0 (Mulders, 2018a). In that paper we used parametric distributions of planet properties to test the methodology and make inferences about the underlying population of Kepler exoplanets using emcee (Foreman-Mackey et al., 2013). In this paper, we will expand epos with simulated planetary systems from the Bern planet population synthesis model. In this section, we will briefly describe how epos works and how we implement simulated observations of planet formation model outcomes.

2.1 Basic Approach

The Exoplanet Population Observation Simulator is a software package designed to take into account detection biases in exoplanet surveys. The basic approach has been outlined in Paper I and summarized here, consisting of 5 steps.

Step 1

Define a distribution of planetary systems. Paper I used analytic forms for planet occurrence rates and planetary architectures. Here, we will use the orbital architectures of the Bern planet population synthesis model.

Step 2

From this derive a transiting planet population using a Monte Carlo approach by assigning a random orientation 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 the geometric transit probability and survey completeness.

Step 4

Compare the detectable planet population with exoplanet survey data using a set of summary statistics of planetary system properties.

Step 5

In Paper I we repeated steps 1-4 until the simulated detectable planet population matches the observed planet population to constrain the intrinsic distribution of planetary systems. Instead, we discuss what parameters and processes may need to be further investigated to yield better matches between modeled and real exoplanet systems.

We do not specifically take into account variations in planet occurrence with stellar properties such as mass, metallicity, and potentially age. Instead, we assume that variations in planet properties will average out across the survey, and that the observed planetary system properties are representative of a $5$ Gyr old, solar mass, solar metallicity star.

2.2 Implementation of Planet Formation Mode

We implement a “planet formation mode” into version v2.0 (Mulders, 2019) of epos 222https://github.com/GijsMulders/epos to simulate observations of a planet formation model. We adjust steps 1,4, and 5 in as follows.

In step 1, we simulate a survey of $n_{\star}$ stars, where a fraction $\eta_{s}$ has a planetary system. Each planetary system is randomly drawn from a population synthesis model, and has properties $\{R_{p},P,\delta i\}_{k}$, where: $R_{p}$ is the planet radius; $P$ is the orbital period; $\delta i$ is the inclination of the planet with respect to the system inclination $i$; and $k$ is an index for each planet in the system counted from inside out. The system inclination, $i$, and longitude of ascending node for each planet’s orbit $\Omega_{k}$, are generated using a Monte Carlo approach as described in Paper I.

We then carry out a synthetic survey by determining which planets are transiting and detectable as in steps 2 and 3 to generate an observable distribution. We describe the comparison with observations, the equivalent of step 4, in Sections 4 and 5. We do not perform the MCMC fitting (step 5) to refine parameters in the models, but highlight which features and processes could be modified to improve the match between the synthetic surveys and that of Kepler.

2.3 Kepler Survey Completeness

The survey completeness is defined as the probability that the Kepler pipeline detects a transiting planet of the same size and orbital period and classifies it as a reliable planet candidate (e.g. vetting). The detection efficiency of the Kepler pipeline for each star is calculated with KeplerPORTs (Burke & Catanzarite, 2017). We use the stellar radii and giant star classification from Gaia (Berger et al., 2018) which yields ($n=122,682$) main sequence stars, mainly of spectral types F, G, and K. The survey detection efficiency is then the average detection efficiency of all main sequence stars. While it is possible to select a narrower spectral type range to select the stars most similar to the sun, in practice this does not affect the survey detection efficiency because the average star in the Kepler survey is already of similar mass and radius to the sun. To optimize the number of detected planets in the sample, we choose the larger stellar sample, which is representative of sun-like stars.

The planet candidate vetting efficiency is calculated based on the procedure in Thompson et al. (2018) assuming a planet candidate is a reliable detection if the disposition score is above 0.9. We parametrize the vetting efficiency as a double broken power-law in period and radius as described in Paper I.

We do not explicitly account for the effects of stellar multiplicity or astrophysical false positives. The Kepler planet host sample has received extensive follow-up observations and analysis (Furlan et al., 2017, 2018). As a result the false positive probability of planet candidates in the DR25 catalog is below $1\%$ (Morton et al., 2016) and the impact of stellar multiplicity on planet occurrence rates is minimal (e.g. Ziegler et al., 2018).

2.4 Parametric Distributions

Here we summarize the parametric distributions from Paper I that we use for the comparison with the Bern planet population synthesis model in Sections 3 and to guide the interpretation of the simulated observables in Section 5.

The parametric distributions define five independent distributions of planet properties: the orbital period, $P$; the planet Radius, $R_{p}$; the mutual inclination333We refer to the inclination between a planet and the system inclination as mutual inclination to discriminate between the inclination of a planet with respect to the observer. of planets in multi-planet systems, $\delta i$; the period ratio between planets, $\mathcal{P}$; and the period of the innermost planet, $P_{0}$. Each distribution was assumed to be independent of the other.

• •

The orbital period distribution is a broken power-law with a break at $P\sim 12$ days. Interior to the break, the occurrence rate of sub-Neptunes rapidly decreases.

• •

The planet radius distribution is also a broken power-law with a break at $R_{p}\sim 3R_{\oplus}$. Above the break, planets become rarer with increasing size.

• •

The mutual inclination distribution is a Rayleigh distribution with a mode at $\delta i=2\arcdeg$.

• •

The period ratio distribution is generated from a log-normal distribution in dimensionless spacing. The median period ratio is $\mathcal{P}=2.1$

• •

The period distribution of the innermost planet is a peaked distribution described by a broken power-law with a break at $P\sim 12$ days.

The purpose of the comparisons with parametric distributions is to develop an intuitive understanding of the major differences between the model outcomes and the intrinsic planet population that may not be immediately apparent from the direct observables. For example, the mutual inclination distribution is not a direct observable, but has to be derived from multi-planet frequencies using a set of assumptions about planet architectures. Of course, some of those assumptions may not be justified by the observations or the planet formation model, and may skew the derived parameters. For example, as Zhu et al. (2018) have recently pointed out, a different set of assumptions may lead to very different conclusions on the mutual inclination distribution. A similar case can be made for the period and radius distributions, which are predicted to be covariant due to the sculpting effect of photo-evaporation and which may bias estimates of $\eta_{\oplus}$ (Lopez & Rice, 2018; Pascucci et al., 2019). Therefore, we will use the parametric distributions for an initial comparison with the models but also perform a synthetic survey to make a direct comparison between the simulated systems and the Kepler exoplanets in Section 5.

3 Population Synthesis Model

In this section we briefly describe the Bern planet population synthesis models. We use an updated version of the planet population synthesis model described in Mordasini (2018) for a solar-mass host star (Emsenhuber et al. in prep). This global model of planet formation and evolution includes the accretion of planetesimals by seed planetary cores, the gravitational interactions between planet cores using the Mercury integrator (Chambers, 1999), gas-driven planet migration, the accretion of gaseous envelopes from an evolving disk, and atmospheric loss through photo-evaporation and impact stripping (Alibert et al., 2005; Fouchet et al., 2012; Mordasini et al., 2012b, a; Fortier et al., 2013; Jin et al., 2014; Benz et al., 2014; Mordasini et al., 2015). We investigate four different sets of models, each with $1$, $20$, $50$, and $100$ seed planetary cores per disk. The models with one core per disk are mainly included for comparison with previous studies (e.g. Mordasini et al., 2009a, b), and can not be directly compared to the multi-planet systems from Kepler.

Each of the four sets consists of an ensemble of stars with a planet-forming disk. All stars are assumed to be single and of solar mass, with a composition following the observed metallicity distribution of field stars. Every disk contains a reservoir of solids in the form of $300$ meter-sized planetesimals. The total mass in planetesimals is calculated from the disk mass and stellar metallicity. A range of initial disk conditions is used (see Mordasini 2018) to simulate the distribution of observed protoplanetary disk properties (e.g. Andrews et al., 2013; Pascucci et al., 2016). The distribution of disk masses spans two orders of magnitude to reflect the observed range in young protoplanetary disk masses (Tychoniec et al., 2018). The total mass in planetesimals per disk ranges between $6$ and $1,500\,M_{\oplus}$, with a mean of $165\,M_{\oplus}$.

Each disk is also seeded with planetary cores that can accrete those planetesimals and collide with other cores to grow into planets. The seed planetary cores are distributed randomly with a uniform probability in log of the distance between $0.06$ and $40$ au, with $1$, $20$, $50$, or $100$ planetary cores per disk. The starting masses of the seed planetary cores are $0.01~{}M_{\oplus}$. In all models, the seed planetary cores constitute a negligible fraction of the solid mass compared to the planetesimal disk mass.

N-body interactions are calculated for $10$ Myr in the 20 and 50-seed core models and for $20$ Myr for the 100 seed core model. This timescale is long enough to resolve system-wide instabilities after disk dispersal and form dynamically stable planetary systems within $1$ au (e.g. Hansen & Murray, 2012). After that, thermodynamic evolution of the atmosphere (cooling and contraction) and atmospheric escape are tracked for $5$ Gyr. Planet radii are calculated self-consistently by solving the internal structure equations both during the formation and evolution phase, based on atmospheric accretion, loss, and cooling (e.g. Mordasini et al., 2012b, a; Jin & Mordasini, 2018). The largest changes in planetary radii occur during the first $\sim 100$ Myr, and the model does not predict significant radius evolution for the age range of stars in the Kepler field.

A thousand systems are simulated in each set of multi-seed models. Each set of the $20$, $50$, and $100$ seed planetary core models yields $983$, $954$, and $978$ surviving planetary systems, respectively. Of each planetary system we use the planet radii, orbital periods, and inclinations ($\{R_{p},P,\delta i\}_{k}$) as input to epos. 30,000 planets are simulated in the single seed core models that yield 29,993 planets. For the single seed model we use the set of planet radii and orbital periods ($\{R_{p},P\}$) as input to epos.

To show the typical type of output of the Bern population synthesis models, the ensemble of orbital periods and planet radii of simulated exoplanets in the 1 and 100 seed-core model are shown in Figure 2. The distribution of planetary system properties of 100 seed-core model are shown in Figure 3.

The simulated ensemble of exoplanets of the 100 seed core model (Fig. 2a) populate the same parameter space as the Kepler exoplanets, ranging in orbital period from a day to a year and in planet size from Mars to Jupiter. A quick comparison between the simulated exoplanets and the parametric distributions derived from Kepler in Paper I show that the model under-predicts the occurrence of planets with sizes between 2-3 $R_{\oplus}$ (mini-Neptunes) and over-predicts the abundance of (rocky) planets at orbital periods less than $\sim{}10$ days.

The model with one seed planetary core per disk shows more structure in the period-radius distribution (Fig. 2b). The radius distribution shows a pronounced radius valley at $\sim 1.6R_{\oplus}$ as seen in the Kepler data (Fulton et al., 2017; Jin & Mordasini, 2018). The period-radius distribution, however, also shows many features not present in the Kepler data: In particular, the model shows a dearth of short-period earth-sized planets which are below the Type-I migration mass; a spike of mini-Neptunes and giants at an orbital period of 2 days corresponding to the inner edge of the computational disk at 0.03 au; several pile-ups (at about 20, 50 and 100 days) caused by different types of migration traps in the disk (Dittkrist et al., 2014); and a spike of Jupiter-sized planets at all periods resulting from runaway gas accretion. Otherwise, the planet radius distribution shows a similar trend as Kepler with a fairly constant occurrence of planets smaller than $3R_{\oplus}$ and a sharp decrease in occurrence for larger planets.

Figure 3 show the typical outcome of planetary system architectures of the Bern model, here for a model with 100-seed cores. Figure 3a shows the orbital inclination distribution of planets versus their semi-major axis. The inclinations show no clear trend with semi-major axis, and the marginalized inclination distribution has the same shape as the Rayleigh distribution used in the parametric solution of Paper I. Figure 3b shows the period ratios, $\mathcal{P}$, of adjacent planet pairs as a function of semi-major axis. It also shows the location of the innermost planet in each planetary system. The simulated period ratio distribution becomes wider at larger semi-major axes, but there is no clear sign of a covariance between period ratio and star-planet separation. There are concentrations of planets at the 2:1 and 3:2 orbital resonances.

The distributions of inclinations, period ratios, and innermost planet locations of the Bern model follow wide but peaked distributions that are similar in shape to the distributions observed or inferred from Kepler. However, there are also a number of key differences between the observations and the simulated systems in the 20, 50, and 100 core models that we will discuss in depth in the next sections.

4 Planet Occurrence Rates

The next step is to compare the model’s planet occurrence rates for different classes of planets to the Kepler data. The occurrence rate in a certain period-radius area is defined as the total number of simulated planets in that area across all simulations divided by the number of simulated stars, multiplied by a factor $\eta$, the fraction of stars that have planetary systems akin to those modeled here. The factor $\eta=0.2$ is chosen such that the total number of detectable planets is roughly equal to that of Kepler. For the simulations with one seed planetary core per disk, we set $\eta=2$ because stars have, on average, more than one planet.

In Figure 4 we show the simulated populations as they would appear when observed in the Kepler survey. Figure 4e shows the exoplanet candidate list color-coded by the survey completeness. The planet occurrence rates are then calculated as the sum of the inverse completeness of each planet in the area divided by the number of stars in the survey, $n_{\star}=122,682$.

Here, we define occurrence rates of different classes of planets that roughly compare to the observed clusters in the Kepler data, Fig. 4e:

• •

Hot Jupiters: $R=5-20~{}R_{\oplus}$, $P=1-10$ days.

• •

Warm Giants: $R=4-14~{}R_{\oplus}$, $P=20-300$ days.

• •

Hot Earths: $R=0.7-2.2~{}R_{\oplus}$, $P=3-200$ days, and below the radius valley.

• •

Mini-Neptunes: $R=1.2-4.0~{}R_{\oplus}$, $P=3-200$ days, and above the radius valley.

The period and radius ranges are chosen to span roughly equal areas in log space, with $d\ln Pd\ln R_{p}\approx 3.1$. Hot and warm giant planets are separated by the “period valley” at $\sim 10$ days (Santerne et al., 2016). Giant planets are separated from sub-Neptunes at $4~{}R_{\oplus}$. Sub-Neptunes are split into Hot Earths and mini-Neptunes at the evaporation valley (Fulton et al., 2017). Bins are mostly bounded by detection limits at the largest periods and the smallest radii where planets have been detected. Most notably these bins exclude the habitable zone, though we will make inferences based on the simulated planet populations in Section 6.

The Kepler planet occurrence rates indicate that mini-Neptunes are the most populous class of planets at $62.3\%$, followed by Hot Earths ($42\%$), Warm Giants ($13\%$), and last and least Hot Jupiters ($0.6\%$). The confidence intervals are based on poisson statistics only, and do not account for systematic uncertainties in the planet occurrence rate calculation.

Planet occurrence rates for the 20, 50, and 100 seed models (Fig. 4bcd) are calculated assuming 1 in 5 stars have a simulated planetary system ($\eta=0.2$). Despite the different number of seed cores, the overall planet occurrence rates as function of orbital period and planet radius are very similar. This indicates that the number of planets formed in each simulation is not dependent on the amount of seed cores but on other properties that do no vary between the different sets of simulation, such as the mass/amount of planetesimals that can be accreted. All three multi-seed models roughly produce the same planet occurrence rates for Hot Jupiters (1-1.5%), Warm Giants (8-11%), mini-Neptunes (13-14%), and Hot Earths (51-60%). The simulated rates approach the derived planet occurrence rates from Kepler within a factor of 2, with the exception of mini-Neptunes which the models under-predicts by a factor 4-5.

The single seed model predicts a similar occurrence rate of Warm Giants as the multi-seed models and as the Kepler data when $\eta=2$ planets each are assigned to the simulated stars. The occurrence rate of mini-Neptunes ($36\%$) is thrice that of the multi-seed models but still twice as low as Kepler. However, the single seed model over-predicts the number of hot Jupiters by more than a factor 10. The model also over-predicts the number of Hot Earths by a factor $2$ at $84\%$. It should be noted that the clustering of rocky planets in period-radius space is very different from the observed smooth distribution of Kepler planets. The occurrence rate of hot Earths is therefore very sensitive to the locations of the bin boundaries in Figure 4.

We also show the planet occurrence rates as a function of planet radius and orbital period in Figure 5. The giant planet occurrence rates in the multi seed-core models (Fig. 5a) show the same increasing trend with orbital period as the Kepler planet occurrence rates (see also Fernandes et al. 2019). The single seed models show a pronounced peak for hot-Jupiters at $P\approx 3$ days, corresponding to a very strong trap for migrating planets at the disk inner edge, that is not seen in the Kepler data. Sub-Neptunes show a flat trend with orbital period, and do not capture the observed decrease in occurrence interior to 10 days (Youdin, 2011; Mulders et al., 2015b) which is likely due to a migration trap at the inner edge of the protoplanetary disk (Lee & Chiang, 2017; Carrera et al., 2019). This indicates that the planet trap in the model is located to close to the star and possibly that too many planets are pushed further in when captured in resonances with other migrating planets, as no decrease at the shortest periods is seen either.

The planet radius distribution shows the same overall trend as the Kepler distribution, with sub-Neptunes occurring roughly ten times as frequent as giant planets (Fig. 5b). The radius distribution among sub-Neptunes, however, has a distinctly different shape than the observed distribution. The multi-seed models predict a much higher fraction of earth-sized planets compared to mini-Neptunes ($2-3R_{\oplus}$) than observed. The single-core model produces a larger amount of mini-Neptunes that fits better to the observed distribution. But it also shows a pronounced peak of giant planets at $\sim 10\,R_{\oplus}$ that is not seen in the Kepler data. The observed break in the radius distribution at $\sim 3R_{\oplus}$ is reproduced by the single seed model, while the multi-seed models show a more gradual decrease in occurrence in the $1-4R_{\oplus}$ range.

The dearth of mini-Neptunes in the simulated multi-seed models highlights a key component of the new Bern planet population synthesis models. As proto-planets dynamically interact and compete for material, their growth process diverges from that of an isolated planet. While multi-planet interactions result in fewer planets undergoing type-II migration and becoming hot Jupiters, it also suppresses the growth of mini-Neptunes to a level that appears inconsistent with the observed abundances. The ratio of mini-Neptunes to rocky planets is larger in the single-seed models compared to the multi-seed models (Fig. 5b). Increasing the gas/dust mass of the disk may not be a solution, as this may trigger more planets to reach runaway growth and form giant planets, potentially overshooting the amount of warm giants and/or hot Jupiters. A more efficient formation pathway for accreting hydrogen atmospheres without ballooning into giant planets appears necessary. This was also pointed out by Suzuki et al. (2018) for planets at larger separations probed by microlensing surveys.

The evaporation valley also appears to get weaker when including more seed planetary cores, either through atmospheric stripping or through dynamical scattering, providing a potential diagnostic for calibrating the amount of dynamical interaction.

5 Planetary System Architectures

Here we inspect the architectures of multi-planet systems formed in the 20, 50, and 100-seed core models and evaluate the ability of the Bern model to reproduce planetary system architectures observed with Kepler. We make two types of comparisons between the model and the data. One method is to put the planetary system properties in context of the debiased parametric distributions derived from Kepler. The other method is to simulate a survey using the planetary systems from the Bern model and compare its observable properties directly to the Kepler data. Both methods have their pros and cons, and we will therefore present them side by side.

The parametric distributions are more intuitive as they directly tell you what the distribution of planet formation model outcomes should be. The major downsides of the parametric models of the exoplanet population, however, is that they make a number of simplifying assumptions that may not be justified by real planetary systems or by the formation model. In particular, the parametric distributions derived in Paper I do not allow for correlations between different parameters. The forward modeling approach is also better for recognizing and identifying clusters and valleys or discontinuities in distributions that are not specifically included in the parametric distributions, such as orbital resonances.

We focus on three properties of the model planetary systems and their associated observable summary statistic.

The inclination distribution of planets within a planetary system (e.g. their mutual inclinations) whose observable is the frequency of multi-planet systems. 2. 2.

The period ratio of adjacent planets in a system, whose (biased) distribution can be observed. 3. 3.

The location of the innermost planet in a planetary system, whose (biased) distribution can also be observed.

An example of these summary statistics for the 100 seed-core model and how they relate to the instrinsic model properties is shown in the appendix. All summary statistics are calculated for orbital periods between $1$–$400$ days and planet radii between $0.5$–$20~{}R_{\oplus}$. These summary statistics mainly trace the dynamical interactions between forming protoplanets and the trapping of migrating planets at the inner edge of the protoplanetary disk.

5.1 Mutual inclinations

The mutual inclinations between planets in multi-planet systems trace the amount of dynamical interactions during the late stages of planet formation. Orbital damping by gas keeps the planets’ mutual inclinations low, while gravitational interactions between (proto)planets tend to increase their mutual inclinations. The Kepler systems, with a typical mutual inclination of a few degrees (e.g. Fabrycky et al., 2014), show evidence of both damping and excitation (e.g. Dawson et al., 2016).

The inclination distributions of the 20, 50, and 100 seed core models are broad but with a clear peak (Fig. 6a). The distribution narrows when a larger number of cores are included in the model, and hence more dynamical interactions take place. The median inclination also increases from $\delta i=0.05\arcdeg$ in the 20 core model to $\delta i=0.4\arcdeg$ in the 50 core model to $\delta i=3\arcdeg$ in the 100 core-model. The inclination distribution of the 100-core model has the same shape as the Rayleigh distribution used in the parametric solution, while the distributions with fewer cores are wider and have a tail towards very low inclinations of $10^{-4}~{}\arcdeg$. These extremely low inclinations indicate that some protoplanets may not experience gravitational interactions if the number of seed cores in the simulation is low.

The median inclination of the 100 core model is close to that derived from Kepler, while the inclinations of the 20 and 50 core models are significantly lower. This is reflected in the simulated observations of the multi-seed models (Fig. 6b) that predict fewer multi-planet systems when more cores are included. The simulated observations of the multi-seed models, however, all predict a larger number of systems with two or more planets per star compared to what is observed with Kepler. The relative frequencies of 2/3 and 3/4 planets-per-system in the simulated observations match the observed relative frequencies.

The predicted ratio of observed single to double planet systems, however, is significantly smaller than observed. This effect is commonly referred to as the Kepler dichotomy as it may hint at an additional population of planet(ary system)s (Johansen et al., 2012; Hansen & Murray, 2013; Ballard & Johnson, 2016; Moriarty & Ballard, 2016). However, the dichotomy may also be an artifact of the planet detection pipeline (Zink et al., 2019) or a complex feature in the architectures of planetary systems (Zhu et al., 2018; Sandford et al., 2019) – though we do not see such a feature arise in the Bern models. We can account for this dichotomy following the procedure in Mulders et al. (2018), i.e. by assigning isotropic mutual inclinations to a fraction of simulated planetary systems. The ratio of 1/2 planets-per-system, would lead to a better match to the data if $f_{\text{iso}}\approx 40\%$ of the systems are assumed to have high mutual inclinations. This value is consistent with the values derived in Mulders et al. (2018); 2019MNRAS.490.4575H, but significantly larger than the $5-19\%$ effect reported by Zink et al. (2019).

The frequencies of systems with 5 or more planets would still be over-predicted in the simulated observations after correcting for the dichotomy. Most likely, the large number of simulated planets close to the star (see also §5.3) increase the number of observable systems with large numbers of planets per system. We will address this in a future iteration of the Bern model.

5.2 Period Ratios

The period ratio between adjacent planets is another tracer of the amount of dynamical excitation and orbital damping during planet assembly. More gravitational interactions between protoplanets lead to wider spacings and fewer orbital resonances, while more damping by the gas leads to smaller spacings and more orbital resonances.

The period ratios of the multi-seed models follow a broad distribution with pronounced peaks at the location of orbital resonances (Fig. 7). The period ratio distribution of the 100 seed-core model is most similar in shape to the Kepler distribution, with a median at $\mathcal{P}\sim 2$ and a long tail towards large period ratios (Fig. 7a). The 20 and 50-core models have a period ratio distribution that peaks at shorter period ratios than Kepler, both with a median at $\mathcal{P}\sim 1.5$. Orbital resonances, indicated with the dashed gray lines, appear more prominent when fewer planetary seed cores are considered, hinting at fewer gravitational interactions between forming planets in those cases.

The overall shape of the period ratio distribution in a simulated survey of the 100 seed-core model, shown as a cumulative distribution in (Fig. 7b), is similar to the observed one. While the 20 and 50 core models have on average smaller period ratios than observed, the 100-core model is a close match to the observed distribution, indicating the role that increased dynamical interactions play in making the synthetic simulated populations more closely match the distribution observed with Kepler. This observational comparison clearly favors the 100 seed-core model. However, the period ratio distribution shows more structure around orbital resonances than the Kepler data, implying that more dynamical interactions between planets are needed to push more planets out of resonances and smooth out these structures.

5.3 Innermost Planet Location

The location of the innermost planet in a planetary system traces the location of the gas disk inner disk edge — which in turn is set by magnetospheric truncation — and how efficiently migrating planets are trapped there. The efficiency at which planets are pushed out of the trap by other planets in the system determines the distribution interior to the planet trap (e.g. Carrera et al., 2019). The distribution of the innermost planet in each system for the Bern models is shown in Figure 8a. The inner planet distributions in the 50 and 100 seed core models show peaks just short-wards of $\sim 0.1$ au. The 100 core model shows a long tail towards smaller star-planet separations, while the 50 seed-core model also displays a tail toward longer periods, more similar to what is inferred from Kepler. The distribution of the innermost planet in the 20-core model does not show a clear peak but a wide distribution instead. The peaks in the location of the innermost planet in the 50 and 100 core model nearly line up with the inferred peak from the Kepler data (Paper I), which is at $0.1$ au or $P=12$ days. However, the simulated distributions are skewed in a way that it over-predicts the number of planets at periods less than $10$ days and under-predicts the number of innermost planets at longer orbital periods.

The overabundance of short-period planets in the 50 Bern models is most apparent from the simulated observations. Figure 8b shows the distribution of the innermost observed planet in each planetary system in the simulated survey compared to the Kepler data. The simulated distributions peak around an orbital period of less than $2$ days, whereas the observed distribution peaks at longer periods ($P\approx 7$ days). This comparison highlights that the number of detected planets at short orbital periods is very sensitive to the location where planets are trapped. Because the transit probability increases rapidly with decreasing orbital period, the number of detectable planets increases to a level where it completely dominates the observational comparison if the occurrence rate at short orbital periods is high.

Combined with the over-abundance of all planets at short orbital periods (Fig. 4a), this suggest that the migration trap at the disk inner edge needs to be placed at a larger distance from the star. The inner edge of the gas disk in the current models is placed at 0.03 au and the innermost planetesimals are located at twice that distance. A planet migrating alone can migrate down to the inner edge of the computational gas disk at 0.03 au. Planets even further in get pushed there by other planets when migrating within a resonant convoy of multiple planets. A disk inner edge truncated at 0.1 au would likely provide a better match between the Bern models and the Kepler exoplanet population. This would facilitate a more in-depth comparison of the orbital period distribution interior to the disk inner edge to constrain the physics of migration and trapping. A disk inner edge at 0.1 au is consistent with observations of the dust and gas in protoplanetary disks (e.g. Salyk et al., 2011; Pinte et al., 2008).

The distribution of orbital periods is generally broad. Although the distribution does not have the exact same shape as inferred from Kepler, it is not clear that a range of disk inner edge locations, as suggested by Lee & Chiang (2017) is needed to explain the observations. Instead, the stochastic nature of dynamical interactions during planet formation may be responsible for creating the broad period distribution interior to the disk inner edge (see also Carrera et al. 2019)

6 Predictions

Population synthesis models make predictions for the occurrence and composition of planets even in regions that are not directly constrained by observations. Those predictions can now be informed and refined by empirical constraints from the observable planet population. Here, we briefly explore what the model that most closely matches the Kepler data predicts for the number and type of planets we expect to find in the habitable zone of a solar-mass star. While the Kepler stars span a range of stellar masses, the average stellar properties are representative of a sun-like star, and we leave an exploration of the spectral type dependence of habitable zone planets for a future paper.

Figure 9 shows the simulated masses and radii of planets in the habitable zone of the 100 seed-core model. The habitable zone is here defined as an orbital period range of $338-788$ days, corresponding to the conservative habitable zone (Kopparapu et al., 2013) around a sun-like star. The populations that match the overall Kepler statistics, i.e. $\eta=0.2$ stars have a planetary system from the simulated set, predicts an occurrence rate of $0.8-1.7\,R_{\oplus}$ planets of $\Gamma_{\oplus}=\frac{d^{2}N}{d\ln R\,d\ln P}=18\%$. This rate is lower than the extrapolated power-law distributions from Paper I at $\Gamma_{\oplus}=53\%$, which we have addressed in Pascucci et al. (2019).

More importantly, the model predicts which planets in this size range are rocky and which planets have substantial hydrogen atmospheres. The fraction of true rocky planets is only $\Gamma_{\oplus}=7\%$, though we caution that this number is likely to change as the models become better calibrated to the observed exoplanet populations. However, future direct imaging missions should contend with the possibility that a large fraction of observable exo-Earths could have substantial hydrogen/helium atmospheres.

7 Results and Discussion

We have shown the need to evaluate planet formation models in a framework that accounts for the current observational biases in exoplanets surveys. This can be achieved with epos both through direct comparisons with parametric distributions of planet properties and occurrence rates, as well as through performing simulated observations of synthetic planet populations.

The Bern planet population synthesis models makes quantitative predictions for the distribution of planet and planetary system properties that we have evaluated with epos using the Kepler exoplanet statistics. We generally find good agreement between the models and data over a large range of planetary (system) properties. The Bern population synthesis models with multiple (20, 50, 100) interacting seed planet cores per protoplanetary disk form co-planar planetary systems of sub-Neptunes at short orbital periods that are typical of the Kepler “super-earths”. In contrast, synthetic populations with one seed planet core per disk produces many sharp features in the planet period-radius diagram that are not seen in exoplanet surveys. The smooth observed distributions indicate that planets rarely form via the growth of isolated oligarchs, but instead grow through the interactions of multiple protoplanetary bodies.

A quantitative comparison also reveals a number of key differences that indicates parameters and processes that could be refined in the model to provide a better match to the Kepler database. These are: increasing the formation efficiency of mini-Neptunes, increasing the amount of gravitational interactions between proto-planets, moving the planet trap at the inner disk edge outward, and decreasing the overall occurrence of planets.

7.1 Formation efficiency of mini-Neptunes

Simulated surveys of the Bern models where 1 in 5 of stars have planetary systems drawn from the synthetic population provides roughly the right proportion of detectable Hot Jupiters ($\sim 1\%$), Warm Giants ($\sim 10\%$), and Hot Earths ($\sim 50\%$). The Bern models under-predict the occurrence rate of mini-Neptunes ($\sim 10\%$ simulated vs $\sim 60\%$ observed) compared to rocky planets and warm giants. This indicates that a larger fraction of planetary cores need to accrete gaseous envelopes than in the current models, but without triggering more runaway growth into giant planets. A better understanding of gas accretion is needed and a more efficient formation channel for planets with substantial gaseous envelopes has to be incorporated in the synthesis model to match the Kepler planets. A similar trend is seen in micro-lensing surveys that probe larger planet-star separations, where the Bern model under-predicts the occurrence of intermediate-mass planets (Suzuki et al., 2018).

7.2 Dynamical excitation of planetary orbits

While the period ratio distribution of the 100 seed core model match closely with the observed distribution, too many planets are in orbital resonances. This may be indicative of a lack of dynamical interactions between planets in the simulations, which may in turn be a result of the short (20 Myr) N-body integration time considered here. Systems could go unstable at a later phase, as was seen in Izidoro et al. (2017), and running the simulations longer could result in a better match to the observed systems. Alternatively, external perturbations from stellar or planetary companions (e.g. Malmberg & Davies, 2009; Hansen, 2017) or a reduction in the amount of gas to reduce orbital damping (Dawson et al., 2016) may be necessary components for a planet population synthesis model.

7.3 Trapping planets at the inner disk edge

The Bern model over-predicts the occurrence of planets at short orbital periods, and the innermost planets of multi-planet systems are located too close to their host stars by a factor four. Future simulations could utilize a larger inner disk gas radius of $0.1$ au where planets are trapped, which would still yield planets inside of that location due to dynamical scattering. Such a radius would be consistent with observations of protoplanetary disks (Salyk et al., 2011) and coincide with the break in exoplanet occurrence rates that was previously inferred to be connected with the gas disk inner edge (Mulders et al., 2015b; Lee & Chiang, 2017).

7.4 Planet occurrence rates

The Bern models predict a larger fraction of stars with planets than observed. The Kepler planet occurrence rates can be matched when only 1 in 5 of stars have planetary systems consistent with the simulated systems. A similar conclusion was reached by Fernandes et al. (2019) based on a comparison with more distant giant planets from radial-velocity surveys. A lower planet formation efficiency would provide a better match to the observed populations, for example if 4 in 5 stars form planets below the detection limits. Such a low planet formation efficiency could be achieved by extending the distribution of protoplanetary disks to include lower mass disk that form planets below the detection limit, or by lowering the formation efficiency of planetesimals such that fewer solids end up in planets. Alternatively, external factors such as binarity (Kraus et al., 2016) or cluster interactions (Cai et al., 2018) could reduce the fraction of stars with planets.

8 Future Outlook

The observational comparisons in this paper provide a means of properly evaluating how well synthetic planet populations match the planetary systems that have been detected to date. By understanding which diagnostics are matched between model and real systems, we are able to constrain the physical processes that occur during planet formation and calibrate the nuisance parameters used in such models. Iterating between these approaches allows us to make, with greater confidence, more quantitative predictions regarding planetary properties, such as bulk or atmospheric compositions, that will be measured by future observational facilities such as JWST and the next generation of telescopes.

We acknowledge a constructive and positive report from an anonymous referee. We would like to thank Bertram Bitsch, Daniel Angerhausen, Rachel Fernandes, Leslie Rogers, and Aaron Hamann for providing feedback on the manuscript and code. 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 NASAs Nexus for Exoplanet System Science (NExSS) research coordination network sponsored by NASA’s Science Mission Directorate. C.M. and A.E. acknowledge the support from the Swiss National Science Foundation under grant BSSGI0$\_$155816 “PlanetsInTime”. Parts of this work have been carried out within the framework of the NCCR PlanetS supported by the Swiss National Science Foundation.

Appendix A Summary Statistics

Figure 10 show the summary statistics of the simulated populations of the 100 seed core Bern model compared with the equivalent observables from Kepler. These statistics are calculated for planets with sizes between $0.5$ and $20\,R_{\oplus}$ and orbital periods between $1$ and $400$ days.

The multi-planet frequency — the number of detectable planets per star — of the 100 seed-core model is shown in Figure 10b. This statistic mainly traces the mutual inclination distribution but is also sensitive to planet spacings, sizes, and orbital periods. Most notably this statistic does not reflect the intrinsic multi-planet frequency (gray dotted line) as typically only a few planets in each multi-planet system are transiting, even in systems with low mutual inclinations.

Figure 10b shows the distribution of period ratios between adjacent observed planets. The distribution before applying detection biases (gray line) peaks at shorter orbital periods than the observable one (purple). Generally, debiasing shifts the intrinsic period ratio distribution to larger values because planet pairs with a planet at large orbital periods are less likely to be observed (e.g. Brakensiek & Ragozzine, 2016). The reason why we observe the opposite trend is the covariance between orbital period ratio and orbital period that can be seen in Figure 3b: simulated planets at larger distances from the star have smaller period ratios. Because those planets also have low detection probabilities, they drop out of the observable sample in higher proportions than short-period, large period ratio planets, biasing the distribution.