Formation of short-period planets by disc migration

ABSTRACT

Protoplanetary discs are thought to be truncated at orbital periods of around 10 d. Therefore, the origin of rocky short-period planets with P < 10 d is a puzzle. We propose that many of these planets may form through the Type-I migration of planets locked into a chain of mutual mean motion resonances. We ran N-body simulations of planetary embryos embedded in a protoplanetary disc. The embryos experienced gravitational scatterings, collisions, disc torques, and dampening of orbital eccentricity and inclination. We then modelled Kepler observations of these planets using a forward model of both the transit probability and the detection efficiency of the Kepler pipeline. We found that planets become locked into long chains of mean motion resonances that migrate in unison. When the chain reaches the edge of the disc, the inner planets are pushed past the edge due to the disc torques acting on the planets farther out in the chain. Our simulated systems successfully reproduce the observed period distribution of short-period Kepler planets between 1 and 2 R_⊕. However, we obtain fewer closely packed short-period planets than in the Kepler sample. Our results provide valuable insight into the planet formation process, and suggests that resonance locks, migration, and dynamical instabilities play important roles in the formation and evolution of close-in small exoplanets.

planets and satellites: dynamical evolution and stability, planets and satellites: formation, planets and satellites: general

1 INTRODUCTION

An important interaction between a protoplanetary disc and a magnetized star is the disruption of the disc at the Alfvén radius (e.g. Long, Romanova & Lovelace 2005), where the magnetic energy density due to the star’s magnetic field equals the kinetic energy density due to the in-falling material. At the same time, the magnetic connection between the star and the disc causes the star’s rotational angular momentum to become ‘locked’ at a point near the Alfvén radius (e.g. Koenigl 1991; Ostriker & Shu 1995; Long et al. 2005). In other words, the inner edge of a protoplanetary disc is expected to be approximately at the point of co-rotation with the host star,

\begin{eqnarray*} r_{\rm in} \approx r_{\rm co} = \root 3 \of {\frac{G M_\star }{\Omega _\star ^2}}. \end{eqnarray*}

(1)

Fig. 1 shows the distribution of rotation periods of 433 young FGK stars in the Orion Nebula (Herbst et al. 2002), showing that most of these stars have rotation periods of around 10 d. Therefore, it is likely that protoplanetary discs are usually truncated at around 10 d. This estimate is also consistent with Keck interferometric observations of FU Orionis, which suggest a disc inner edge at 0.07 ± 0.02 au, or at around 10 d (Eisner & Hillenbrand 2011). Indeed, the occurrence rate of sub-Neptunes appears to drop sharply for P < 10 d (e.g. Youdin 2011; Mulders, Pascucci & Apai 2015; Petigura et al. 2018; Hsu et al. 2019). How these planets formed and how they obtained those short-orbital periods, is an interesting puzzle, and is the subject of this investigation.

Figure 1.

Rotational periods of 433 young FGK stars in the Orion Nebula (top), showing a sharp peak at ∼10 d (Herbst et al. 2002). The disc is thought to stop at the co-rotation point. Yet, the Kepler catalogue (middle) contains many superEarths interior to 10 d (Borucki et al. 2010). We propose that these short-period planets may form when a group of protoplanets become locked into a chain of mutual mean motion resonances (bottom) which then migrates past the inner edge of the disc.

For a single sub-Neptune planet experiencing Type-I migration, disc torques are expected to simply deposit the planet at the inner edge of the disc. This can occur because of the loss of Lindblad torques just inside the disc (Goldreich & Tremaine 1980), or a large co-rotation torque or the reflection of waves may halt migration at or just outside the inner edge of the disc (e.g. Tanaka, Takeuchi & Ward 2002; Masset, D’Angelo & Kley 2006b; Terquem & Papaloizou 2007; Tsang 2011). In a recent work, Lee and Chiang (2017) argued that this feature disqualifies Type-I migration as a mechanism for forming short-period planets as the period distributions that they obtained did not match observations.

In this work we show that this simple view of migration is not a complete story. Sub-Neptune planets typically form in multiple planet systems where dynamical effects between planets are important (e.g. Terquem & Papaloizou 2007). In this work we show that the combination of Type-I migration with resonance locks of multiple sub-Neptune planets not only produces short-period planets, but also predicts a period distribution that is consistent with that observed in the Kepler field.

Planet formation is thought to begin with the formation of small 1–100 km bodies called planetesimals (Chiang & Youdin 2010; Johansen et al. 2014). Planetesimals first experience superexponential runaway growth (Greenberg et al. 1978; Wetherill & Stewart 1989; Kokubo & Ida 1996) followed by a period of slower oligarchic growth, which ends when the oligarchs reach their isolation mass (Kokubo & Ida 1998, 2000; Thommes, Duncan & Levison 2003; Chambers 2006). The formation of oligarchs may overlap with, or be followed by, pebble accretion. Pebble accretion is likely to be particularly important for the formation of giant planet cores at large orbital separations (Lambrechts & Johansen 2012; Johansen & Lambrechts 2017). Once the planet reaches a certain size, the Lindblad and co-rotation torques cause the planet to migrate inward (Paardekooper et al. 2010). In the process of migration, the planets capture each other into chains of mean motion resonances with several planets each (Terquem & Papaloizou 2007). While the disc is present, the resonant chains tend to be stabilized by eccentricity and inclination damping from the disc. After the disc dissipates, these resonant chains often become dynamically unstable, which leads to planet mergers and a dynamically hotter system (Cossou et al. 2014; Izidoro et al. 2017).

