Formation of short-period planets by disk migration

TL;DR

This study proposes that short-period planets form through Type-I migration and resonance locking in protoplanetary disks, successfully reproducing observed period distributions and highlighting the roles of migration and dynamical interactions.

Contribution

The paper introduces a new model combining N-body simulations and Kepler observation modeling to explain the formation of short-period planets via resonance chains and disk migration.

Findings

Resonance chains migrate inward together in simulations.

Simulated period distribution matches Kepler observations between 1-2 R⊕.

Fewer closely packed planets are produced than observed.

Abstract

Protoplanetary disks are thought to be truncated at orbital periods of around 10 days. Therefore, origin of rocky short period planets with $P < 10$ days 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 disk. The embryos experienced gravitational scatterings, collisions, disk 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 disk, the inner planets are pushed past the edge due to the disk 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_\oplus$. 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 the the formation and evolution of close-in small exoplanets.

Tables

Table 1. Table 1: We ran planet formation simulations for three values of the disk metallicity Z 𝑍 Z . All models had a stellar mass of 1 ​ M ⊙ 1 subscript 𝑀 direct-product 1M_{\odot} and 125 embryos with a mass of 0.4 ​ M ⊕ 0.4 subscript 𝑀 direct-sum 0.4M_{\oplus} each. The embryos start at 1 au and the separations are chosen to follow the solid surface density Σ solid = Z ​ Σ gas subscript Σ solid 𝑍 subscript Σ gas \Sigma_{\rm solid}=Z\,\Sigma_{\rm gas} . The disk metallicity affects the initial location of embryos a 0 subscript 𝑎 0 a_{0} , as well as the disk opacity, temperature, and the strength and direction of disk torques.

Model	$Z$	$N_{\mathrm{emb}}$	$M_{\mathrm{emb}}$	$a_{\mathrm{in}}$	$a_{\mathrm{out}}$	$N_{\mathrm{runs}}$
Z05	0.005	125	0.4$M_{\oplus }$	1 AU	7.6 AU	200
Z10	0.010	125	0.4$M_{\oplus }$	1 AU	6.0 AU	200
Z20	0.020	125	0.4$M_{\oplus }$	1 AU	4.9 AU	200

Equations

$$r_{\mathrm{in}}\approx r_{\mathrm{co}}=\sqrt[3]{\frac{G{M}_{\star }}{{\Omega }_{\star }^2}}.$$

$${\Gamma }_{\mathrm{tot}}={\Gamma }_{\mathrm{Lind}}+{\Gamma }_{\mathrm{CoRot}}.$$

$${\log }_{10}\left(\frac{{\dot{M}}_{\mathrm{disk}}}{M_{\odot }/\mathrm{yr}}\right)=-8-1.4{\log }_{10}\left(\frac{t_{\mathrm{disk}}+10^5\mathrm{yr}}{10^6\mathrm{yr}}\right).$$

$${\dot{M}}_{\mathrm{disk}}$$

$$T$$

$${\Sigma }_{\mathrm{solid}}=Z{\Sigma }_{\mathrm{gas}},$$

$$P(A)=P_{\mathrm{tran}}(A)\cdot P_{\mathrm{pipe}}(A)$$

$$P(A,B)=P_{\mathrm{tran}}(A,B)\cdot P_{\mathrm{pipe}}(A)\cdot P_{\mathrm{pipe}}(B)$$

$$\dot{a}=\frac{-9}{2}{a}^{-11/2}{G}^{1/2}\frac{m{R}_{\star }^5}{M_{\star }^{1/2}{Q}_{\star }^{\prime }},$$

$${\Gamma }_{\mathrm{tot}}={\Gamma }_{\mathrm{L}}{\Delta }_{\mathrm{L}}+{\Gamma }_{\mathrm{C}}{\Delta }_{\mathrm{C}}$$

$${\Delta }_{\mathrm{L}}$$

$${\Delta }_{\mathrm{C}}$$

$$P_{\mathrm{e}}=\frac{1+{\left(\frac{e}{2.25h}\right)}^{1.2}+{\left(\frac{e}{2.84h}\right)}^6}{1-{\left(\frac{e}{2.02h}\right)}^4}.$$

$$\frac{{\Gamma }_{\mathrm{L}}}{{\Gamma }_0/{\gamma }_{\mathrm{eff}}}$$

$$\frac{{\Gamma }_{\mathrm{C}}}{{\Gamma }_0/{\gamma }_{\mathrm{eff}}}$$

$$x=-\frac{\partial \ln {\Sigma }_{\mathrm{gas}}}{\partial \ln r},\beta =-\frac{\partial \ln T}{\partial \ln r}.{\Gamma }_0={\left(\frac{q}{h}\right)}^2{\Sigma }_{\mathrm{gas}}{r}^4{\Omega }_k^2,$$

