Breaking the Chains: Hot Super-Earth systems from migration and disruption of compact resonant chains

TL;DR

This study uses N-body simulations to explore the formation and evolution of hot super-Earth systems, revealing that most resonant chains become unstable after gas dissipation, aligning with observed system architectures.

Contribution

The paper introduces a comprehensive simulation framework that models the formation, migration, and post-dissipation evolution of super-Earth systems, highlighting the role of dynamical instability.

Findings

Most resonant chains become unstable after gas dissipation.

Stable resonant systems constitute less than 25% of all systems.

The observed high rate of single-planet systems suggests a significant dynamical instability in super-Earth systems.

Abstract

"Hot super-Earths" (or "Mini-Neptunes") between 1 and 4 times Earth's size with period shorter than 100 days orbit 30-50\% of Sun-like type stars. Their orbital configuration -- measured as the period ratio distribution of adjacent planets in multi-planet systems -- is a strong constraint for formation models. Here we use N-body simulations with synthetic forces from an underlying evolving gaseous disk to model the formation and long-term dynamical evolution of super-Earth systems. While the gas disk is present, planetary embryos grow and migrate inward to form a resonant chain anchored at the inner edge of the disk. These resonant chains are far more compact than the observed super-Earth systems. Once the gas dissipates resonant chains may become dynamically unstable. They undergo a phase of giant impacts that spreads the systems out. Disk turbulence has no measurable effect on the outcome. Our simulations match observations if a small fraction of resonant chains remain stable, while most super-Earths undergo a late dynamical instability. Our statistical analysis restricts the contribution of stable systems to less than $25\%$. Our results also suggest that the large fraction of observed single planet systems does not necessarily imply any dichotomy in the architecture of planetary systems. Finally, we use the low abundance of resonances in Kepler data to argue that, in reality, the survival of resonant chains happens likely only in $\sim 5\%$ of the cases. This leads to a mystery: in our simulations only 50-60\% of resonant chains became unstable whereas at least 75\% (and probably 90-95\%) must be unstable to match observations.

Figures

Equations

$${\dot{M}}_{\text{gas}}=3\pi \alpha h^2{r}^2{\Omega }_k{\Sigma }_{gas},$$

$$T=h^2\frac{G{M}_{\odot }}{r}\frac{\mu }{\mathcal{R}}$$

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

$$\Re =\tanh \left(\frac{r-0.1}{0.005}\right).$$

$$x=-\frac{\partial ln{\Sigma }_{gas}}{\partial lnr},\beta =-\frac{\partial lnT}{\partial lnr}.$$

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

$$\begin{gathered} {\Delta }_{\text{L}}=[P_{\text{e}}+\frac{P_{\text{e}}}{|P_{\text{e}}|}\times \{0.07\left(\frac{i}{h}\right)+0.085{\left(\frac{i}{h}\right)}^4 \\ {-0.08\left(\frac{e}{h}\right){\left(\frac{i}{h}\right)}^2\}]}^{-1}, \end{gathered}$$