Unlike the single planet case, a resonant chain typically does not stop at the edge of the disc as the outer planets continue to experience inward torques, which force the entire chain inward. Supersonic corrections to the classical Type-I migration torques predict that a planet at the inner edge of the disc may be pushed into the cavity through resonant interaction with planets farther out (Brasser et al. 2018). Fig. 1 illustrates an extension of that idea to an extended resonant chain. A similar result was also seen in fig. 3 of Carrera et al. (2018). In this paper we investigate whether the ‘migrating resonant chain’ mechanism could be responsible for the majority of the observed short-period planet population.

1.1 Transit selection effects

A key part of this investigation is an accurate model of detection biases in transit surveys in general, and the Kepler catalogue (Borucki et al. 2010) in particular. Planets with shorter periods have a higher transit probability, and the detection efficiency of transit surveys depends on the transit depth, number of transits available, and stellar noise. Most authors who tackle this problem attempt to ‘reverse’ the process to obtain planet occurrence rates (Howard et al. 2012; Fressin et al. 2013; Petigura, Howard & Marcy 2013a; Petigura, Marcy & Howard 2013b)

Here we take a different approach. Rather than relying on reverse-engineered occurrence rates, we forward model the selection effects in the Kepler catalogue. In other words, we take the planetary systems that come out our simulations and we apply selection effects and compare that directly against the Kepler catalogue. This approach has many advantages. First, there is a large variance in the occurrence rates reported in the literature and often suggests very different size and period dependence. A review of some of these differences was presented in the SAG 13 report.¹ These differences arise due to different assumptions regarding completeness, whether they account for false positives, and even variations in the statistical methodology. Choosing among the published occurrence rates can lead to different results. But while inferred occurrence rates vary, the final Kepler catalogue does not. So we go directly to the Kepler catalogue for our data comparison, eliminating a middle step (occurrence rates) that is uncertain and out of our control. An additional benefit of a forward modelling is the ability to use the latest Kepler catalogue (DR25), which is also the most robust, the most homogeneous, and best understood catalogue in terms of completeness.

Finally, some of our analysis relies on the joint detection probability of two planets in the same system in order to obtain a distribution of period ratios (Section 3.4). This type of information simply is not available in any occurrence rate because it depends on the mutual inclinations between planet orbits — two planets with a high mutual inclination will almost never be detected together regardless of their individual detection probabilities. This type of geometric effect can and is included in our forward model.

This paper is organized as follows. In Section 2 we describe our planet formation model, initial conditions, as well as our forward model of Kepler observations. In Section 3 we present our results. We discuss the implications of our results in Section 4, and conclude in Section 5.

2 METHODS

2.1 Formation model

We use some of the planet formation simulations published in Carrera et al. (2018). A full description of the planet formation model is included in this publication. We run N-body simulations using the mercury code with the hybrid integrator (Chambers 1999) and a user-defined force to compute disc torques. The formulas for the disc torques are included in Appendix  A, and were implemented into mercury by Izidoro et al. (2017). Here we simply note that, in Type-I migration, there are two sources of disc torques acting on a planet: Lindblad torque and co-rotation torque,

\begin{eqnarray*} \Gamma _{\rm tot} = \Gamma _{\rm Lind} + \Gamma _{\rm CoRot}. \end{eqnarray*}

(2)

For typical disc parameters, Lindblad torque is negative and co-rotation torque is positive. Lindblad torque normally exceeds the co-rotation torque so that planet migration is usually inward. However, the strength of the two torques depend differently on the planet mass and the temperature profile of the disc (see Appendix  A). As general rule, a steep temperature profile increases the strength of the co-rotation torque more than the Lindblad torque. For a steady-state accretion disc, where the density and temperature gradients add up to −3/2 (e.g. equation 4), outward migration requires dln T/dln r < −0.87 (Brasser, Bitsch & Matsumura 2017).

In all simulations the stellar mass is 1 M_⊙ and we use the protoplanetary disc model of Bitsch et al. (2015). The model provides formulas for the disc temperature as a function of metallicity (Z) and accretion rate |$\dot{M}_{\rm disc}$|⁠. As in Carrera et al. (2018), we tie the accretion rate to the age of the disc t_disc using

\begin{eqnarray*} \log _{\rm 10}\left(\frac{\dot{M}_{\rm disc}}{\mathrm{M}_\odot \, \mathrm{yr}^{-1}} \right) = - 8 - 1.4\, \log _{\rm 10}\left(\frac{t_{\rm disc}+ 10^5 \, {\rm yr}}{10^6 \, {\rm yr}} \right). \end{eqnarray*}

(3)

In this formula, t_disc = 0 corresponds to a stellar age of t_⋆ = 10⁵ yr, which roughly corresponds to the moment when the disc becomes gravitationally stable (Bitsch et al. 2015). Our simulations all begin at t_disc = 1 Myr, which roughly corresponds to the time-scale of embryo formation at 1 au (Kokubo & Ida 2000). We assume that by that time the disc has been truncated at 0.1 au (P ∼ 10 d), and that no further stellar spin evolution takes place.

Given the disc temperature and |$\dot{M}_{\rm disk}$|⁠, we use Shakura and Sunyaev (1973) α-viscosity and hydrostatic equilibrium to solve for the gas surface density Σ

\begin{eqnarray*} \dot{M}_{\rm disc} &=& 3 \pi \alpha H^2 \Omega \Sigma , \end{eqnarray*}

(4)

\begin{eqnarray*} T &=& \left(\frac{H}{r} \right)^2 \frac{G M_\odot }{r} \frac{\mu }{\mathcal {R}}, \end{eqnarray*}

(5)

where G is the gravitational constant, μ = 2.3 is the molecular mass, and |$\mathcal {R}$| is the gas constant.

Finally, we assume that disc photoevaporation becomes an important process at t_disc = 5 Myr and it dissipates the disc quickly (e.g. Gorti, Dullemond & Hollenbach 2009). To simulate that, we fix the disc temperature and reduce its mass exponentially with an e-folding time-scale of 10⁴ yr for 10 e-folds. At 5.1 Myr, we remove the disc entirely for computational convenience.