$${\gamma }_{\mathrm{eff}}$$

$$Q$$

$$\chi$$

$$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_{\mathrm{s}}=\frac{1.1}{{\gamma }_{\mathrm{eff}}^{1/4}}\sqrt{\frac{q}{h}}.$$

$$F(p)$$

$$G(p)$$

$$K(p)$$

$${\mathbf{a}}_{\mathrm{mig}}=-\frac{\mathbf{v}}{t_{\mathrm{mig}}},$$

$$t_e$$

$$t_i$$

$$t_{\mathrm{wave}}$$

$${\mathbf{a}}_e=-2\frac{(\mathbf{v}\cdot \mathbf{r})\mathbf{r}}{r^2{t}_e},{\mathbf{a}}_i=-\frac{v_z}{t_i}\mathbf{k}$$

Full text

††pubyear: 2019††pagerange: Formation of short-period planets by disk migration–A

1 Introduction

An important interaction between a protoplanetary disk and a magnetised star is the disruption of the disk at the Alfvén radius (e.g. Long et al., 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 disk 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 disk is expected to be approximately at the point of co-rotation with the host star,

[TABLE]

Figure 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 days. Therefore, it is likely that protoplanetary disks are usually truncated at around 10 days. This estimate is also consistent with Keck interferometric observations of FU Orionis, which suggest a disk inner edge at $0.07\pm 0.02$ AU, or at around 10 days (Eisner & Hillenbrand, 2011). Indeed, the occurrence rate of sub-Neptunes appears to drop sharply for $P<10$ days (e.g. Youdin, 2011; Mulders et al., 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.

For a single sub-Neptune planet experiencing Type-I migration, disk torques are expected to simply deposit the planet at the inner edge of the disk. This can occur because of the loss of Lindblad torques just inside the disk (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 disk (e.g. Tanaka et al., 2002; Masset et al., 2006b; Terquem & Papaloizou, 2007; Tsang, 2011). In a recent work, Lee & 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 super-exponential 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 et al., 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 disk is present, the resonant chains tend to be stabilised by eccentricity and inclination damping from the disk. After the disk 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 disk, 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 disk may be pushed into the cavity through resonant interaction with planets farther out (Brasser et al., 2018). Figure 1 illustrates an extension of that idea to an extended resonant chain. A similar result was also seen in Figure 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 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 that tackle this problem attempt to “reverse” the process to obtain planet occurrence rates (Howard et al., 2012; Petigura et al., 2013b, a; Fressin et al., 2013).

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 suggest very different size and period dependence. A review of some of these differences was presented in the SAG 13 report111https://exoplanets.nasa.gov/exep/exopag/sag/#sag13. 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 organised 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 disk torques. The formulas for the disk 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 disk torques acting on a planet: Lindblad torque, and co-rotation torque,

[TABLE]

For typical disk parameters, Lindblad torque is negative and co-rotation torque is positive. Lindblad torque normally exceeds the co-rotation torque so that migration is usually inward. But outward migration is also possible. For a steady-state accretion disk, where the negative density and temperature gradients add up to 3/2, outward migration requires $d\ln T/d\ln r<-0.87$ (Brasser et al., 2017).

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

[TABLE]

In this formula, $t_{\rm disk}=0$ corresponds to a stellar age of $t_{\star}=10^{5}$ yr, which roughly corresponds to the moment when the disk becomes gravitationally stable (Bitsch et al., 2015). Our simulations all begin at $t_{\rm disk}=1$ Myr, which roughly corresponds to the timescale of embryo formation at 1 AU (Kokubo & Ida, 2000).

Given the disk temperature and $\dot{M}_{\rm disk}$, we use Shakura & Sunyaev (1973) $\alpha$-viscosity and hydrostatic equilibrium to solve for the gas surface density $\Sigma$

[TABLE]

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

Finally, we assume that disk photoevaporation becomes an important process at $t_{\rm disk}=5$ Myr and it dissipates the disk quickly (e.g. Gorti et al., 2009). To simulate that, we fix the disk temperature and reduce its mass exponentially with an e-folding timescale of $10^{4}$ yr for 10 e-folds. At 5.1 Myr, we remove the disk 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 disk torques. In other words, the surface density of embryos reflects the solid surface density of the disk at $t_{\rm disk}\approx 0$. Therefore we set

[TABLE]

where $\Sigma_{\rm gas}$ is computed at $t_{\rm disk}=0$ (Equation 3). The full list of simulations is shown in Table 1. We vary the disk metallicity across $Z$ = 0.5%, 1%, and 2%. For each metallicity we run 200 simulations. In each simulation the innermost embryo is placed at 1 AU. We divide the disk into 125 radial bins starting at 1 AU, such that each bin has a solid mass of 0.4 $M_{\oplus}$, and we place a 0.4 $M_{\oplus}$ 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_{\rm Hill}$. Classical models of oligarchic growth predict embryo separations of 5-10 $R_{\rm Hill}$ (e.g. Kokubo & Ida, 2000). Therefore, these simulations effectively model the final stages of embryo formation.

The total solid mass available to make planets is 50 $M_{\oplus}$. All embryos are given an initial eccentricity of $e=0.002$ and inclination $I=0.10^{\circ}$, as well as random mean anomaly, argument of preicenter, and longitude of ascending mode, chosen uniformly between $0^{\circ}$ and $360^{\circ}$.

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 days, 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 et al. (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% and 50%. This is a slightly wider range than the one observed in the solar system.

Finally, for planets beyond 10 days, 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$ day 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 ExoplanetsSysSim222https://github.com/eford/ExoplanetsSysSim.jl 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 onto 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 catalog. 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. Figure 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.

[TABLE]

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,

[TABLE]

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 Figure 3,

For a given disk model, run 200 simulations to generate 200 simulated planetary systems. 2. 2.

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

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

Use the simulated exoplanet catalog to compute planet occurrence rates in different period bins (Section 3.3) 5. 5.

Repeat steps 2-4 two hundred times.

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 procedures also ensures that each system is selected many times, and is 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

Figure 4 shows six snapshots from one of our simulations. The simulation is for a disk with a metallicity of $Z=1$% (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 scale height of the disk (grey region). As the embryos grow into planets, dynamical separations increase and torques from the disk 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.

The planets locked into a resonant chain migrate together as a group as the disk evolves. Eventually, the innermost planet reaches a zero torque migration region (a planet trap Masset et al., 2006a) at the inner edge of the disk. If the planet were alone, that is where migration would stop, as was done in the model of Lee & 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 disk edge and inside the disk 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 disk is no longer sufficient to move the planet chain significantly.

While planets are embedded in the disk, the dampening of eccentricity and inclination helps stabilise the planet orbits. In the bottom two frames of Figure 4 we see what happens to the system after the disk 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).

Figure 5 shows the same simulation from the point of view of the evolving disk structure. The figure shows the net torque experienced by a planet, as a function of the planet’s mass and its location in the disk. 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 disk. As the disk 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).

As shown in Figure 5, inside a few hundred days, at any given point in time the sign of the torque depends mainly on mass, and not so much on period. But as the disk 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

Figure 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 days are likely to have extended gaseous envelopes (Carrera et al., 2018).

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

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 disk models in Table 1. Figure 7 shows the “observed” period distribution of our simulated systems after taking into account observational effects.

For short-period planets ($P<10$ days) we applied a size cut-off of $R_{\rm p}=1-2R_{\oplus}$ 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 days. 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$ days as rocky cores.

Case B

Treat all planets with $P>10$ days as Jupiter-radius planets.

This provides the minimum and maximum detection efficiencies for planets beyond 10 days. For planets with $P>10$ days, the size cut-off is $R_{\rm p}>1R_{\oplus}$ for both Kepler and simulated planets. In Figure 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 disks. A distribution of disk metallicities might reproduce the exoplanet period distribution better than any single disk model alone.

3.4 Period ratios

When it comes to the orbital period ratios of short period planets, our simulation results are somewhat ambiguous. Figure 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.

Looking at more distant neighbours, Figure 9 shows the period ratios between any short period planet and any planet with $10\;{\rm d}<P<100\;{\rm 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 disk metallicities) might resolve most of the discrepancy.

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 (Figure 8) — our simulations predict fewer planet pairs with period ratios $P_{\rm out}/P_{\rm in}\lesssim 2$ than is observed in the Kepler catalogue. We have identified three possible effects that might be responsible for this discrepancy:

• •

One possible explanation is that too many of our systems experience a dynamical instability after the disk dissipates. In dynamical instability, the gravitational interactions between planets cause the orbits to evolve chaotically, leading to orbit crossing and giant impacts. The last two frames of Figure 4 illustrate this concept. After these collisions, the system is left with a smaller number of more widely separated planets. Figure 10 shows the distribution of period ratios just before the disk dissipates ($t=5$ Myr) and long after the disk 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 there is a correlation between the frequency of dynamical instabilities and some other property of the planetary system, like the masses of the planets. Our simulations under-produce Earth-mass planets (Figure 6) and it stands to reason that low-mass systems should be stable at closer orbital separations. Indeed, as shown in Figure 11, Kepler systems with smaller planets 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.

Future work should investigate 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 disk properties, such as the total amount of mass in embryos, and the dynamical history of the resulting planetary system.

4.2 Ultra short period planets

Ultra short period planets are generally defined as planets with periods less than 1 or 2 days. As Figure 6 shows, our simulations generally struggle to produce planets with periods shorter than two days, 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 disk at 10 days. But Figure 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 day. It is possible that the population of USPs with $P<1$ day come from the minority of disks where the inner edge is at $P\sim 1$ day.

• •

We also did not take into account the effect of stellar tides. Stellar tides become efficient at very short orbital periods ($P\sim 1$ day) 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 planets semimajor axis to decay at a rate

[TABLE]

where $m$ is the planet mass, $M_{\star}$ and $R_{\star}$ are the stellar mass and radius, 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 et al., 2008; Schlaufman et al., 2010; Penev et al., 2012).

5 Conclusion

The magnetic coupling between protoplanetary disks and their host stars implies that most protoplanetary disks are likely to be truncated at the point of co-rotation with the star’s rotational period, or at about 10 days. We propose that short period planets ($P<10$ days) may form through the Type-I migration of a resonant chain of planets. The process of super-Earth 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 disk, 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 (Figure 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$ days were also broadly consistent with observation, but the exact distribution depended greatly on the disk 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 ultra short period planets. We speculate that the observed period ratios may be partially tied to the frequency of dynamical instabilities. As for the under-production of ultra short period planets, we suspect that stellar tides, and the fact that some disks probably have an inner edge closer to $P\sim 1$ day, 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 the the formation and evolution of short period planets.

Acknowledgements

D.C. acknowledges Angie Wolfgang who provided helpful discussions and insight on modelling transit detection biases. D.C.’s research was supported by an appointment to the NASA Postdoctoral Program within NASA’s Nexus for Exoplanet System Science (NExSS), administered by Universities Space Research Association under contract with NASA. E.B.F. 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 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.

Appendix A Disk torques

This appendix is a slightly condensed rewrite of Appendix A of Carrera et al. (2018) with the correction of two minor typos. Following Paardekooper et al. (2010, 2011), all the formulas below include a smoothing of the planet potential with a smoothing length of $b=0.4h$ where $h=H/r\approx 0.05$. In Type-I migration, a planet experiences a negative Lindblad torque $\Gamma_{\rm L}$ and a positive co-rotation torque $\Gamma_{\rm C}$. The total torque on the planet is given by

[TABLE]

The expressions for $\Delta_{\rm L}$ and $\Delta_{\rm C}$ were produced by Cresswell & Nelson (2008); Coleman & Nelson (2014); Fendyke & Nelson (2014) and collected by Izidoro et al. (2017),

[TABLE]

where $e$ and $i$ are the planet orbital eccentricity and inclination, where $e_{\rm f}=0.5h+0.01$, and

[TABLE]

The expressions for $\Gamma_{\rm L}$ and $\Gamma_{\rm C}$ were derived by Paardekooper et al. (2010, 2011)

[TABLE]

where $\xi=\beta-(\gamma-1)x$ is the negative of the entropy slope, $\gamma=1.4$ is the adiabatic index, $x$ is the negative of the surface density profile, $\beta$ is the temperature gradient, and $\Gamma_{0}$ is a scaling factor,

[TABLE]

and $q$ is the planet-star mass ratio. Next, we define the effective adiabatic index $\gamma_{\rm eff}$

[TABLE]

where $\rho$ is the gas volume density, $\kappa$ is the opacity and $\sigma$ is the Stefan-Boltzmann constant. Next, $p_{\nu}$ and $p_{\chi}$ 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_{\rm s}}$,

[TABLE]

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

[TABLE]

Equations 10-27 are enough to compute the total torque $\Gamma_{\rm tot}$. To implement disk torques as an external force in an N-body code, we follow Papaloizou & Larwood (2000) and Cresswell & Nelson (2008) in defining the planet migration timescale as $t_{\rm mig}=-L/\Gamma_{\rm tot}$, where $L$ is the planet’s orbital angular momentum. The force term becomes,

[TABLE]

where $\mathbf{v}$ is the planet’s instantaneous velocity. To implement eccentricity and inclination damping, we need their damping timescales ${\rm t_{e}}$ and ${\rm t_{i}}$. Expressions for these timescales have been derived by Papaloizou & Larwood (2000) and Tanaka & Ward (2004), and were later modified by Cresswell & Nelson (2006, 2008),

[TABLE]

where $m$ and $a$ are the planet’s mass and semimajor axis. Similar to Equation 28, we implement eccentricity damping and inclination damping, respectively, as the forces

[TABLE]

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.