$$P_{\text{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}.$$

$${\Delta }_{\text{C}}=\exp \left(\frac{e}{e_{\text{f}}}\right)\left\{1-\tanh \left(\frac{i}{h}\right)\right\},$$

$$e_{\text{f}}=0.5h+0.01.$$

$${\Gamma }_{\text{C}}={\Gamma }_{\text{c,hs,baro}}F(p_{\nu })G(p_{\nu })+(1-K(p_{\nu })){\Gamma }_{\text{c,lin,baro}}+$$

$${\Gamma }_{\text{c,hs,ent}}F(p_{\nu })F(p_{\chi })\sqrt{G(p_{\nu })G(p_{\chi })}+$$

$$\sqrt{(1-K(p_{\nu }))(1-K(p_{\chi })}{\Gamma }_{\text{c,lin,ent}}.$$

$${\Gamma }_{\text{L}}=(-2.5-1.7\beta +0.1x)\frac{{\Gamma }_{\text{0}}}{{\gamma }_{\text{eff}}},$$

$${\Gamma }_{\text{c,hs,baro}}=1.1\left(\frac{3}{2}-x\right)\frac{{\Gamma }_{\text{0}}}{{\gamma }_{\text{eff}}},$$

$${\Gamma }_{\text{c,lin,baro}}=0.7\left(\frac{3}{2}-x\right)\frac{{\Gamma }_{\text{0}}}{{\gamma }_{\text{eff}}},$$

$${\Gamma }_{\text{c,hs,ent}}=7.9\xi \frac{{\Gamma }_{\text{0}}}{{\gamma }_{\text{eff}}^2},$$

$${\Gamma }_{\text{c,lin,ent}}=\left(2.2-\frac{1.4}{{\gamma }_{\text{eff}}}\right)\xi \frac{{\Gamma }_0}{{\gamma }_{\text{eff}}}$$

$$p_{\nu }=\frac{2}{3}\sqrt{\frac{r^2{\Omega }_k}{2\pi \nu }{x}_s^3},$$

$$x_{\mathrm{s}}=\frac{1.1}{{{\gamma }_{\text{eff}}}^{1/4}}\sqrt{\frac{q}{h}}.$$

$$p_{\chi }=\sqrt{\frac{r^2{\Omega }_k}{2\pi \chi }{x}_s^3},$$

$$\chi =\frac{16\gamma (\gamma -1)\sigma T^4}{3\kappa {\rho }^2(hr{)}^2{\Omega }_k^2},$$

$$Q=\frac{2\chi }{3h^3{r}^2{\Omega }_k}$$

$${\gamma }_{\text{eff}}=\frac{2Q\gamma }{\gamma Q+\frac{1}{2}\sqrt{2\sqrt{({\gamma }^2{Q}^2+1{)}^2-16Q^2(\gamma -1)}+2{\gamma }^2{Q}^2-2}}.$$

$$F(p)=\frac{1}{1+{\left(\frac{p}{1.3}\right)}^2,}$$

$$G(p)=\{\begin{array}{ll}\frac{16}{25}{\left(\frac{45\pi }{8}\right)}^{\frac{3}{4}}{p}^{\frac{3}{2}},& \text{if p<45π8}\\ 1-\frac{9}{25}{\left(\frac{8}{45\pi }\right)}^{\frac{4}{3}}{p}^{-\frac{8}{3}},& \text{otherwise}.\end{array},$$

$$K(p)=\{\begin{array}{ll}\frac{16}{25}{\left(\frac{45\pi }{28}\right)}^{\frac{3}{4}}{p}^{\frac{3}{2}},& \text{if p<45π28}\\ 1-\frac{9}{25}{\left(\frac{28}{45\pi }\right)}^{\frac{4}{3}}{p}^{-\frac{8}{3}},& \text{otherwise}.\end{array}.$$

$$t_m=-\frac{L}{{\Gamma }_{tot}}.$$

$$\begin{gathered} t_e=\frac{t_{wave}}{0.780}(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), \end{gathered}$$

$$\begin{gathered} t_i=\frac{t_{wave}}{0.544}(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)), \end{gathered}$$

$$t_{wave}=\left(\frac{M_{\odot }}{m_p}\right)\left(\frac{M_{\odot }}{{\Sigma }_{gas}{a}^2}\right)h^4{\Omega }_k^{-1},$$

Full text

††pubyear: 2016††pagerange: Breaking the Chains: Hot Super-Earth systems from migration and disruption of compact resonant chains–Breaking the Chains: Hot Super-Earth systems from migration and disruption of compact resonant chains

1 Introduction

Among the thousands of confirmed exoplanets, hot super-Earths or mini-Neptunes – with radii between 1 and ${\rm 4R_{\oplus}}$ (1 $<M_{\oplus}$< 20), orbiting very close to their host stars – form by far the largest population (Mayor et al., 2011; Borucki et al., 2010, 2011; Lissauer et al., 2011b; Schneider et al., 2011; Howard et al., 2012; Fabrycky et al., 2012; Petigura et al., 2013; Howard, 2013; Dong & Zhu, 2013; Batalha et al., 2013; Howard, 2013; Fressin et al., 2013; Mullally et al., 2015). Statistical studies suggest that about 1 out of 3 solar type stars (FGK spectral types) host a super-Earth with orbital period shorter than 100 days (Mayor et al., 2011; Howard et al., 2012; Fressin et al., 2013; Petigura et al., 2013). Yet, close-in super Earths are often found in compact multi-planet systems (e.g. Lissauer et al. (2011a)). Their eccentricities and mutual orbital inclinations are estimated to statistically concentrate around low and moderate values (e $\lesssim$ 0.1-0.2; i $\lesssim$ 5 deg; Mayor et al. (2011); Lissauer et al. (2011a); Fang & Margot (2012)).

A fundamental open question in planet formation is: Where and how did hot super-Earth systems form and dynamically evolve?

Current models on the formation of systems of close-in super-Earths can be divided into two main categories (for reviews see Raymond et al., 2008; Raymond & Morbidelli, 2014; Morbidelli & Raymond, 2016): (1) in-situ accretion (Raymond et al., 2008; Hansen & Murray, 2012, 2013; Chiang & Laughlin, 2013; Hansen, 2014; Ogihara et al., 2015a, b) or (2) assembly of planets at moderate or larger distances from the star followed by inward gas-driven migration (Terquem & Papaloizou, 2007; Ida & Lin, 2008, 2010; McNeil & Nelson, 2010; Hellary & Nelson, 2012; Cossou et al., 2014; Coleman & Nelson, 2014, 2016). In-situ accretion models virtually comes in two flavors: (a) standard in-situ accretion models which invokes high mass disks from the beginning to allow multi-Earth masses planets to form in the inner regions (Hansen & Murray, 2012, 2013); and (b) drift-then-assembly model which proposes that these planets form in the innermost regions of the disk in consequence of a local concentration of pebbles or small planetesimals drifting inward due to gas drag (Chatterjee & Tan, 2014; Boley & Ford, 2013; Boley et al., 2014; Chatterjee & Tan, 2015; Hu et al., 2016).

The standard in-situ accretion model suffers from too many fundamental issues to be plausible (see discussions in Raymond & Cossou (2014), Schlichting (2014), Schlaufman (2014), Izidoro et al. (2015), Ogihara et al. (2015a), and Chatterjee & Tan (2015)). For example, it ignores the effects of planet-disk gravitational interaction. Planets forming in-situ grow extremely fast because of short dynamical timescales and the required abundant amount of mass in the inner regions of the disk (Hansen & Murray, 2012, 2013). They tend to reach masses large enough (Hansen & Murray, 2013; Hansen, 2014) to sufficiently perturb the surrounding gas in timescales much shorter than the expected lifetime of protoplanetary disks around young stars (Ogihara et al., 2015a). Planet-disk gravitational interaction leads to angular momentum transfer between the planet and the disk (see recent review by Baruteau et al. (2014)) which typically causes orbital radial migration(Goldreich & Tremaine, 1979, 1980; Lin & Papaloizou, 1979, 1986; Ward, 1986; Artymowicz, 1993; Ward, 1997; Tanaka et al., 2002; Sari & Goldreich, 2004), as well as eccentricity and inclination damping of the planets’ orbit (Papaloizou & Larwood, 2000; Goldreich & Sari, 2003; Tanaka & Ward, 2004). Ignoring planet-disk gravitational interaction (migration and orbital tidal damping) in models of in-situ formation is not self-consistent. Moreover, if planets eventually migrate during the gas disk phase they would not be forming truly in-situ (Ogihara et al., 2015a).

The drift-then-assembly models (Chatterjee & Tan, 2014; Boley & Ford, 2013; Boley et al., 2014; Chatterjee & Tan, 2015; Hu et al., 2016) are quite promising but what actually happens near the inner edge of the disk, which is dangerously close to the sublimation line of silicates is not clear (Morbidelli et al., 2016). Instead, the formation of massive objects beyond the snowline within the lifetime of the disk seems to be generic from theoretical considerations (e.g. Morbidelli et al., 2015). Migration seems to be a generic process as well. So, we focus in this paper on the scenario where planets are assembled at moderate or larger distances from the star and then moved close to the star by gas-driven migration.

Earth-mass planets typically are not able to open a gap in the gaseous disk (e.g. Papaloizou & Lin (1984); Crida et al. (2006)) and migrate in the type-I regime (e.g. Ward (1997); Kley & Nelson (2012)). Sophisticated hydrodynamical simulations including thermodynamical effects show that type-I migration is very sensitive to the disk properties. Planets may either migrate inward or outward depending on the combination of different torques from the disk (Paardekooper & Mellema, 2006, 2008; Baruteau & Masset, 2008; Paardekooper & Papaloizou, 2008; Kley et al., 2009). Indeed, the Lindblad torque tend to push the planet inwards (e.g. Ward (1986, 1997)) but depending on the planet mass the gas-flowing in coorbital motion with the planet may exert a strong torque capable of stopping or even reversing the direction of type-I migration (Kley & Crida (2008); Paardekooper et al. (2010, 2011). There are locations within the disk where the net torque is zero (Lyra et al 2010; Horn et al 2012; Cossou et al 2013; Pierens et al 2013; Bitsch et al 2013, 2014, 2015). However, as the disk dissipates and cools, thermodynamics and viscous effects becomes less important and planets are released to migrate inward (Lyra et al., 2010; Horn et al., 2012; Bitsch et al., 2014). It is hard to imagine planets completely escaping inward type-I migration, although disk winds may be a possible explanation of a global suppression of type-I migration in the inner parts of the disk (Ogihara et al., 2015a; Suzuki et al., 2016).

One of the main criticisms of the inward migration model for the origins of close-in super-Earths comes from the fact that many planet pairs are near but not exactly in first order mean motion resonances (Lissauer et al. (2011b); Fabrycky et al. (2014)). There is a prominent excess of planet pairs just outside first order mean motion resonances (Fabrycky et al., 2014). Planet migration models predicts that as the disk dissipates planets should migrate inward and pile up in long chains of mean motion resonances (Terquem & Papaloizou (2007); Raymond et al. (2008); McNeil & Nelson (2010); Rein (2012); Rein et al. (2012); Horn et al. (2012); Ogihara & Kobayashi (2013); Haghighipour (2013); Cossou et al. (2014); Raymond & Cossou (2014); Ogihara et al. (2015a); Liu et al. (2015, 2016)), in stark contrast to the observations. There are three reasons not to reject the migration model. First, a fraction of planet pairs are indeed in first order mean motion resonance (Lissauer et al. (2011b); Fabrycky et al. (2014)). For example, the recently discovered Kepler 223 system presents a very peculiar orbital configuration. Planets in this system are locked in a chain of resonances which mostly likely could be explained by convergent migration during the gas disk phase (Mills et al., 2016). The TRAPPIST-1 system is another example of planetary system with multiple planets in a resonant chain (Gillon et al., 2017). Second, a number of mechanisms have been proposed for accounting for the excess of planet pairs outside mean motions resonances which could consistently operate with the inward-migration model. Among them are star-planet tidal dissipation (Papaloizou & Terquem, 2010; Papaloizou, 2011; Lithwick & Wu, 2012; Delisle et al., 2012; Batygin & Morbidelli, 2013; Delisle et al., 2014; Delisle & Laskar, 2014), planet scattering of leftover planetesimals (Chatterjee & Ford, 2015), turbulence in the gaseous disk (Rein, 2012; Rein et al., 2012) (but see Section 4), interaction with wake excited by other planets (Baruteau & Papaloizou, 2013) and the effects of asymmetries in the structure of the protoplanetary disk (Batygin, 2015). Third, and most importantly, planets that form in resonance may not remain in resonance. Rather, planets in resonance can become unstable when the gas disk dissipates (Terquem & Papaloizou, 2007; Ogihara & Ida, 2009; Cossou et al., 2014). The systems that survive instabilities are not in resonance. To summarize, the lack of resonance between planet pairs should not be taken as evidence against inward-migration (Goldreich & Schlichting, 2014).

In this paper we assume that hot super-Earth systems form by inward gas-driven migration. Our numerical simulations model the dynamical evolution of Earth-mass planets in evolving protoplanetary disks (Williams & Cieza, 2011). Our nominal simulations include the effects of type-I migration and orbital eccentricity and inclination damping due to the gravitational interaction with the gas. We have also tested in our model the effects of stochastic forcing from turbulent fluctuations in the disk of gas. After the gas disk’s dissipation, our simulations were continued up to 100 Myr. The goal of this study is to help elucidating the following question: is the inward migration model for the origins of hot super-Earths systems consistent with observations?

This paper is structured as follows. In Section 2 we describe our methods and disk model. In Section 3 we present the results of simulations of our fiducial model. In Section 4 we describe our turbulent model and the results of simulations including these effects. In Section 5 we discuss about the role of dynamical instabilities post-gas disk dissipation and the final dynamical architecture of planetary systems produced in our simulations. In section 6 we compare the results of our fiducial and turbulent models with observations. In Section 7 we compare our results with other models in the literature. In Section 8 we discuss about our main results. Finally, in Section 9 we summarize our conclusions.

2 Methods

We use N-body numerical simulations to study the dynamical evolution of multiple Earth-mass planets in evolving protoplanetary disks. We also follow the subsequent phase of dynamical evolution of formed planetary systems post gas-disk dissipation. During the gaseous phase, we mimic the effects of the disk of gas on the planets by applying artificial forces onto the planets (or protoplanetary embryos). These forces were calibrated from truly hydrodynamical simulations. In this section we describe our gas disk model, followed by the details of our prescription for type-I migration, eccentricity and orbital inclination damping, and finally we explain how we set the initial distribution of protoplanetary embryos in the system. We also performed simulations testing the effects of stochastic forcing from turbulent fluctuations in the disk.

To perform our simulations we use an adapted version of Mercury (Chambers, 1999). In all our simulations, collisions are considered perfect merging events that conserve linear momentum. During the gas disk phase, our simulations adapted the global timestep to reduce the integration time. Every 1000 timesteps, the timestep was re-evaluated and, if necessary, decreased to be at most 1/25th of the orbital period at the perihelion of the innermost planet. While this technique is not strictly symplectic, we saw no difference in outcome when using it (although it significantly sped up the simulations).

2.1 Disk-Model

The initial structure of a protoplanetary disk can be derived from the radial disk temperature, the gas surface density and viscosity profiles. To model the disk’s structure and evolution we incorporated the 1D disk model fits derived by Bitsch et al. (2015) into our numerical integrator. There are two major advantages in using this approach rather than calculating the evolution of a viscous disk. The first one is that these fits have been calibrated from sophisticated 3D hydrodynamical numerical simulations including effects of viscous heating, stellar irradiation and radial diffusion. The second one is that this approach is computationally cheaper –and for the purposes of this work– more versatile and robust than having to solve a 1D disk evolution model to account for the disk evolution (e.g. Coleman & Nelson (2014)).

From the standard parametrized accretion rate on the star, given by

[TABLE]

and the hydrostatic equilibrium equation

[TABLE]

it is straightforward to determine the disk surface density $\Sigma_{gas}$ using the disk temperature profile given in Bitsch et al. (2015). In Eq. (1) $\alpha$ is the dimensionless ${\rm\alpha}$-viscosity (Shakura & Sunyaev, 1973), h is the disk aspect ratio, r is the heliocentric distance, and ${\Omega_{k}=\sqrt{GM_{\odot}/r^{3}}}$ is the keplerian frequency. T is the disk temperature at the midplane, G is the gravitational constant, ${M_{\odot}}$ is the stellar mass, and ${\mu}$ is the mean molecular weight. In all our simulations the central star is one Solar mass and ${\rm\mu=2.3gmol^{-1}}$. The disk age (or alternatively ${\rm\dot{M}_{gas}}$) is approximated by the follow equation from Hartmann et al. (1998) and modified by Bitsch et al. (2015)

[TABLE]

We do not recalculate the disk structure following the fits in Bitsch et al. (2015) at every timestep of the numerical integrator. Instead, we solve the disk structure every 500 years. Since the disk structure changes on a longer timescale this approach does not affect the validity of our conclusions and allows us to save substantial computational time.

Using Eq. 1-3 we determine the disk temperature using the temperature profiles fits from the Appendix of Bitsch et al. (2015) which are given for different disk metallicities and different regimes of accretion onto the star (or ages of the disk). Following Bitsch et al. (2015) our disk opacity is the same used in Bell & Lin (1994). In our fiducial simulations the disk metallicity is set 1% and the disk ${\rm\alpha}$-viscosity is set $\alpha=5.4\times 10^{-3}$. In this work we do not explore the effects of these parameters. However, our disk presents all main characteristics expected for a proto-planetary disk, with an inner edge and temporary outward migration zones. None of the results that we will obtain will be dependent on specific characteristics of this disk (e.g. the specific location of an outward migration zone), so we expect that they are fairly robust. The issue of the disk’s lifetime will be discussed in Section 5.

The gas disk lifetime in our simulations is set to 5.1 Myr. As discussed in (Bitsch et al., 2015), after the gas accretion rate onto the star drops below $10^{-9}$ the gas density becomes so low that the disk can be evaporated in very short timescales. Thus, as stressed in Bitsch et al. (2015) these fits should not be used to track the disk evolution beyond disk ages corresponding to ${\rm\dot{M}_{gas}=10^{-9}~{}M_{\odot}/yr}$. In our simulations, we allow the disk evolve from $t_{disk}=0$ to 5 Myr (${\rm\dot{M}_{gas}=10^{-9}}~{}M_{\odot}/yr$ ) then we freeze the disk structure at 5 Myr and we exponentially decrease the surface density using an e-folding timescale of 10 Kyr. After 100 Kyr (at $t_{disk}=5.1$ Myr) the disk is assumed to instantaneously dissipate. This allows a smooth transition from the gas-disk to the gas-free phase.

As the disk gets older and thinner low mass planets may eventually be able to open a gap in the disk (e.g. Crida et al., 2006). In our simulations, at the very late stages of the disk the disk’s aspect ratio is about $\sim$0.03 near 0.1 AU. For a 10 Earth mass planet in the very inner regions of the disk the gap should not be fully open yet. Thus, for simplicity we have neglected effects of type-II migration in our simulations.

Magnetohydrodynamic simulations of disk-star interaction suggest that a young star with sufficiently strong dipole magnetic field is surrounded by a low-density gas cavity (Romanova et al., 2003; Flock et al., 2017; Bouvier et al., 2007). In this context, the inner edge of the circumstellar disk typically corresponds to the approximate location where the angular velocity of the star equals to the keplerian orbital velocity, and migration should not continue within the cavity except in unusual circumstances (Romanova & Lovelace, 2006). As the stellar spin rate evolves, the location of the inner edge would evolve as well. The observed stellar spin rate is between 1-10 days (e.g. Bouvier, 2013), which suggest that the corotation radius is 0.01-0.2AU. The orbital period distribution of the innermost Kepler planets is consistent with this. We have included this characteristic of disks in our simulations considering that fixing the inner edge of the disk at 0.1 AU is reasonable. In all our simulations the disk extends from 0.1 to 100 AU. At the inner edge of the disk, the surface density is artificially changed to create a planet trap at about $\sim$0.1 AU. This is done by multiplying the surface density by the following rescaling factor:

[TABLE]

2.2 Disk-planet interaction: Type-I migration

Based on the underlying disk profile we calculate type-I migration, eccentricity and orbital eccentricity damping. Our simulations start with planets that migrate in the type-I regime. The negative of the surface density profile and temperature gradients are given by

[TABLE]

Following (Paardekooper et al., 2010, 2011) and assuming a gravitational smoothing length for the planet potential of ${\rm b=0.4h}$, the total torque from the gas experienced by a type-I migrating planet can be expressed by

[TABLE]

where ${\rm\Gamma_{L}}$ is the Lindblad torque and ${\rm\Gamma_{C}}$ is the corotation torque from the gravitational interaction of the planet with the gas flowing around its orbit. The total torque that a planet feels also depends on its orbital eccentricity and inclination (Bitsch & Kley, 2010, 2011; Cossou et al., 2013; Pierens et al., 2013; Fendyke & Nelson, 2014). To account for this we calculate ${\rm\Gamma_{\text{tot}}}$ with two rescaling functions to reduce the Lindblad and corotation torques accordingly to the planet’s eccentricity and orbital inclination (Cresswell & Nelson, 2008; Coleman & Nelson, 2014). The reduction of the Lindblad torque can be expressed as

[TABLE]

where

[TABLE]

The reduction factor of the coorbital torque ${\rm\Delta_{C}}$ is simply given by

[TABLE]

where e and i are the planet orbital eccentricity and inclination, respectively. ${e_{\text{f}}}$ is defined in Fendyke & Nelson (2014) as

[TABLE]

Under the effects of thermal and viscous diffusion the co-orbital torque is written as:

[TABLE]

The formulae for ${\Gamma_{\text{L}}}$, ${\Gamma_{\text{c,hs,baro}}}$, ${\Gamma_{\text{c,lin,baro}}}$, ${\Gamma_{\text{c,hs,ent}}}$, and ${\Gamma_{\text{c,lin,ent}}}$ are

[TABLE]

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

where ${\xi=\beta-(\gamma-1)x}$ is the negative of the entropy slope with $\gamma=$1.4 being the adiabatic index. The scaling torque ${\Gamma_{0}=(q/h)^{2}\Sigma_{\text{gas}}r^{4}\Omega_{k}^{2}}$ is defined at the location of the planet. The planet-star mass ratio is represented by q, h is the disk aspect ratio, ${\Sigma_{gas}}$ is the surface density and ${\Omega_{k}}$ is the planet’s Keplerian orbital frequency.

Thermal and viscous diffusion effects contribute differently to the different components of the coorbital torque. For example, the barotropic part of the coorbital torque is not affected by thermal diffusion while the entropy related part is affected by both thermal and viscous diffusion. The parameter governing viscous saturation is defined by

[TABLE]

where ${x_{\rm s}}$ is the non-dimensional half-width of the horseshoe region,

[TABLE]

The effects of thermal saturation at the planet location are controlled by

[TABLE]

where ${\rm\chi}$ is the thermal diffusion coefficient which reads as

[TABLE]

where $\rho$ is the gas volume density, $\kappa$ is the opacity and $\sigma$ is the Stefan-Boltzmann constant. The other variables are defined before.

Finally we need to set

[TABLE]

to define the effective $\gamma$, used in Eq. 18, as

[TABLE]

The functions F, G, and K in Eq. 11 are given by

[TABLE]

[TABLE]

and

[TABLE]

Note that the $p$ takes the form of $p_{\nu}$ (Eq. 16) or $p_{\chi}$ (Eq. 18) as defined above.

Following Papaloizou & Larwood (2000) and Cresswell & Nelson (2008) we define the planet’s migration timescale as

[TABLE]

With this definition of $t_{m}$ the respective timescale for a planet on circular orbit to reach the star is $t_{m}/2$. In Eq. 26, the quantity L is a planet’s orbital angular momentum and ${\rm\Gamma_{tot}}$ is the type-I torque defined in Eq. 6 .

To account for the effects of eccentricity and inclination damping we follow the classical formalism of Papaloizou & Larwood (2000) and Tanaka & Ward (2004) modified by Cresswell & Nelson (2006, 2008). Eccentricity and inclination damping timescales are given by ${\rm t_{e}}$ and ${\rm t_{i}}$, respectively. They are defined as

[TABLE]

and

[TABLE]

where

[TABLE]

and ${\rm M_{\odot}}$, ${\rm a_{p}}$, ${\rm m_{p}}$, ${\rm i}$, ${\rm and~{}e}$ are the solar mass and the embryo semimajor axis, mass, orbital inclination, and eccentricity, respectively.

Using the previously defined timescales, the artificial accelerations to account for type-I migration, eccentricity and inclination damping included in the equations of motion of the planetary embryos in our simulations are namely

[TABLE]

[TABLE]

and

[TABLE]

where ${\mathbb{k}}$ is the unit vector in the z-direction. Eq. 30-32 are given in Papaloizou & Larwood (2000) and Cresswell & Nelson (2006, 2008).

2.3 A Migration map

Combining our disk model and type-I migration recipe we can build a migration map showing the migration rate and direction as a function of a planet’s semi-major axis and mass.

Figure 1 presents an evolving migration map of our chosen disk. It shows the direction and relative speed of migration of planets on circular and coplanar orbits at different locations within the disk. The direction of migration is represented by the color. A negative torque (reddish to black) implies inward migration while a positive one (orange to light yellow) represents outwards migration. The grey lines represent locations of zero torque. The regions enclosed by the grey lines are areas of positive torque where planets migrate outward. The strong positive torque at about 0.1 AU is a consequence of the imposed planet trap at the disk inner edge (Masset et al., 2006).

3 Fiducial Model

We performed 120 simulations of our fiducial model. Here we included the effects of type-I migration, eccentricity and inclination damping forces as described in Section 2. In later sections we will present simulations that test the effects of turbulence.

3.1 Initial conditions

Our simulations start from a population of 20-30 planetary embryos distributed beyond 5 AU. This inner edge was chosen to approximately correspond to the location of the water snow line at the start of the simulations, for our chosen disk model (Bitsch et al., 2015). The embryos’ initial masses are randomly and uniformly selected in the range from 0.1 to 4.5 Earth masses. The total mass in planetary embryos in each simulation is about 60 Earth masses. Adjacent planetary embryos were initially spaced by $\sim$5 mutual hill radii $R_{H,m}$(Kokubo & Ida, 2000), where

[TABLE]

In Eq. 33 $a_{i}$ and $a_{j}$ are the semi-major axes of the planetary embryos $i$ and $j$, respectively. Analogously, their masses are given by $M_{i}$ and $M_{j}$. In all simulations the time which planetary embryos start to evolve in the disk corresponds to ${t_{disk}}$= 0 yr.

3.2 Fiducial Model: dynamical evolution

Figures 2 and 3 show the dynamical evolution of two characteristic simulations during the 5 Myr gas phase. The gray lines in each panel delimit the outward migration regions shown in Figure 1. The imposed planet trap at the disk inner edge at $\sim$ 0.1 AU is also evident. Embryos migrate inwards and converge to form resonant chains. Resonant planet-pairs often become unstable and collide. As planetary embryos collide and grow (and the disk evolves) some enter the outward migration regions as it is possible to see in the panels corresponding to 0.6 and 0.75 Myr. Planets inside the outward migration region slowly migrate inwards, once the outward migration region also moves inward and shrinks (Lyra et al., 2010). Between 1 and 4 Myr the outward migration regions have evolved such that planets larger than 5 Earth masses migrate inward to the inner regions of the disk. This deposits planets in long resonant chains at the inner edge of the disk. Small planetary embryos –typically smaller than 1 Earth mass migrate very slowly and stay beyond 1 AU. At the end of the gas disk-phase planetary systems exhibits compact resonant configurations with 5-10 planets inside 0.5 AU.

Fig. 4 and 5 show the dynamical evolution during and after the gas-disk phase for the simulations shown in Fig. 2 and 3. Recall that the gas lasts 5.1 Myr and the entire simulations last 100 Myr. The simulations from Figures 4 and 5 behaved in a similar fashion during the gas disk phase, but their later evolutionary paths are contrasting examples. The system from Figure 4 remained stable after the gas disk dissipated, surviving in a long resonant chain containing 6 planets interior to 1 AU. The system from Fig. 5 underwent a series of instabilities that led to collisions and consequently, a planetary system that is dynamically less compact and more dynamically excited than the one from Fig. 4. The final orbital eccentricities of planets shown in Figure 4 are less than 0.05 and their orbital inclinations are smaller than 0.1 degree. The final eccentricities of planets in Figure 5 are as high as 0.1 and their inclinations are as high as 6 degrees. The most massive planets were somewhat larger in the unstable simulation as well ($\sim 18{\rm\,M_{\oplus}}$ for the simulation from Fig. 5 vs. $11{\rm\,M_{\oplus}}$ for the simulation from Fig. 4).

4 Turbulent Disks

Hydrodynamical instabilities in gaseous protoplanetary disk generate turbulence in the disk and transport angular momentum. Instabilities include the Rossby-wave instability (Lovelace et al., 1999), the global baroclinic instability (Klahr & Bodenheimer, 2003), the Kelvin–Helmholtz instability generated during dust vertical sedimentation towards the disk midplane (Johansen et al., 2006), and the vertical shear instability (Nelson et al., 2013; Stoll & Kley, 2014). Another potentially important source for the observed gas accretion rate on young stars is turbulence driven by the magnetorotational instability (Balbus & Hawley, 1998, MRI). In a sufficiently ionized and magnetized disk the MRI generates magnetohydrodynamic turbulence which leads to outward angular momentum transport (Brandenburg et al., 1995; Armitage, 1998). MRI turbulence produces large-scale, axisymetric and long-lived density and pressure perturbations in the disk (Hawley et al., 1996). Here, we perform simulations testing the role of this kind of turbulence for the formation of close-in super-Earths by inward migration. We assume that turbulence operates at levels consistent with estimates from 3D magnetohydrodynamic simulations  (e.g. Laughlin et al., 2004; Baruteau & Lin, 2010). However, one should also note that recent studies have proposed that a number of non-ideal effects can suppress magneto-rotational turbulence in a large region of the disk (Turner et al., 2014).

The effects of stochastic forcing on planet migration has been studied by several authors (Nelson & Papaloizou, 2003; Papaloizou & Nelson, 2003; Winters et al., 2003; Laughlin et al., 2004; Nelson, 2005; Ogihara et al., 2007; Adams et al., 2008; Rein & Papaloizou, 2009; Lecoanet et al., 2009; Nelson & Gressel, 2010; Baruteau & Lin, 2010; Pierens et al., 2011; Ketchum et al., 2011; Horn et al., 2012; Pierens et al., 2012; Rein, 2012). However, the effects of turbulence for the origins of close-in super-Earths remains to be carefully addressed.

To model the turbulent motion of the gas, which essentially corresponds to density fluctuations in the gaseous disk, we use the model by (Laughlin et al., 2004) as modified by Ogihara et al. (2007) and Baruteau & Lin (2010). The specific turbulent force is given by

[TABLE]

where

[TABLE]

Note that our $\Gamma$ is larger than that in Laughlin et al. (2004) by a factor of $\sim$20 as found by (Baruteau & Lin, 2010). The potential induced by the turbulent perturbations consists of the sum of N independent, scaled wave-like modes as

[TABLE]

where a single oscillation mode is defined as

[TABLE]

In Eq. 35, according to Baruteau & Lin (2010) the strength of the potential is given by

[TABLE]

Also, as noted in Eq. 36, a specific mode is determined by the azimuthal wavenumber $m$, which we randomly sort with a logarithmic distribution between 1 and 96, the center of its initial radial location $r_{c}$, and azimuthal phase $\phi_{c}$. Only modes with wavenumber $m$ smaller than 6 are considered (Ogihara et al., 2007). To sample $r_{c}$ we use a log-normal distribution to select values between ${\rm r_{in}}$ and ${\rm r_{out}}$. In our simulations ${\rm r_{in}=0.1AU}$ and ${\rm r_{out}=25AU}$. The outer edge of the turbulent region (${\rm r_{out}}$) roughly corresponds to the initial location of the outermost planetary embryo of the system. Thus, our planet formation/migration region is significantly wider than that in Baruteau & Lin (2010) and Ogihara et al. (2007). Because of this we assume the existence of a larger number of wave-like modes N in the disk than most of previous studies. We set ${\rm N=125}$, similar to Horn et al. (2012). $\phi_{c}$ is sorted with a uniform distribution between 0 and 2$\pi$. The dimensionless parameter $\xi$ is sorted according to a Gaussian distribution with a unitary standard deviation and mean-value zero. Each mode has an radial extent of $\sigma=\pi r_{c}/4m$. The planet coordinates are represented by the radial distance $r$ and its azimuthal coordinate $\theta$. $\Omega_{c}$ is the Keplerian frequency calculated at $r_{c}$. Still, turbulent fluctuations appear and disappear in the disk. To account for this phenomenon in Eq. 35, the wave mode lifetime is defined by $\Delta t=0.2\pi r_{c}/mc_{s}$ (Baruteau & Lin, 2010), where $c_{s}$ is the local sound speed. Thus, a given mode $m$ may only exist from its birth time $t=t_{0}$ to $\tilde{t}=t-t_{0}~{}=~{}\Delta t$. If $\tilde{t}>\Delta t$ then a new mode is created to take the extinguished one’s place.

Finally, according to Ogihara et al. (2007) the radial, azimuthal and vertical components of the artificial force to account for the effects of turbulence may be written as

[TABLE]

[TABLE]

and

[TABLE]

For simplicity, in our model we assume that the there is no feedback from the stochastic density/pressure fluctuations in the disk on our disk structure model. Basically, we assume the same underlying disk model than that used for our fiducial simulations.

4.1 Simulations

We performed 120 simulations with initial conditions identical to those from our 120 fiducial simulations but with the effects of turbulence. To illustrate the dynamical evolution of the planets in one characteristic simulation including turbulent effects we used the same initial distribution of planetary embryos as in Figure 4. The result of this simulation is shown in Figure 6.

Figure 6 shows that turbulence is important in the outermost parts of the disk, visible by the “random walks” in semi-major axis of planets beyond $\sim$5 AU. Turbulence also tends to increase the orbital eccentricity of planets in this region (Ogihara et al., 2007; Rein, 2012). However, as planets migrate inward the effects of turbulence weakens. When planets reach the inner edge of the disk the turbulence is essentially negligible. It is easy to understand this result by inspecting our turbulent model. In our model, the extent of the fluctuation mode scales with $r_{c}$ (radial coordinate of the mode). In addition, assuming an aspect ratio of 3% for the inner regions of the gaseous disk and a given $m$, the mode lifetime scales as $\Delta t\sim r_{c}^{\frac{3}{2}}$. Thus, modes in the inner regions of the disk (e.g. inside 1 AU) are also short-lived and this implies that they have modest to negligible contributions in disturbing the orbit of those planets (see also discussion on the effects of using large wavenumbers in Ogihara et al. (2007)). Modes generated farther out, on the other hand, may be excessively far from planets already reaching the inner regions of the disk, to be able to strongly interact with them. In addition, it is also important to recall that the turbulence strength scales with the local gas surface density. At one given location, as the disk dissipates its effects becomes relatively weaker. There is thus no significant difference between the results of our fiducial model and the results from that including effects of turbulence. Section 7 will discuss in some details why our results are different from others (Adams et al., 2008; Rein, 2012) in this respect.

In the next two Sections we perform a careful statistical analysis of the results of our fiducial and turbulent models. We will interpret the origins of their eventual differences and compare these results with observations and other works in the literature.

5 Simulation outcomes

Our simulations follow a bifurcated evolutionary path. During the gas disk phase, planetary embryos grow and migrate inward to the inner edge of the disk. The planets settle into long chains of mean motion resonances. It is important to recall that the resonant chains are established very early, so a disk with a reduced lifetime would not help in avoiding these configurations. Resonant chains are typically established in less than 1.5 Myr (e.g. Figures 4-6). After the gas disk dissipates, a large fraction of super-Earth systems undergo a dynamical instability. The planets’ orbits cross, leading to a phase of collisions that destroys the resonant chain. However, a fraction of resonant chains remain stable and never undergo a late phase of collisions.

We now analyze more in detail the evolution and outcome of simulations. We present both the fiducial and turbulent sets of simulations.

Figure 7 shows the cumulative distribution of the time of the last collision in our simulations during two different epochs: from 0 to 5 Myr and from 5 to 100 Myr. During the gas disk phase most systems have their last collisions after 1 Myr (left panel). This is expected since this corresponds to when most planets are approaching the disk inner edge (see for example Figs 5 and 6), and have already reached relatively more compact configurations. This generates dynamical instabilities and collisions. Figure 7 also shows that during the gas disk phase the cumulative distributions of our fiducial and turbulent models are similar. The cumulative distributions of the last collision epoch grows broadly at constant rate from 1 to 4.5 Myr, both in our fiducial and turbulent model simulations. Thus, from 1.5 to 5 Myr there is no preferential time for last collision take place and every forming planetary system exhibit at least a few collisions during the gas disk phase. Because the disk is still present, the eccentricities and inclinations are damped again after each collision and the resonant chain is recovered due to the effect of residual migration.

After gas dispersal most late collisions occur during the first 20 Myr (right panel of Fig. 7). Dynamical instabilities events tend to start as soon as the gas goes away (or when the gas becomes sufficiently rarefied). After about 20 Myr the rate of collisions start to drop. Also, the cumulative distributions do not reach 1. This indicates that not all planetary systems underwent instabilities after the gas disk phase (or did not experience any collision). About 60% of the fiducial simulations and 50% of the turbulent simulations experienced dynamical instabilities and at least one collision after the gas disk phase. Comparatively, only 20% of the in-situ formation simulations of (Ogihara et al., 2015a) were dynamically unstable after gas dispersal. This is probably a consequence of the very small final number of planets in their planetary systems. We also recognize that the fraction of simulations presenting dynamical instabilities may increase if our simulations were integrated for longer than 100 Myr. We are limited in the sense of extending the integration time of these simulations to Gyr timescales because of the very small timestep necessary to resolve the orbits of planets that reach the inner orbit of the disk and the consequent long CPU time demanded. However, we do not expect a linear growth with time of the number of unstable systems.

Figure 8 presents representative simulated systems at 5 and 100 Myr. For reference, we also show selected observed planetary systems. Lines connect planets belonging to the same planetary system. More massive planets tend to park preferentially at the inner edge of the disk (Figure 8, left panel). This is because more massive planets simply migrate faster than small ones and just scatter outward or collide with small ones during their radial excursion to the innermost regions of the disk (Izidoro et al., 2014). However, we stress that there is no dramatic mass-ranking in our model, in contrast with systems of super-Earths produced by in-situ accretion (Ogihara et al., 2015a). This is a strong argument that favours the migration model over the in-situ model. After the gas dissipates, instabilities largely erase the mass-ranking (Figure 8 - middle panel).

5.1 Architecture of planetary systems before gas dispersal

Our simulations can be separated into two groups: those that underwent a late dynamical instability after the gas disk dissipated, and those that remained dynamically stable. We often refer to these as stable and unstable systems. It is important to keep in mind that the unstable systems started out as resonant chains. Thus, it is worth first to investigate how unstable planetary systems compare to stables ones before the gas dispersal, namely at 5 Myr. The question we want to address is the following: Is there any systematic difference in the architecture of unstable and stable systems before gas dispersal? To answer this question we now separate the simulations of our fiducial and turbulent models in stable and unstable groups. Thus, we are left with four sets of simulations which we naturally nominate: fiducial-stable, fiducial-unstable, turbulent-stable and turbulent-unstable. The results of our analysis for each of these groups are shown in Figure 9.

We first focus on the results of our fiducial model. Observing the period ratio distribution of adjacent planet pairs in Figure 9 (top-plot) it is clear that at 5 Myr planetary systems are found in compact first order mean motion resonances (seen as the vertical lines in the plot). Yet, a more clinical analysis reveals that unstable planetary systems are slightly more dynamically compact than stable ones (at 5 Myr). Figure 9 (middle-plot) also shows that at 5 Myr planetary systems have typically multiple planets. Interestingly, the typical number of planets of stable chains is smaller than the number of planets of unstable ones. Figure 9 also show how the masses of planets in stable and unstable system compare to each other at 5 Myr. The mass cumulative distributions of Figure 9 (bottom-plot) show that planets in stable systems are predominantly more massive than those in unstable simulations. Finally, note that these trends are robust since they are observed also in the turbulent sets of simulations.

The results of Figure 9 are quite intuitive since we expect more compact systems to be more prone to exhibit dynamical instabilities than more spread out ones. However, it is important to recall that although there are some quantitative differences, planetary systems shown in Figure 9 are qualitatively similar in some key aspects. The most notable one is that, before gas dispersal, planets pairs are essentially found in compact resonant configurations. However, this is about to change for unstable systems. While stable systems configurations remains essentially unchanged during 100 Myr, in the next Section we show how dynamical instabilities sculpt unstable planetary systems.

5.2 The importance of dynamical instabilities

In this section we analyse the architecture of unstable planetary systems at two epochs: before gas dispersal (at 5Myr) or, equivalently, “before instability” and at 100 Myr (“after instability”). Our main goal is to compare how dynamical instabilities shape the architecture of planetary systems. Note that the dynamical architectures of stable systems are essentially identical at 5 and 100 Myr.

Figure 10 shows the number of planets in unstable planetary systems before and after gas dispersal, or equivalently before and after the dynamical instability phase. Given how compact resonant chains are, they typically contain 6-10 planets inside 0.5 AU. However, given that instabilities spread out the system and reduce the number of planets (due to collisions), the typical unstable system contains 2-5 planets inside 0.5 AU.

There is no glaring difference between the populations of resonant chains in the fiducial and turbulent simulations. Yet roughly 10% more fiducial simulations were unstable at later times (Figure 7). This difference probably comes from a subtle difference in some characteristic of the resonant chains (e.g. number of planets in the chain, commensurabilities and masses of planet-pairs (Matsumoto et al., 2012), amplitude of libration of resonant angles (Adams et al., 2008), etc). Perhaps part of this different is also due to small number statistics.

Figure 11 shows the period ratio distribution in systems that underwent dynamical instabilities after the gas dispersed. We also add the sample of planet candidates from the Kepler mission (Borucki et al., 2011; Batalha et al., 2013; Rowe et al., 2014). To better compare our simulations with observations we applied a simple filter to the observations and our simulations. We only included planets with orbital radii smaller than 0.5 AU. In the Kepler data we included planets with orbital period shorter than 130 days and radii smaller than 4 Earth radii. The motivation for choosing these cutoffs comes from the completeness of Kepler observations (Silburt et al., 2015). We advance to the reader that to account for the effects of inclination distribution we will perform simulated observations in Section 6.

At the end of the gas disk phase at 5 Myr, planets are universally found in chains of mean motion resonances, seen as the vertical jumps in Fig. 11. Resonant chains are far more compact than the observed systems. These are the “before instability” systems (see also the stable systems in Figure 9). There is little difference between the fiducial and turbulent simulations. Instabilities break the resonant chains created during the gas disk phase, promote collisions and scattering events that reduce the number of planets in the system (Figure 10). This tends to produce planetary systems with planets far more apart from each other, and also on orbits with higher eccentricities and orbital inclinations (compare Figures 4 and 5). This same trend is observed in the turbulent simulations. After dynamical instabilities, our planetary systems are less compact than the observed systems, at least for period ratios smaller than 3111The tail of the period ratio distribution probably suffer from selection bias effects.. This suggest that dynamical instabilities play a crucial role in sculpting system of super-Earths. (see also Pu & Wu, 2015, who proposed that the Kepler systems were sculpted by instabilities but without starting from resonant chains)

Figure 12 shows the cumulative distribution of semi-major axis, mass, orbital eccentricity and mutual inclination of simulated planets inside 0.5 AU. Cumulative distributions of semi-major axis (left-hand upper panel of Figure 12) are broadly identical before and after instabilities. Nevertheless, it should be natural to expect that, before dynamical instabilities, simulations have a much smaller fraction of planets far out than afterwards, where planets have been scattered everywhere. However, this effect is not evident in the cumulative distributions because it is normalized to account for planets only inside 0.5 AU. The only difference is that the, before instability, systems show a pileup of planets at about 0.15 AU, near the disk inner edge, whereas this edge has been wiped out in the after instability systems. There is no significant difference between fiducial and turbulent simulations.

The top right panel of Fig. 12 shows the cumulative mass distributions. Planets before instability are far less massive than afterwards. The median planet mass after instability is about 3 times larger. This is simply due to the fact that the unstable systems underwent a late phase of collisions222 Note that, in the inner regions, collisions are much more common than ejections. It is an issue of Safronov number (Safronov, 1972).. There is again no difference between turbulent and fiducial simulations.

The bottom panels of Fig. 12 show the eccentricity and mutual inclination distributions of simulated planets. The orbital distribution of Kepler planets was derived in a series of statistical studies (Lissauer et al., 2011a; Kane et al., 2012; Tremaine & Dong, 2012; Figueira et al., 2012; Fabrycky et al., 2014; Plavchan et al., 2014; Ballard & Johnson, 2014; Van Eylen & Albrecht, 2015). To represent the eccentricity distribution of observations derived from statistical studies we used a Rayleigh distribution with with $\sigma_{e}=.1$ (e.g. Moorhead et al. (2011)). For the inclination distribution we used a Rayleigh distribution with $\sigma_{i}=1.5^{\circ}$ (e.g. Fang & Margot (2012)).

The eccentricities and inclinations of planets before instability (or equivalently of planets in stable systems) are extremely low. These resonant chains have very low eccentricities and are close to perfectly coplanar. As expected, the eccentricity and inclination distributions of the planets are strongly affected by dynamical instabilities. Compared to the eccentricity distribution inferred from statistical analysis, planets after instabilities are in better agreement with the assumed values but the difference is still noticeable (Fig. 12; bottom-left). Of course, we have used a single distribution to represent the expected quantities while the real data may require more than a single distribution to fit the data (Lissauer et al., 2011a). The inclination distribution of planets in our simulations is also quite different from the suggested by previous studies (Fig. 12; bottom-right). In Section 6.2, we will discuss how the inclination distribution of planets in our simulations compare to the distribution inferred from statistical analysis.

6 Matching the observed Kepler planets

In our simulations super-Earth systems follow a typical evolutionary path. Planets grow, migrate inward and pile up into resonant chains far more compact than the observed ones (see Fig. 11). A substantial fraction of these resonant chains become unstable, causing their planetary systems to spread out dynamically.

Following the previous section, we divide our simulations into two batches: those that did not undergo instabilities after the dissipation of the disk and those that did. Stable simulations remain in compact resonant chains whereas the unstable systems have undergone a late phase of accretion and spread out considerably. Note that for this analysis we are using both the fiducial set of simulations and the set that included turbulence, given that we could find no significant difference in outcomes of these simulations (Section 5).

6.1 The period ratio distribution

Figure 11 shows the period ratio distribution of adjacent planet pairs in unstable systems before and after dynamical instabilities. Before instabilities planetary systems are clearly far more compact than the observed Kepler systems. Similarly to planet pairs of stable systems, before instability planet pairs are locked in resonant chains, seen as vertical lines in Figure 11 (see Figure 9 for stable systems). After instabilities, the unstable systems have spread out compared with those that remain stable, and are no longer preferentially found in resonance. In fact, they are are modestly more spread out than observed Kepler systems.

We now test the effect of planet mass on the planets’ orbital spacing. We divide our simulations (unstable systems) into three groups by total mass $M_{tot}=M_{1}+M_{2}$, where $M_{1}$ and $M_{2}$ are the masses of adjacent planets pairs. The low-mass planet pairs have $M_{tot}<25{\rm\,M_{\oplus}}$, the high-mass pairs have $M_{tot}>34.5{\rm\,M_{\oplus}}$, and the medium-mass pairs lie in between. These boundaries were chosen to put a roughly equal number of planet pairs in each bin.

Figure 13 (right panel) shows that, after instabilities, planets are spaced by mutual Hill radii and not by period ratio. Higher-mass planet pairs are systematically more widely-spaced than lower-mass pairs in terms of the period ratio of adjacent planets (left panel). However, pairs of planets with different masses have very similar distributions when measured in mutual Hill radii. They also provide a good match to the Kepler systems, for which we assigned masses using the probabilistic mass-radius relation of Wolfgang et al. (2016): $M=2.7\left(R/R_{\oplus}\right)^{1.3}$. This fits nicely with the results of Fang & Margot (2012), who inferred that Kepler systems are typically spaced by roughly $20\pm 10$ mutual Hill radii. It also emphasizes that the Kepler systems appear to be the result of dynamical instabilities (see also Cossou et al., 2014; Pu & Wu, 2015).

We have shown that the unstable planet pairs are more spread out than the observed ones, and that the stable planet pairs are more compact than the observed ones (see Fig. 11). It follows that a mixture of the two populations may match observations. Indeed, several chains of three or more resonant planets have been identified, such as Kepler-223 (Mills et al., 2016), Kepler-80 (MacDonald et al., 2016), GJ 876 (Rivera et al., 2010) and TRAPPIST-1 (Gillon et al., 2017). While dynamical instabilities can generate resonances (Raymond et al., 2008), the delicate architecture of resonant chains indicate that they are signposts of stable systems. Matching observations thus requires that a fraction of resonant chains remain stable after dissipation of the gas disk.

We attempt to match the observed period ratio distribution using a mixture of stable and unstable systems from our simulations. There are two challenges in this exercise. First, our simulations provide planet masses but most Kepler observations only provide planetary radii. Several studies have used mass constraints to produce mass-radius relationships for transiting planets (e.g. Lissauer et al. (2011b); Fang & Margot (2012); Weiss & Marcy (2014)). As above, we adopt the probabilistic study of Wolfgang et al. (2016). We find that most of our simulated planets are significantly more massive than those inferred for the Kepler systems. As we impose a cutoff of $4R_{\oplus}$ in the Kepler data, the maximum mass inferred using the Wolfgang mass-radius relationship is $16.4{\rm\,M_{\oplus}}$ and the median planet mass in the sample is $6.3{\rm\,M_{\oplus}}$. The median mass of unstable planets at the end of our simulations is $14{\rm\,M_{\oplus}}$.333It’s interesting to note that the median mass of the planets trapped in resonant chains is $5.6{\rm\,M_{\oplus}}$. This is comparable to the inferred value for Kepler sample of $6.3{\rm\,M_{\oplus}}$. One might also speculate that the accretion of planets during late instabilities is inefficient due to collisional erosion and mass loss (see e.g. Leinhardt & Stewart, 2012; Stewart & Leinhardt, 2012). That remains an interesting avenue for future study.

The second challenge in this exercise is understanding observational bias. The transit probability scales linearly with the orbital radius such that close-in planets are far easier to detect (see Charbonneau et al., 2007; Winn, 2010; Wright & Gaudi, 2013). For a perfectly coplanar planetary system aligned with the observer, every planet transits, although it remains a higher probability than the outer planets will be missed. But for a system with many planets on strongly inclined orbits relative to the observer, none transits. In a three planet system, if the middle planet does not transit but the inner and outer planets do, the inferred period ratio is $P_{3}/P_{1}$ rather than $P_{2}/P_{1}$ or $P_{3}/P_{2}$, and is pushed to a much higher value. This spreading to higher period ratios is a function of the mutual inclinations among planets.

We attempt to quantify observational selection effects by performing simulated observations of our planet pairs. We wrote a simple code to observe each of our simulated systems from a large number of vantage points evenly spaced on the celestial sphere. For a given line of sight we kept track of which planets were detected. We then assembled the detected planetary systems from all line of sights.

Figure 14 (left panel) shows the effect of observational bias on our simulated systems. The very low-inclination resonant stable systems were barely affected by observational bias. The plane of these systems is so thin that they are almost always either detected or not (depending on the line-of-sight). The higher-inclination unstable systems (e.g. compare Figures 4 and 5) are more strongly affected by observational bias. Naturally, the orbital configuration of adjacent planet pairs in these systems are systematically shifted to larger period ratios, as planets are eventually missed in transit.

We performed a simple experiment to determine the best-fit combination of simulations to match observations. We tested the effect of the mixing ratio of stable and unstable systems on how well the simulations match observations. We restricted ourselves to planet pairs in which each planet was less massive than the Kepler cutoff of $16.4{\rm\,M_{\oplus}}$, which corresponds to the mass of a $4R_{\oplus}$ planet with the Wolfgang et al. (2016) mass-radius relation. This is our upper size cutoff for Kepler planets. We also restricted this analysis to systems with $P_{2}/P_{1}\leq 3$, to reduce the contribution from systems with hard-to-characterize missed planets. We then generated different samples of simulated planets by varying the fraction of stable systems $F$ included in the sample. We performed K-S tests to roughly judge the goodness of fit for each sample. We note that even though our simulated observations contained thousands of planet pairs, the effective number of points used to generate $p$ values was limited by a combination of the number of simulated planet pairs with appropriate masses and the Kepler sample.

Figure 14 (right panel) shows that our simulated observations roughly match the Kepler sample if less than 25% of planet pairs come from stable simulations (for a probability $p\geq 10\%$ that the two samples are consistent with having been drawn from the same distribution). This can be interpreted as an indicator of the fraction of observed systems that did not undergo an instability, for which the planets survived in a resonant chain. The abundance of resonant chains among known Kepler systems is perhaps 5% (Fabrycky et al., 2014). We expect this to correspond to the contribution from stable systems. Our simulations are indeed consistent with that value.

We expect that this analysis was further affected by the fact that our simulated planets were in general far more massive than the Kepler planets. After an instability, the mutual inclinations between the planets’ orbits naturally correlates with the strength of the gravitational scattering, i.e., the planets’ masses. Pairs of low-mass planets in systems with other, high-mass planets, may therefore have higher mutual inclinations than they otherwise would. This “inclination inflation” should have the effect of pushing period ratios to higher values. We expect that modestly lower-mass systems would therefore have lower mutual inclinations compared with our simulations. Lower-mass systems would appear more compact and thus require a smaller contribution of stable systems to match observations.

To test the other extreme, we performed the same exercise as above without taking observational bias into account. For the same restrictions as above ($M<16.4{\rm\,M_{\oplus}}$, $P_{2}/P_{1}\leq 3$), the best match to observations was for the smallest contribution from stable systems. However, up to nearly 20% of stable systems were allowed while maintaining $p\geq 10\%$.

We conclude that our simulations can indeed provide an adequate match to the period ratio distribution of Kepler planet pairs. We can firmly constrain the abundance of stable resonant chains to contribute less than $\sim 25\%$ of planet pairs.

6.2 The Kepler dichotomy

The Kepler super-Earth sample is bimodal (Lissauer et al., 2011a; Fabrycky et al., 2014): stars tend to either have one or many super-Earths (Fang & Margot, 2012; Ballard & Johnson, 2014). It has been proposed that Kepler dichotomy is a signature of planet-planet scattering (Johansen et al., 2012), in-situ growth close-in (Moriarty & Ballard, 2016), or instabilities produced by spin-orbit (mis-)alignment (Spalding & Batygin, 2016).

We performed synthetic observations to determine whether our simulated planetary systems are consistent with the Kepler dichotomy. We considered viewing angles spanning from 30 degrees above the initial $i=0$ plane to 30 degrees below, with even sampling in azimuth. For each viewing angle we determined the number of planets that transited for each of our stable and unstable simulations.

From viewing angles where at least one transiting planet was detected, Figure 15 shows the distribution of the number of transiting planets. The difference between stable and unstable systems is striking. Stable systems have such low mutual inclinations that it is common to detect high-N systems. Only 18% of detections were of a single planet in transit, and 66% of detections had $N\geq 3$. In contrast, for unstable systems 78% of outcomes were single-planet detections and only 7% of detections had 3 or more planets.

If we combine stable and unstable systems in a 1-to-9 ratio (i.e., with 90% unstable and 10% stable), we naturally obtain a “dichotomy” that is almost identical to the observed Kepler dichotomy (Fig. 15). The significant mutual inclinations in the unstable systems produce a peak at $N=1$, while the very low mutual inclination stable systems provide a long tail to high-N. Low-multiplicity systems $N=1-2$ are dominated by unstable systems whereas high-multiplicity systems ($N\geq 4$) are more often stable. The Kepler-223 (Mills et al., 2016) and TRAPPIST-1 (Gillon et al., 2017) multi-resonant super-Earth systems appear to be good examples of high-N stable systems.

This would suggest that the Kepler dichotomy is simply an observational artifact. Our simulated super-Earth systems naturally produce a spike of apparent singleton planets. However, each of those systems contains at least one – and in some cases many more – additional planets within 0.5 AU and beyond. As we showed above this same sample of simulations matches the observed period ratio distribution (see Figure 14). Of course, we do not claim to match all the details of the Kepler sample, but our results strongly suggest that there is no need to invoke special evolutionary histories for single super-Earth systems. If, on the other hand, singleton super-Earth systems turn out to have a high false positive rate, this analysis would need to be re-visited.

Previous studies had mixed success in matching the Kepler dichotomy. Johansen et al. (2012) drew from the observed period ratio distribution and varyied the width of a Gaussian distribution of mutual inclinations. Starting from triple-planet systems, they were unable to reproduce the dichotomy, quantified as the relative abundance of observed triple-, double-, and single-planet systems. Lissauer et al. (2011a) and Ballard & Johnson (2014) were equally unable to match the Kepler dichotomy simply because their single component-model of planetary architectures underestimates the number of single planet systems.

Building on the work of Johansen et al. (2012), we tested whether we could match the dichotomy with a single inclination distribution (single component-model of system architecture). We generated synthetic planetary systems as follows. The closest planet was placed between 0.05 and 0.1 AU, and subsequent planets were spaced by drawing a period ratio evenly between 1.5 and 3. Systems extended out to 0.5 AU, naturally providing a wide range of planet multiplicities. The planets’ orbital inclinations were drawn from a Rayleigh distributions with $\sigma$ varying from 1 to 10 degrees. The results of this experiment matched qualitatively those of Johansen et al. (2012). Figure 16 (left panel) shows that no Rayleigh distribution matches the Kepler multiplicity distribution. The best result is for $\sigma=4^{\circ}$, which provides an acceptable match for systems with 1, 2, or 3 planets. However, the Rayleigh distribution dramatically underproduces systems with $N>3$.

The reason our simulations match the dichotomy can be inferred from their mutual inclination distribution (Figure 16, right panel). The small contribution of near-coplanar stable systems provides large $N$ systems whereas the broad inclination distribution of unstable systems creates a peak at low $N$.

Fang & Margot (2012) were able to match the dichotomy by introducing an additional parameter: the multiplicity distribution. This extra parameter renders the problem much simpler. For instance, consider perfectly co-planar planetary systems. The observed distribution can be retrieved in a straightforward way if the multiplicity distribution of these systems matches the observed distribution (to within a small observational bias). Indeed, assuming certain statistical distributions for the number of planets and mutual orbital inclinations of planet pairs, Fang & Margot (2012) were able to match both the number of multiple and single-planet Kepler systems for some combinations of parameters.

Hansen & Murray (2013) performed in-situ growth simulations starting from a set of initial conditions (planetary embryos) reflecting a putative radial distribution of mass in solids. Similar to previous studies, they found that the number of single planet systems are more common in the Kepler data than in their simulated population. In fact, to match the Kepler multiplicity distribution in-situ growth simulations also require a very specific mix of protoplanetary systems produced from simulations with quite different initial conditions. Simulations starting with a distribution of protoplanetary embryos derived from a shallow disk can account for the observed single planet population while simulations starting with sufficiently steep radial mass distributions of solids tend to produce multiple transiting planets. A mix of planetary systems produced from steep and shallow disks allows one to build a good match to the Kepler planet multiplicity (Moriarty & Ballard, 2016). Indeed, the radial distribution of mass in solids in real protoplanetary disks may vary from disk to disk but in-situ accretion simulations also predict that the radial distribution of planets in multi-planet systems should be mass ranked. More specifically, they predict that planet mass should decrease with semi-major axis (Ogihara et al., 2015a; Moriarty & Ballard, 2016). This expected signature seems to be at odds with the Kepler data.

Our simulations match the observed Kepler dichotomy with no free parameters or ad-hoc assumptions. Even the relative abundance of stable and unstable systems was chosen to be the one that already matches the observed period ratio distribution. The mutual inclination distribution in our simulations is clearly not a simple function like the ones used in statistical studies. We encourage statistical studies to connect with dynamical models to include more realistic inputs.

There is one significant difference between the population of planets that matches the dichotomy (Fig. 15) and the one that matches the period ratio distribution (Fig. 14). In constructing our sample to match the period ratio distribution, we selected the lowest-mass planet pairs from our simulations because those were the most reasonable match to the observed planets. However, in each simulated system there is generally a mix of planet pairs of different masses, so we could not select just the lower-mass pairs when performing simulated observations. However, the inclinations in unstable systems are generated by dynamical instabilities after gas dispersal. This is a mass-dependent process, as the degree of excitation scales with the planets’ escape velocities. We do find lower mutual inclinations for lower-mass planet pairs. For whole systems of lower-mass planets one might expect systematically lower mutual inclinations. This would tend to create a smaller peak at $N=1$. While we have not performed this exercise (because a new set of simulations would be required), we can infer the expected outcome. We expect our simulations to be able to match the multiplicity distribution for $N\geq 2$ but underestimate the Kepler peak at $N=1$. This would lead us to predict that the observed Kepler singletons include a certain percentage of planets that we cannot explain. These might be false positives or perhaps true singletons generated by another mechanism (see, e.g., Izidoro et al., 2015).

7 Comments on different disk turbulence effects in different works

To compare how the effects of turbulence affect the outcome of our simulations we calculated a diversity of parameters at two different ages of our simulations. This was done first at the disappearance of the disk and then after 100 Myr of integration in a gas free scenario. We found that essentially there are no differences between the results of simulations assuming turbulent or non-turbulent disk, both at the end of the disk lifetime and in the aftermath of the planet instability. Here we compare our results with other previously published in the literature.

Adams et al. (2008) studied the effects of disk turbulence in the context of gas giant exoplanets using a stochastic pendulum model. These authors conclude that resonances should be rare in turbulent systems. This result is different from what we find here since our fiducial and turbulent models produced very similar results. So let’s interpret the origin of this difference. First of all, we have to recall that our study is dedicated to low mass planets while these authors have focused on the effects of disk turbulence for gas giant planets. Yet, we have used quite different approaches. For example, in stark contrast with our simulations, the model by Adams et al. (2008) assumes that the planet pairs are initially already locked in mean motion resonance. More importantly, they neglect the effects of planet migration, gas disk dispersal, tidal damping of orbital eccentricity and inclination while applying the turbulence forces and checking on whether the rogue resonant pair will survive in resonance or not.

Rein (2012) numerically studied the effect of stochastic migration of planets in systems extracted from Kepler data. As transit-surveys only provide the radius but not mass of these planets he used the mass-radius relationship of Fabrycky et al. (2012) to estimate the planetary masses. The stellar mass and planet periods were also taken from the Kepler Objects of Interest (KOI) catalog. Planets’ eccentricities and orbital inclinations were set to zero. To mimic type-I migration and eccentricity damping Rein (2012) used simple migration timescales of about $10^{3}$ to $10^{4}$ years. In his prescription for the stochastic forcing he assumed that the turbulent strength is a small fraction ($\sim 10^{-5}-10^{-6}$) of the gravitational force from the central Rein & Papaloizou (2009).

The main result in Rein (2012) is that it is possible to reproduce the period ratio distribution of close-in super-Earths if the effect of stochastic forcing is included in the simulations during migration. However, the results of our simulations are very different from his. We can understand the differences from at least three sources. Let’s first focus on the gas disk phase.

The first great difference between our model and Rein’s is how synthetic forces were implemented. In Rein (2012)’s simulations only the outermost planet of the system is allowed to type-I migrate and feel the damping of eccentricity and inclination. Yet, in his model all planets feel stochastic forcing. The motivation for this choice is unclear. In our simulations all planets felt stochastic forcing and type-I, eccentricity and inclination damping self-consistently. Another important issue is that Rein’s simulations were only integrated for $10^{4}$ orbital periods of the outermost planet of the system. This integration time may not guarantee capture in resonance for a large fraction of the planet pairs in the Kepler catalog (see Figure 1 in Rein (2012)). Basically, the approach in Rein (2012) strongly favors the production of dynamically relaxed planetary systems. While Rein’s simulations had short integration time, our protoplanetary disks lived for about 5 Myr. Thus, planets in our simulations had the time to migrate to the inner edge of the disk and pile up in long chain of resonances before the gas dissipated (e.g. Figure 6). In our simulations the disk surface density decreases with increasing time according to the disk model derived from hydrodynamical simulations (Bitsch et al., 2015). Finally, the third reason which may explain the difference between the results of these models comes from the stochastic forcing model. We used a more sophisticated prescription than Rein (2012) to mimic the effects of turbulence in the disk of gas. For example, Rein’s simulations scale the stochastic force using simply a fraction of the relative gravitational force from the central star while our model scales the strength of the turbulence taking into account the disk surface density, aspect ratio and the distance between the center of the gas fluctuation and the planet (Laughlin et al., 2004; Ogihara et al., 2007). The stochastic kicks that each planet feels in Rein’s simulations are uncorrelated because their model is based on a first order Markov process. This means that even in the case when planets are very close to each other they may feel very different stochastic kicks (Rein & Papaloizou, 2009). Moreover, the amplitude of the turbulence forcing in his simulations is purposely chosen to produce results consistent with observations (Batygin & Adams, 2017). The model of Laughlin et al. (2004) is more robust in this sense.

It may be easier to understand our results if we compare them with those in Ogihara et al. (2007) where the stochastic forcing model is similar. In our fiducial simulations we set $\alpha=5.4\times 10^{-3}$ (alpha-viscosity parameter), thus the turbulence strength parameter $\gamma$ is about $\sim 2.5\times 10^{-4}$ inside 1 AU (see Eq. 38). Recall that in our model $\gamma$ comes from Baruteau & Lin (2010) where it is obtained from calibration with three-dimensional magnetohydrodynamics calculations. Still, the typical values of $\alpha$ in MRI active zones is in the range of $10^{-3}$ and $10^{-1}$ (Fromang & Papaloizou, 2007) which is in great agreement with our chosen value. Curiously, the simulations in Ogihara et al. (2007) show that the effects of turbulence were only pronounced in cases where $\gamma\geq 10^{-1}$. Smaller $\gamma$ produced results very similar to those where turbulence was not included. We confirm this result with our simulations. However, it is not clear if $\gamma\gtrsim 0.1$ would be consistent with results from MHD simulations of very turbulent disks.

It is also interesting to note that our results seem to agree well with conclusions from a very recent paper by Batygin & Adams (2017) where is derived an analytic criterion for turbulent disruption of mean motion resonances. According to these authors –at the inner regions of the disk– only planet pairs with mass ratios smaller than $(m_{1}+m_{2})/M_{\odot}\lesssim 10^{-5}$ should be susceptible to disruption of resonant configurations. Given that about 90% (50%) of our planet pairs have combined masses larger than 6${\rm\,M_{\oplus}}$ (10${\rm\,M_{\oplus}}$), this fits nicely with our results.

8 Why do so many resonant chains go unstable?

Figure 14 shows that our simulations can match the sample of Kepler planets if at least 75% of resonant chains go unstable and at most 25% remained stable. Among the Kepler multiple-planet systems, only a handful have been characterized as resonant chains (e.g. Mills et al. (2016)). There is an additional excess of planet pairs just exterior to the 3:2 and 2:1 resonances (Fabrycky et al., 2014). All told, it appears that $\sim 5\%$ of Kepler planet pairs are in resonance. If we equate that occurrence rate with the probability of a given system remaining stable in our simulations, it follows that $\sim 5\%$ of resonant chains remain stable. This yields an acceptable match to the Kepler period ratio distribution (Fig. 14).

Only 50-60% of resonant chains were unstable in our fiducial and turbulent sets of simulations (Fig. 7). This is a far cry from the roughly 95% required to match both the period ratio distribution and abundance of observed resonances in the Kepler systems.

How can we explain the deficit of unstable systems in our simulations? Our simulations are idealized and are missing several physical effects. For instance, we have only used a crude description of the conditions at the inner edge of the disk, which is critical in anchoring resonant chains. We did not include the effects of general relativity, which causes close-in orbits to precess. We also did not include tidal interactions between the planets and the star. As young stars may dissipate more strongly than previously thought (Mathis, 2015), it is possible that tides play an important role (Bolmont et al, in preparation). For simplicity, in our simulations the inner edge of the disk remains fixed at $\sim$0.1AU during all disk lifetime. However, in reality, the inner edge of the disk should move as gas the disk and star evolve, due to the balance between the stellar magnetic pressure and the viscous torque of the disk (e.g. Koenigl, 1991). Liu et al. (2017) showed that under certain conditions two planets can experience divergent migration as the disk inner edge moves outward. This mechanism could be also another important inductor of dynamical instabilities.

We have neither considered gas accretion onto planetary cores (e.g. Ginzburg et al., 2016), nor mass loss during planetary collisions (e.g. Inamdar & Schlichting, 2016), nor effects of fast-rotating stars (Spalding & Batygin, 2016). Also, we have not included a reservoir of small planetary bodies, as pebbles or planetesimals, in the disk. Left-over planetesimals in the disk could be another potential trigger of later dynamical instabilities in these systems (Chatterjee & Ford, 2015). We also did not explore the effects of the gas disk viscosity, the initial distribution of planetary embryos and photo-evaporation in the disk. All these are certainly interesting routes for future studies.

We should also recall that because we overestimate the masses of the real Kepler super-Earths, our simulations may simply underestimate the rate of instabilities. It is also possible that our simulations were stopped prematurely. The number of instabilities was still increasing at 100 Myr, even though the rate was decreasing. It is not out of the question that within 5 Gyr all of our resonant chains would have become unstable. Indeed, for a similar-aged parent star a system of close-in super-Earths is dynamically far older than the Earth orbiting the Sun, because the super-Earths have completed 1-2 orders of magnitude more orbital periods. Of course, if resonant chains simply undergo later instabilities, then the abundance of resonances should decrease with the stellar age. Kepler-223, the best-characterized resonant chain to date, is a relatively old star (age of $\sim 6$ Gyr Mills et al., 2016). However, more data is needed to test this idea.

Given the successes of our model, we consider this shortcoming to be important. We believe that we are missing a trigger to explain why so many resonant chains go unstable.

9 Conclusions

We have simulated the migration and growth of system of close-in super Earths in an evolving gaseous disk model (Bitsch et al., 2015). We found that stochastic forcing from disk turbulence have no measurable effect on the growth of close-in super-Earths.

We propose that systems of super-Earths follow a standard evolutionary path. Embryos grow in the outer disk and migrate inward due to torques from the gaseous disk (Figs 1-3). When the first embryo reaches the inner edge of the disk its migration is stopped by the planet disk-edge interaction (Masset et al., 2006) and other embryos migrate into a resonant chain with up to 10 or more close-in planets (Figs. 4-6) . The configuration of resonant chains are far more compact than the observed Kepler systems (Fig. 11). As the gas disk dissipates, about 50-60% of our resonant chains become unstable and undergo a late phase of giant collisions. This spreads their orbits out and spaces them by mutual Hill radii rather than by orbital period ratio (Fig. 13, left panel). Our simulations match the period ratio distribution of the observed Kepler planet pairs if 75-100% of resonant chains go unstable. Taking into account the abundance of resonances among Kepler planets (Fabrycky et al., 2014), we expect that the true instability rate of resonant chains is roughly 95%. Our results also suggest that the large number of detected single planet systems is simply an observational effect due to the mutual orbital inclination of planets in Kepler systems. Thus, although planets have been essentially observed either in single or multiple planet systems our results suggest that this does not necessarily imply the existence of any dichotomy in the architecture of planetary systems.

Finally, when comparing our simulations with the Kepler systems, we are left with a mystery: why does it appear that so many resonant chains go unstable?

Acknowledgements

We thank the anonymous referee for his/her constructive comments that helped us to improve the manuscript. A. I, S. N. R., A. M., A. P. and F. H. thank the Agence Nationale pour la Recherche for support via grant ANR-13-BS05-0003- 01 (project MOJO). Computer time for this study was provided by the computing facilities MCIA (Mésocentre de Calcul Intensif Aquitain) of the Université de Bordeaux and of the Université de Pau et des Pays de l’Adour. A. I thanks also partial financial support from FAPESP (Proc. 16/19556-7 and 16/12686-2) during the final preparation of this paper. B. B. thanks the Knut and Alice Wallenberg Foundation for their financial support.