2.2 Initial conditions

We assume that planetesimals form early and do not experience significant migration either due to aerodynamic drag or to disc torques. In other words, the surface density of embryos reflects the solid surface density of the disc at t_disc ≈ 0. Therefore we set

\begin{eqnarray*} \Sigma _{\rm solid} = Z \,\, \Sigma _{\rm gas}, \end{eqnarray*}

(6)

where Σ_gas is computed at t_disc = 0 (equation 3). The full list of simulations is shown in Table 1. We vary the disc metallicity across Z  = 0.5 per cent, 1 per cent, and 2 per cent. For each metallicity we run 200 simulations. In each simulation the innermost embryo is placed at 1 au. We divide the disc into 125 radial bins starting at 1 au, such that each bin has a solid mass of 0.4 M_⊕, and we place a 0.4 M_⊕ embryo in the middle of each bin. This works out to a dynamical separation of 7.3 mutual Hill radii for the first two embryos in our baseline model (Z01). As the semimajor axes increase, the dynamical separations drop down to 0.4 R_Hill. Classical models of oligarchic growth predict embryo separations of 5–10 R_Hill (e.g. Kokubo & Ida 2000). Therefore, these simulations effectively model the final stages of embryo formation.

Table 1.

We ran planet formation simulations for three values of the disc metallicity Z. All models had a stellar mass of 1 M_⊙ and 125 embryos with a mass of 0.4 M_⊕ each. The embryos start at 1 au and the separations are chosen to follow the solid surface density |$\Sigma _{\rm solid} = Z\, \Sigma _{\rm gas}$|⁠. The disc metallicity affects the initial location of embryos a₀ as well as the disc opacity, temperature, and the strength and direction of disc torques.

Model	Z	Nemb	Memb	ain	aout	Nruns
Z05	0.005	125	0.4M⊕	1 au	7.6 au	200
Z10	0.010	125	0.4M⊕	1 au	6.0 au	200
Z20	0.020	125	0.4M⊕	1 au	4.9 au	200

The total solid mass available to make planets is 50 M_⊕. All embryos are given an initial eccentricity of e = 0.002 and inclination I = 0.10° as well as random mean anomaly, argument of preicentre, and longitude of ascending mode, chosen uniformly between 0° and 360°.

2.3 Planet radius

Before we can simulate Kepler observations, we need to assign each planet a physical radius, as that determines the transit depth. For short-period planets, with periods less than 10 d, we assume that any primordial atmosphere has been fully evaporated and the planet is rock–iron core (e.g. Carrera et al. 2018). To compute the radius of the rock–iron core we use the planet structure model of Zeng, Sasselov & Jacobsen (2016). Naturally, the planet radius depends partly on the iron mass fraction, which is unknown for exoplanets and varies widely in the Solar system. To account for this source of uncertainty, every time we compute a core radius we randomly pick an iron mass fraction between 20 per cent and 50 per cent. This is a slightly wider range than the one observed in the Solar system.

Finally, for planets beyond 10 d, the atmosphere may represent a very large fraction of the planet’s transit radius. Rather than modelling atmosphere processes which have many uncertainties, we simply consider two limiting cases: we either treat all P > 10 d planets as bare rocky cores or treat them all as Jupiter-radius planets. This gives us a range from the smallest and largest transit radii that are consistent with each simulation.

2.4 Simulating transit observations

We use the ExoplanetsSysSim² package to simulate transit observations of planetary systems by the Kepler spacecraft. The package computes the geometric transit probability, as well as the detection efficiency of the Kepler pipeline taking into account the stellar noise, transit depth, and number of transits. To obtain an accurate simulation of the number of transits, we use the observing window function for each Kepler star, including any gaps in the observations, such as when a star falls on to a bad pixel on the Kepler CCD.

After a simulated planetary system has been selected and planetary radii have been assigned (Section 2.3), we randomly select a random G-type star from the Kepler input catalogue. We put all the planets around that star, preserving all of their orbital properties (semimajor axis, eccentricity, inclination, longitude of ascending node, argument of periastron, and mean anomaly). Preserving the orbital angles allows us to accurately model the joint probability that any given pair of planets both transit.

To simulate observed occurrence rates (Section 3.3), we only need to consider each plant in isolation. Fig. 2 shows how we compute total probability that a given planet is detected. First, we compute the sky-averaged probability that a planet transits. This quantity is primarily a function of semimajor axis, eccentricity, and stellar radius. Then we compute the probability that, if the planet transits, it would be detected by the Kepler pipeline. This quantity is a function of the transit depth, number of transits, and the amount of stellar noise. The product of these two values gives us a total detection probability, which is proportional to the occurrence rate.

\begin{eqnarray*} P(A) = P_{\rm tran}(A) \cdot P_{\rm pipe}(A) \end{eqnarray*}

(7)

Figure 2.

Kepler detection model for a sample G-dwarf in the Kepler target list. The total probability that a planet is detected (bottom) is the product of the probability that the planet transits (top) and the detection efficiency of the Kepler pipeline (middle). These vary from star to star as they depend on the stellar noise and the observation window for each star. At longer orbital periods, the position of the planet within its orbit significantly affects the number of transits that occur within the observing window for that star, which affects the detection efficiency.

For computing the joint probability that both planets in a planet pair are detected (Section 3.4), we need to take into account the relative orientations and shapes of the two planet orbits. Using that information, we compute a sky-averaged probability that both planets in a pair transit. To compute the detection efficiency of the Kepler pipeline, we make the simplifying assumption that the pipeline detection efficiencies are independent, so we compute those values as in the single-planet case. Finally, the total probability of detection is the product of all three values,

\begin{eqnarray*} P(A,B) = P_{\rm tran}(A,B) \cdot P_{\rm pipe}(A) \cdot P_{\rm pipe}(B) \end{eqnarray*}

(8)

2.5 Bootstrap procedure

If we were to simply apply transit observations (Section 2.4) to the simulated systems, we would obtain estimated occurrence rates, but we would not have confidence intervals on those quantities. To obtain error bars, we use the bootstrap process sketched in Fig. 3:

• For a given disc model, run 200 simulations to generate 200 simulated planetary systems.

• Randomly select, with replacement, 200 systems from the list of simulations.

• Assign to each system a G-type Kepler target list, assign each planet an iron fraction between 20 per cent and 50 per cent, and simulate Kepler observations (Sections 2.3 and 2.4).

• Use the simulated exoplanet catalogue to compute planet occurrence rates in different period bins (Section 3.3)

• Repeat steps 2–4 200 times.

Figure 3.

Sketch of our bootstrap method. We run 200 planet formation simulations with the same disc parameters. We randomly select 200 systems (with replacement) and assign each system a random G-type star from the Kepler input catalogue, including the radius, stellar noise, and observation window for that star. We simulate the observations of these systems to produce a simulated exoplanet catalogue. This is 1 bootstrap sample. We draw 200 bootstrap samples to obtain a distribution of simulated exoplanet catalogues that we then compare with the Kepler exoplanet catalogue.

Each iteration of steps 2–4 produces a possible exoplanet catalogue and list of occurrence rates. Comparing across all 200 bootstraps samples, we can obtain medians and confidence intervals on all the occurrence rates. The confidence intervals thus produced naturally take into account any uncertainties due to small number statistics that may occur when only a few simulations produce planes in a certain period range. This bootstrap procedure also ensures that each system is selected many times and observed many times around different Kepler targets, with slightly different transit depths due to different iron fractions.

While the procedure described above is expressed in terms of single planet occurrence rates, the procedure for computing the frequency of planet pairs and period ratios is entirely analogous.

3 RESULTS

3.1 Example simulation

Fig. 4 shows six snapshots from one of our simulations. The simulation is for a disc with a metallicity of Z = 1 per cent (Model Z10 in Table 1). The embryos (blue dots) are initially arranged in near-flat near-circular orbits from 1 to 6 au. Early on, the embryos experience strong dynamical interactions and collisions, occasionally reaching inclinations as high as the scaleheight of the disc (grey region). As the embryos grow into planets, dynamical separations increase and torques from the disc dampen the planets’ orbital inclinations and eccentricities. This dampening, along with migration, causes the newly forming planets to settle into long chains of mutual mean motion resonances.

Figure 4.

Snapshots from a simulation from Model Z10. The protoplanetary disc begins to dissipate at 5 Myr. The grey region shows the scaleheight of the disc as atan(H/r). The run begins with 125 embryos, each with a mass of 0.4 M_⊕, at a disc age of 1 Myr. The embryos have strong dynamical interactions that excite their eccentricities and inclinations. As the embryos merge, disc torques dampen eccentricities and inclinations and cause planets to migrate inward. At t = 2 Myr the planets have settled into a chain of mean motion resonances, and a migration trap separates the planets into two groups. After 3 Myr the inner chain has pushed past the inner edge of the disc. The snapshot at 5 Myr is right before the disc dissipates. The last two panels show the effect of the dynamical instabilities that occur after the dampening effect of the gas disc has been removed at 5.1 Myr. This run produces two planets with P < 10 d.

The planets locked into a resonant chain migrate together as a group as the disc evolves. Eventually, the innermost planet reaches a zero torque migration region (a planet trap; Masset et al. 2006a) at the inner edge of the disc. If the planet were alone, that is where migration would stop, as was done in the model of Lee and Chiang (2017). However, the presence of nearby planets significantly changes the dynamics at this critical junction. The other planets in the chain continue to migrate inward, and in doing so, they force the inner planet past the disc edge and inside the disc inner cavity. This process continues, with more and more planets being pushed past the edge, until it no longer can. At some point the torque on the planets that remain in the disc is no longer sufficient to move the planet chain significantly.

While planets are embedded in the disc, the dampening of eccentricity and inclination helps stabilize the planet orbits. In the bottom two frames of Fig. 4 we see what happens to the system after the disc dissipates. With the dampening removed, the resonant chains eventually become dynamically unstable and break apart, leading to a new wave of collisions and a dynamically hotter system (see Izidoro et al. 2017; for a detailed discussion).

Fig. 5 shows the same simulation from the point of view of the evolving disc structure. The figure shows the net torque experienced by a planet, as a function of the planet’s mass and its location in the disc. The colour scale shows the strength and direction of the net torque on the planet, and the region of positive torque (i.e. outward migration) is marked by a black boundary. Changes in the torque map are caused by changes in the thermal structure of the disc. As the disc evolves, it cools, and the steepness of the temperature gradient (which determines the relative strength of the co-rotation torque) also changes (Bitsch et al. 2015).

Figure 5.

Evolution of the migration map for the protoplanetary disc of Model Z10. The disc dissipates after 5.1 Myr. The colourbar shows the torque on the planet, normalized by the effective adiabatic index γ_eff and |$\Gamma _0 = (q/h)^2\, \Sigma \, r^4 \Omega ^2$|⁠, where q is the planet–star mass ratio, and h is the disc aspect ratio. The black curve marks the zero-torque boundary. A positive torque leads to inward migration and negative torque leads to outward migration. However, when a group of planets are locked into mean motion resonances, the resonant chain moves as a group, and the direction and speed of migration results from the net torque across all the planets. The white circles represent the planets in the sample simulation shown in Fig. 4.

As shown in Fig. 5, in a few 100 d, at any given point in time the sign of the torque depends mainly on mass, and not so much on period. But as the disc evolves and cools, the mass range that experiences a net outward torque decreases. The figure also shows the somewhat complex nature of migration for a resonant chain. In general, at any one time some planets in the chain experience negative torques and others experience positive torques, and it is the balance of torques on all the planets that determines the direction and speed of migration of the chain as a whole.

3.2 Planet radii

Fig. 6 shows the final periods and core radii of the simulated planets (blue) after taking into account the probability of transit and detection efficiencies of the Kepler pipeline. The vertical lines result from the fact that each planet is sampled many times and each time it receives a random iron fraction, which affects the planet’s transit radius. It is worth highlighting that this plot only shows the size of the planet’s core, but planets beyond 10 d are likely to have extended gaseous envelopes (Carrera et al. 2018).

Figure 6.

Periods and sizes of Kepler planets (black) and simulated planets (blue) after applying detection biases. For simulated planets, only the size of the rocky core is shown. The vertical lines reflect the fact that each planet is sampled many times, and assigned a random iron fraction of 20–50 per cent. Models Z05, Z10, and Z20 have disc metallicities of Z = 0.5 per cent, 1 per cent, and 2 per cent, respectively. The grey region marks short-period planets, which are likely to be photoevaporated, while the rest may have gaseous envelopes (Carrera et al. 2018).

The figure also shows that the disc metallicity clearly has an effect on the final sizes of the planets and the frequency of short-period planets. The simulations with low-metallicity discs (Model Z05) start out with more sparsely separated embryos that experience fewer collisions and take longer to reach the inner edge of the disc.

3.3 Period distribution

One of the most important tests for a formation model for short-period planets is its ability to reproduce the observed orbital period distribution. In this investigation we model planet radii and observational effects with the goal of making this particular test as robust as possible for each of the disc models in Table 1. Fig. 7 shows the ‘observed’ period distribution of our simulated systems after taking into account observational effects.

Figure 7.

The resonant–chain model (colour) predicts a planet period distribution that is consistent with that observed in the Kepler catalogue (black). The solid point is the median of 200 boostrap resamplings and the error bars cover the middle 90 per cent interval. The predicted period distribution takes into account the transit probability and the detection efficiency of the Kepler pipeline (Section 2.4). Because planets beyond 10 d may have large envelopes which increase the planet’s radius and detection efficiency, we consider two limiting cases: We either treat all planets as rocky cores (blue) or treat all planets beyond 10 d as Jupiters (red). Both the Kepler and simulated planets received the same size and period cut-offs (see the main text).

For short-period planets (P < 10 d) we applied a size cut-off of R_p = 1–2 R_⊕ for both Kepler and simulated planets, since our model is geared for the production of small planets. The figure shows that the proposed resonant–chain predicts a short-period planet distribution that is very consistent the one observed in the Kepler field.

The figure also extends the period distribution to orbital periods longer than 10 d. These are certainly not short-period planets, and it is likely that many of them possess extended gaseous envelopes (Carrera et al. 2018). For these planets, modelling the planet atmosphere is fraught with peril. Uncertainty in the planet radius translates into uncertainty in the implied occurrence rates, as larger transit depths are easier to detect by the Kepler pipeline. Instead, our approach is to consider two limiting cases:

Case A Treat all planets with P > 10 d as rocky cores.

Case B Treat all planets with P > 10 d as Jupiter-radius planets.

This provides the minimum and maximum detection efficiencies for planets beyond 10 d. For planets with P > 10 d, the size cut-off is R_p > 1 R_⊕ for both Kepler and simulated planets. In Fig. 7, results assuming that planets are rocky cores are shown in blue, and results assuming Jupiter-radius planets are shown in red. For the low-metallicity model (Z05) the rocky-core and Jupiter-radius models are both consistent with the data. For the higher metallicity models (Z10 and Z15), larger radii (red) are favoured in the last two period bins.

Finally, the planets in the Kepler catalogue will certainly have come for a diverse group of protoplanetary discs. A distribution of disc metallicities might reproduce the exoplanet period distribution better than any single disc model alone.

3.4 Period ratios

When it comes to the orbital period ratios of short-period planets, our simulation results are somewhat ambiguous. Fig. 8 shows the distribution of period ratios between any two short-period planets in the same system. In most of our models, the ‘observed’ period ratio between two short-period planets is typically larger than observed in the Kepler field. The main exception is Model Z05, which also has larger uncertainties due to small number statistics. It is possible that the apparently large period ratios might be due to higher mutual inclinations between planets, as that might make neighbouring planets less likely to both transit.

Figure 8.

Cumulative distribution of period ratios between short period (P < 10 d) planets. Black: Kepler catalogue of 1–2 R_⊕ planets. Blue: Simulated planets after applying Kepler detection biases. All short-period planets are treated as rocky cores with no gaseous envelope. The light blue bands cover the middle 90 per cent of bootstrap resamplings, and the solid line is the median. The models have disc metallicities of Z = 0.5 per cent, 1 per cent, and 2 per cent, respectively.

Looking at more distant neighbours, Fig. 9 shows the period ratios between any short-period planet and any planet with 10 d < P < 100 d. While none of the models fit the Kepler data perfectly, it seems likely that some combination of these models (i.e. a distribution of disc metallicities) might resolve most of the discrepancy.

Figure 9.

Cumulative distribution of period ratios between one long-period (P > 10 d) and one short-period (P < 10 d) planet. Black: Kepler catalogue of 1–2 R_⊕ planets. Colour: Simulated planets after applying Kepler detection biases. Short-period planets are treated as rocky cores. As limiting cases, planets beyond 10 d are treated as either bare rocky cores (blue) or as Jupiter-radius planets (red). The coloured bands cover the middle 90 per cent of bootstrap resamplings. The models have disc metallicities of Z = 0.5 per cent, 1 per cent, and 2 per cent, respectively.

4 DISCUSSION

4.1 Planet separations

Perhaps the most notable discrepancy between our simulation results and Kepler observations is in the separations between short-period planets (Fig. 8) – our simulations predict fewer planet pairs with period ratios P_out/P_in ≲ 2 than is observed in the Kepler catalogue. We have identified three possible effects that might be responsible for this discrepancy:

• One possibility is that our simulations may incorrectly predict the fraction of planetary systems that have a dynamical instability after the disc dissipates. A planetary system is said to be stable if orbits never cross. In a dynamical instability, gravitational interactions between planets cause the orbits to evolve chaotically, leading to orbit crossing and giant impacts. The last two frames of Fig. 4 illustrate this concept. After these collisions, the system is left with a smaller number of more widely separated planets. Fig. 10 shows the distribution of period ratios just before the disc dissipates (t = 5 Myr) and long after the disc dissipates (t = 25 Myr). At 5 Myr, most planet pairs are at or near mean motion resonances and have small period ratios. Over the next 20 Myr, most simulations have a dynamical instability that breaks these resonant chains and causes planets to collide. If our runs overestimate the frequency of instabilities, they may underestimate the number of systems with small period ratios.

• A related possibility is that our runs underproduce low-mass systems. Our simulations produce a somewhat narrow range of planet masses (Fig. 6), which corresponds to the limited range of initial conditions that we explored. Since low-mass planetary systems can be dynamically stable for smaller period ratios (Chambers, Wetherill & Boss 1996), a set of runs that produces too few low-mass systems would be expected to also underproduce systems with small period ratios. As seen in Fig. 11, Kepler systems with smaller planets do seem to prefer smaller period ratios.

• Finally, it is possible that there is a correlation with planet mass that is unrelated to dynamical instabilities. Perhaps the most natural alternative is that more massive planets may be more likely to be trapped at the 2:1 mean motion resonance.

Figure 10.

Histogram of the period ratio distribution of adjacent planets at two points in time. The disc lifetime is 5.1 Myr. The grey lines show the locations of the 2:1, 3:2, and 4:3 mean motion resonances. At t = 5 Myr (orange) many planet pairs are in or near mean motion resonances. At t = 25 Myr (blue) most systems have had collisions, leading to larger separations. This plot does not take observational effects into account. Those effects cause a small tilt toward smaller period ratios.

Figure 11.

Distribution of period ratios for small Kepler planets with P < 200 d inside two radius bins. Period ratios less than two preferentially come from planetary systems with small radius planets. This may be due to a combination of small planets being more stable at small orbital periods, or more massive planets preferentially being trapped at the 2:1 mean motion resonance.

Future work should explore how the dynamical history of planetary systems determines the observed distribution of period ratios. In addition, it will be important to explore possible connections between disc properties, such as the total amount of mass in embryos, and the dynamical history of the resulting planetary system (e.g. Izidoro et al. 2019).

4.2 Ultrashort-period planets

Ultrashort-period planets are generally defined as planets with periods less than 1 or 2 d. As Fig. 6 shows, our simulations generally struggle to produce planets with periods shorter than two d, and none of our simulations produced planets with periods less than a day. There are a couple of factors that we did not take into account that may be responsible for USPs:

• All of our simulations had the inner edge of the disc at 10 d. But Fig. 1 shows that there is in fact a distribution of stellar rotation rates (and hence, co-rotation radii) all the way down to about 1 d. It is possible that the population of USPs with P < 1 d come from the minority of discs where the inner edge is at P ∼ 1 d.

• We also did not take into account the effect of stellar tides. Stellar tides become efficient at very short orbital periods (P ∼ 1 d) and some authors have suggested that they may be responsible for the population of USPs (e.g. Lee & Chiang 2017).

Future investigations should investigate how these factors couple with the migration of resonant chains. Stellar tides pose a difficult challenge because the strength of the tides is not well understood. The tide raised by a planet on its host star causes the planet's semimajor axis to decay at a rate

\begin{eqnarray*} \dot{a} = \frac{-9}{2} a^{-11/2}G^{1/2} \frac{mR_\star ^5}{M_\star ^{1/2}Q_\star ^{\prime }}, \end{eqnarray*}

(9)

where m is the planet mass, M_⋆ and R_⋆ are the stellar mass and radius, respectively, and |$Q_\star ^{\prime }$| is the effective tidal quality factor (e.g. Goldreich & Soter 1966). Estimates of the quality factor |$Q_\star ^{\prime }$| vary by several orders of magnitude and depend strongly on the assumed properties of the planet and forcing period (e.g. Matsumura, Takeda & Rasio 2008; Schlaufman, Lin & Ida 2010; Penev et al. 2012).

5 CONCLUSION

The magnetic coupling between protoplanetary discs and their host stars implies that most protoplanetary discs are likely to be truncated at the point of co-rotation with the star’s rotational period, or at about 10 d. We propose that short-period planets (P < 10 d) may form through the Type-I migration of a resonant chain of planets. The process of superEarth formation from a population of planetesimals and embryos naturally leads to the formation of long chains of planets that are locked into mutual mean motion resonances (e.g. Terquem & Papaloizou 2007; Cossou et al. 2014; Izidoro et al. 2017; Carrera et al. 2018). When a migrating chain reaches the inner edge of the disc, the innermost planet experiences zero or even outward torques, but the aggregate inward torque from the outer planets in the chain force all planets into smaller orbital periods.

After a careful treatment of selection biases, using a forward model of the transit probabilities, and the detection efficiency of the Kepler pipeline (Sections 1.1 and 2.4), we showed that the ‘migrating chain’ model predicts period distributions that are consistent with those observed in the Kepler field (Fig. 7). Therefore, we conclude that a majority of short-period planets may well be produced through the migration of a resonant chain of planets. The observed period ratios between short-period planets and planets with P > 10 d were also broadly consistent with observation, but the exact distribution depended greatly on the disc properties.

The main area where the model needs improvement is that the model predicts slightly larger period ratios between short-period planets than what is observed, and the models generally struggle to produce ultrashort-period planets. We speculate that the observed period ratios may be partially tied to the frequency of dynamical instabilities. As for the underproduction of ultrashort-period planets, we suspect that stellar tides, and the fact that some discs probably have an inner edge closer to P ∼ 1 d, probably play a role.

Our results provide valuable insight into the planet formation process and suggests that resonance locks, migration, and dynamical instabilities play important roles in the formation and evolution of short-period planets.

ACKNOWLEDGEMENTS

DC acknowledges Angie Wolfgang who provided helpful discussions and insight on modelling transit detection biases. DC’s research was supported by an appointment to the National Aeronautics and Space Administration (NASA) Post-doctoral Program within NASA’s Nexus for Exoplanet System Science (NExSS), administered by Universities Space Research Association under contract with NASA. EBF acknowledges support from NASA Exoplanet Research Program award NNX15AE21G. The results reported herein benefitted 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. The Center for Exoplanets and Habitable Worlds is supported by the Pennsylvania State University, the Eberly College of Science, and the Pennsylvania Space Grant Consortium. We gratefully acknowledge support from National Science Foundation (NSF) grant MRI-1626251. This research or portions of this research were conducted with Advanced CyberInfrastructure computational resources provided by The Institute for CyberScience at The Pennsylvania State University (http://ics.psu.edu), including the CyberLAMP cluster supported by NSF grant MRI-1626251. AI thanks financial support provided by Fundação de Amparo à Pesquisa do Estado de São Paulo (FAPESP) via grants #2016/12686-2 and #2016/19556-7. This paper includes data collected by the Kepler mission. Funding for the Kepler mission is provided by the NASA Science Mission directorate.

Footnotes

1

https://exoplanets.nasa.gov/exep/exopag/sag/#sag13

2

https://github.com/eford/ExoplanetsSysSim.jl

APPENDIX A: DISC TORQUES

This appendix is a slightly condensed rewrite of appendix  A of Carrera et al. (2018) with the correction of some typos. Following Paardekooper et al. (2010); Paardekooper, Baruteau & Kley (2011); all the formulas below include a smoothing of the planet potential with a smoothing length of b = 0.4 h, where h = H/r ≈ 0.05. In Type-I migration, a planet experiences a negative Lindblad torque Γ_L and a positive co-rotation torque Γ_C. The total torque on the planet is given by

\begin{eqnarray*} \Gamma _{\rm tot} = \Gamma _{\rm L} \Delta _{\rm L} + \Gamma _{\rm C} \Delta _{\rm C} . \end{eqnarray*}

(A1)

The expressions for Δ_L and Δ_C were produced by Cresswell and Nelson (2008); Coleman and Nelson (2014); Fendyke and Nelson (2014); and collected by Izidoro et al. (2017),

\begin{eqnarray*} \Delta _{\rm L} &=& \left[ P_{\rm e} + \frac{P_{\rm e}}{|P_{\rm e}|} \times \left\lbrace 0.07 \left(\frac{i}{h}\right) + 0.085\left(\frac{i}{h}\right)^4\right.\right.\nonumber \\ && \left.\left. - \,0.08\left(\frac{e}{h} \right) \left(\frac{i}{h} \right)^2 \right\rbrace \right]^{-1} \end{eqnarray*}

(A2)

\begin{eqnarray*} \Delta _{\rm C} &=& \exp \left(\frac{e}{e_{\rm f}} \right) \left\lbrace 1 - \tanh \left(\frac{i}{h} \right) \right\rbrace, \end{eqnarray*}

(A3)

where e and i are the planet orbital eccentricity and inclination, where e_f = 0.5h + 0.01 and

\begin{eqnarray*} P_{\rm e} = \frac{ 1 + \left(\frac{e}{2.25\,h}\right)^{1.2} + \left(\frac{e}{2.84\,h}\right)^6 }{ 1 - \left(\frac{e}{2.02\,h}\right)^4 }. \end{eqnarray*}

(A4)

The expressions for Γ_L and Γ_C were derived by Paardekooper et al. (2010, 2011):

\begin{eqnarray*} \frac{\Gamma _{\rm L}}{\Gamma _0/\gamma _{\rm eff}} &=& (-2.5 -1.7\beta + 0.1x), \end{eqnarray*}

(A5)

\begin{eqnarray*} \frac{\Gamma _{\rm C}}{\Gamma _0/\gamma _{\rm eff}} &=& 1.1\left(\frac{3}{2}-x\right) F(p_\nu)~G(p_\nu)\nonumber \\ &&+\, 0.7\left(\frac{3}{2}-x\right) (1 - K(p_\nu)) \nonumber \\ &&+\, \frac{7.9 \xi }{\gamma _{\rm eff}} F(p_\nu) F(p_\chi) \sqrt{G(p_{\rm \nu }) G(p_\chi)}\nonumber \\ &&+\,\left(2.2 - \frac{1.4}{\gamma_{\rm eff}}\right)\xi \sqrt{(1 - K(p_\nu))~(1 - K(p_\chi))} \,, \end{eqnarray*}

(A6)

where ξ = β − (γ − 1)x is the negative of the entropy slope, γ = 1.4 is the adiabatic index, x is the negative of the surface density profile, β is the temperature gradient, and Γ₀ is a scaling factor,

\begin{eqnarray*} x = - \frac{\partial {\rm ln} ~\Sigma _{\rm gas}}{\partial {\rm ln}~r}, ~~~ \beta = - \frac{\partial {\rm ln} ~T}{\partial {\rm ln}~r}. ~~~ \Gamma _0 = \left(\frac{q}{h}\right)^2 \Sigma _{\rm gas} r^4 \Omega _k^2, \end{eqnarray*}

(A7)

and q is the planet–star mass ratio. Next, we define the effective adiabatic index γ_eff

\begin{eqnarray*} \gamma _{\rm eff} = \frac{2~Q \gamma }{ \gamma Q + \frac{1}{2} \sqrt{2 \sqrt{(\gamma ^2Q^2+1)^2-16 Q^2(\gamma -1)}+2 \gamma ^2 Q^2 - 2} }, \nonumber\\ \end{eqnarray*}

(A8)

\begin{eqnarray*} Q &=& \frac{2 \chi }{3~h^3 r^2 \Omega _k}, \end{eqnarray*}

(A9)

\begin{eqnarray*} \chi &=& \frac{16~\gamma ~(\gamma -1)~\sigma T^4}{3 \kappa ~\rho ^2~(hr)^2~\Omega _k^2}, \end{eqnarray*}

(A10)

where ρ is the gas volume density, κ is the opacity, and σ is the Stefan–Boltzmann constant. Next, p_ν and p_χ are parameters that govern the viscous saturation and the thermal saturation, respectively. They are defined in terms of the non-dimensional half-width of the horseshoe region x_s,

\begin{eqnarray*} p_\nu = \frac{2}{3}\sqrt{\frac{r^2\Omega _k}{2\pi \nu }x_s^3}, ~~~ p_\chi = \sqrt{\frac{r^2\Omega _k}{2\pi \chi }x_s^3}, ~~~ x_{\rm s} = \frac{1.1}{\gamma _{\rm eff}^{1/4}}\sqrt{\frac{q}{h}}. \end{eqnarray*}

(A11)

Finally, we give the expression for functions F, G, and K,

\begin{eqnarray*} F(p) &=& \frac{1}{1+ \left(\frac{p}{1.3}\right)^2} \end{eqnarray*}

(A12)

\begin{eqnarray*} G(p) = \left\lbrace \begin{array}{l{@}{\quad}l}\frac{16}{25} \left(\frac{45\pi }{8}\right)^{3/4} p^{3/2} & {\rm if} p \lt \sqrt{\frac{8}{45\pi }}\\ 1 - \frac{9}{25} \left(\frac{8}{45\pi }\right)^{4/3} p^{-8/3} & {\rm otherwise} \end{array} \right. \end{eqnarray*}

(A13)

\begin{eqnarray*} K(p) &=& \left\lbrace \begin{array}{l{@}{\quad}l}\frac{16}{25} \left(\frac{45\pi }{28} \right)^{3/4} p^{3/2} & {\rm if} p \lt \sqrt{\frac{28}{45\pi }}\\ 1 - \frac{9}{25} \left(\frac{28}{45\pi } \right)^{4/3} p^{-8/3} & {\rm otherwise} \end{array} \right. \end{eqnarray*}

(A14)

Equations (A1)–(A14) are enough to compute the total torque Γ_tot. To implement disc torques as an external force in an N-body code, we follow Papaloizou and Larwood (2000); Cresswell and Nelson (2008) in defining the planet migration time-scale as t_mig = −L/Γ_tot, where L is the planet’s orbital angular momentum. The force term becomes,

\begin{eqnarray*} \mathbf {a}_{\rm mig} = -\frac{\mathbf {v}}{t_{\rm mig}}, \end{eqnarray*}

(A15)

where |$\mathbf {v}$| is the planet’s instantaneous velocity. To implement eccentricity and inclination damping, we need their damping time-scales t_e and t_i. Expressions for these time-scales have been derived by Papaloizou and Larwood (2000); Tanaka and Ward (2004), and were later modified by Cresswell and Nelson (2006, 2008),

\begin{eqnarray*} t_e = \frac{t_{\rm wave}}{0.780} \left[ 1-0.14\left(\frac{e}{h}\right)^2 + 0.06\left(\frac{e}{h}\right)^3+ 0.18\left(\frac{e}{h}\right)\left(\frac{i}{h}\right)^2 \right];\nonumber\\ \end{eqnarray*}

(A16)

\begin{eqnarray*} t_i &=& \frac{t_{\rm wave}}{0.544} \left[ 1-0.3\left(\frac{i}{h}\right)^2 + 0.24\left(\frac{i}{h}\right)^3+ 0.14\left(\frac{e}{h}\right)^2\left(\frac{i}{h}\right) \right]; \nonumber\\ \end{eqnarray*}

(A17)

\begin{eqnarray*} t_{\rm wave} &=& \left(\frac{\mathrm{M}_{\odot }}{m}\right) \left(\frac{\mathrm{M}_{\odot }}{\Sigma _{\rm gas} a^2}\right) h^4 \Omega _k^{-1}, \end{eqnarray*}

(A18)

where m and a are the planet’s mass and semimajor axis. Similar to equation (A15), we implement eccentricity damping and inclination damping, respectively, as the forces

\begin{eqnarray*} \mathbf {a}_e = -2~\frac{(\mathbf {v \cdot r})\mathbf {r}}{r^2 t_e}, ~~~ \mathbf {a}_i = -2~\frac{v_z}{t_i}\mathbf {k} \end{eqnarray*}

(A19)

where |$\mathbf {r}$| and |$\mathbf {v}$| are the position and velocity vectors of the planet, v_z is the z component of the planet’s velocity, and |$\mathbf {k}$| is the unit vector in the z direction. The formula for |$\mathbf {a}_i$| includes a factor of 2, which corrects a typo in Carrera et al. (2018); Cresswell and Nelson (2006, 2008).