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

A. Izidoro, Masahiro Ogihara, Sean N. Raymond, Alessandro Morbidelli, A.

Pierens, Bertram Bitsch, C. Cossou, F. Hersant

► To cite this version:

A. Izidoro, Masahiro Ogihara, Sean N. Raymond, Alessandro Morbidelli, A. Pierens, et al.

HAL Id: hal-01489242

https://hal.science/hal-01489242v1

Submitted on 23 Aug 2025

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Breaking the chains: hot super-Earth systems from migration and disruption of compact resonant chains

Andre Izidoro,<fnref n="1" />,<fnref n="2" />* Masahiro Ogihara,<fnref n="3" /> Sean N. Raymond,<fnref n="1" /> Alessandro Morbidelli,<fnref n="4" /> Arnaud Pierens,<fnref n="1" /> Bertram Bitsch,<fnref n="5" /> Christophe Cossou<fnref n="6" /> and Franck Hersant<fnref n="1" />

Accepted 2017 May 17. Received 2017 May 17; in original form 2017 March 9

## ABSTRACT

'Hot super-Earths' (or 'mini-Neptunes') between one and four times Earth's size with period shorter than 100 d orbit 30-50 per cent of Sun-like stars. Their orbital configuration - measured as the period ratio distribution of adjacent planets in multiplanet systems - is a strong constraint for formation models. Here, we use N-body simulations with synthetic forces from an underlying evolving gaseous disc to model the formation and long-term dynamical evolution of super-Earth systems. While the gas disc is present, planetary embryos grow and migrate inward to form a resonant chain anchored at the inner edge of the disc. 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. Disc 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 per cent. 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 ~5 per cent of the cases. This leads to a mystery: in our simulations only 50-60 per cent of resonant chains became unstable, whereas at least 75 per cent (and probably 90-95 per cent) must be unstable to match observations.

**Key words:** methods: numerical – planets and satellites: dynamical evolution and stability – planets and satellites: formation – planet–disc interactions – protoplanetary discs.

## 1 INTRODUCTION

Among the thousands of confirmed exoplanets, hot super-Earths or mini-Neptunes – with radii between 1 and 4 R⊕ (1 < M⊕ < 20), orbiting very close to their host stars – form by far the largest population (Borucki et al. 2010, 2011; Lissauer et al. 2011b; Mayor et al. 2011; Schneider et al. 2011; Fabrycky et al. 2012; Howard et al. 2012; Batalha et al. 2013; Dong & Zhu 2013; Fressin et al. 2013; Howard 2013; Mullally et al. 2015; Petigura, Howard & Marcy 2013). Statistical studies suggest that about one out of three solar-

type stars (FGK spectral types) host a super-Earth with orbital period shorter than 100 d (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 multiplanet 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^\circ$; Lissauer et al. 2011a; Mayor et al. 2011; 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, Barnes & Mandell 2008; Raymond & Morbidelli 2014;

<footnote n="1">¹Laboratoire d'astrophysique de Bordeaux, Univ. Bordeaux, CNRS, B18N, allée Geoffroy Saint-Hilaire, F-33615 Pessac, France</footnote>

<footnote n="2">²UNESP, Univ. Estadual Paulista – Grupo de Dinâmica Orbital & Planetologia, Guaratinguetá, CEP 12516-410 São Paulo, Brazil</footnote>

<footnote n="3">³Division of Theoretical Astronomy, National Astronomical Observatory of Japan, 2-21-1, Osawa, Mitaka, Tokyo 181-8588, Japan</footnote>

<footnote n="4">⁴Laboratoire Lagrange, UMR7293, Université Côte d'Azur, CNRS, Observatoire de la Côte d'Azur, Boulevard de l'Observatoire, F-06304 Nice Cedex 4, France</footnote>

<footnote n="5">⁵Lund Observatory, Department of Astronomy and Theoretical Physics, Lund University, Box 43, SE-22100 Lund, Sweden</footnote>

<footnote n="6">⁶IAS Institut d'Astrophysique Spatiale, Université Paris Sud, Bâtiment 121, F-91405 Orsay, France</footnote>

<footnote>* E-mail: izidoro.costa@gmail.com</footnote>
Morbidelli & Raymond 2016): (1) *in situ* accretion (Raymond et al. 2008; Hansen & Murray 2012, 2013; Chiang & Laughlin 2013; Hansen 2014; Ogihara, Morbidelli & Guillot 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 come in two flavours: (a) standard *in situ* accretion models that invoke high-mass discs 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 disc in consequence of a local concentration of pebbles or small planetesimals drifting inward due to gas drag (Boley & Ford 2013; Boley, Morris & Ford 2014; Chatterjee & Tan 2014, 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; Schlaufman 2014; Schlichting 2014; Chatterjee & Tan 2015; Izidoro et al. 2015; Ogihara et al. 2015a). For example, it ignores the effects of planet–disc gravitational interaction. Planets forming *in situ* grow extremely fast because of short dynamical time-scales and the required abundant amount of mass in the inner regions of the disc (Hansen & Murray 2012, 2013). They tend to reach masses large enough (Hansen & Murray 2013; Hansen 2014) to sufficiently perturb the surrounding gas in time-scales much shorter than the expected lifetime of protoplanetary discs around young stars (Ogihara et al. 2015a). Planet–disc gravitational interaction leads to angular momentum transfer between the planet and the disc (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, Takeuchi & Ward 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–disc 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 disc phase they would not be forming truly *in situ* (Ogihara et al. 2015a).

The drift-then-assembly models (Boley & Ford 2013; Boley et al. 2014; Chatterjee & Tan 2014, 2015; Hu et al. 2016) are quite promising but what actually happens near the inner edge of the disc, 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 disc seems to be generic from theoretical considerations (e.g. Morbidelli et al. 2015). Migration seems to be a generic process as well. Therefore, 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 disc (e.g. Papaloizou & Lin 1984; Crida, Morbidelli & Masset 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 disc properties. Planets may either migrate inward or outward depending on the combination of different torques from the disc (Paardekooper & Mellema 2006, 2008; Baruteau & Masset 2008; Paardekooper & Papaloizou 2008; Kley, Bitsch & Klahr 2009). Indeed, the Lindblad torque tends to push the planet inward (e.g. Ward 1986; Ward 1997) but depending on the planet

mass the gas-flowing in co-orbital 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; Paardekooper, Baruteau & Kley 2011). There are locations within the disc where the net torque is zero (Lyra et al 2010; Horn et al 2012; Bitsch et al. 2013, 2014, 2015; Cossou et al 2013; Pierens et al 2013). However, as the disc dissipates and cools, thermodynamics and viscous effects become less important and planets are released to migrate inward (Lyra, Paardekooper & Mac Low 2010; Horn et al. 2012; Bitsch et al. 2014). It is hard to imagine planets completely escaping inward type-I migration, although disc winds may be a possible explanation of a global suppression of type-I migration in the inner parts of the disc (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 predict that as the disc 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; Horn et al. 2012; Rein 2012; Rein et al. 2012; Haghighipour 2013; Ogihara & Kobayashi 2013; Cossou et al. 2014; Raymond & Cossou 2014; Liu et al. 2015; Ogihara et al. 2015a; Liu, Zhang & Lin 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 disc 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). Secondly, a number of mechanisms have been proposed for accounting for the excess of planet pairs outside mean-motion resonances which could consistently operate with the inward-migration model. Among them are star–planet tidal dissipation (Papaloizou & Terquem 2010; Papaloizou 2011; Delisle et al. 2012; Lithwick & Wu 2012; Batygin & Morbidelli 2013; Delisle & Laskar 2014; Delisle, Laskar & Correia 2014), planet scattering of leftover planetesimals (Chatterjee & Ford 2015), turbulence in the gaseous disc (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 disc (Batygin 2015). Thirdly, and most importantly, planets that form in resonance may not remain in resonance. Rather, planets in resonance can become unstable when the gas disc 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 discs (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 disc of gas. After the gas
disc'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 disc 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 after gas disc 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 discs. We also follow the subsequent phase of dynamical evolution of formed planetary systems after gas disc dissipation. During the gaseous phase, we mimic the effects of the disc of gas on the planets by applying artificial forces on to the planets (or protoplanetary embryos). These forces were calibrated from truly hydrodynamical simulations. In this section we describe our gas disc 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 disc.

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 disc phase, our simulations adapted the global time-step to reduce the integration time. Every 1000 time-steps, the time-step 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 Disc model

The initial structure of a protoplanetary disc can be derived from the radial disc temperature, the gas surface density and viscosity profiles. To model the disc's structure and evolution, we incorporated the 1D disc 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 disc. 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 disc evolution model to account for the disc evolution (e.g. Coleman & Nelson 2014).

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

$$ \dot{M}_{\text{gas}} = 3\pi\alpha h^2 r^2 \Omega_k \Sigma_{\text{gas}}, \quad (1) $$

and the hydrostatic equilibrium equation

$$ T = h^2 \frac{GM_{\odot}}{r} \frac{\mu}{R} \quad (2) $$

it is straightforward to determine the disc surface density $\Sigma_{\text{gas}}$ using the disc temperature profile given in Bitsch et al. (2015). In equation (1) $\alpha$ is the dimensionless $\alpha$-viscosity (Shakura & Sunyaev 1973), $h$ is the disc aspect ratio, $r$ is the heliocentric distance and $\Omega_k = \sqrt{GM_{\odot}/r^3}$ is the Keplerian frequency. $T$ is the disc temperature at the mid-plane, $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 $\mu = 2.3 \text{ g mol}^{-1}$. The disc age (or alternatively $\dot{M}_{\text{gas}}$) is approximated by the follow equation from Hartmann et al. (1998) and modified by Bitsch et al. (2015),

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

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

Using equations (1)–(3), we determine the disc temperature using the temperature profile fits from the appendix of Bitsch et al. (2015) which are given for different disc metallicities and different regimes of accretion on to the star (or ages of the disc). Following Bitsch et al. (2015), our disc opacity is the same used in Bell & Lin (1994). In our fiducial simulations, the disc metallicity is set 1 per cent and the disc $\alpha$-viscosity is set $\alpha = 5.4 \times 10^{-3}$. In this work, we do not explore the effects of these parameters. However, our disc presents all main characteristics expected for a protoplanetary disc, 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 disc (e.g. the specific location of an outward migration zone), so we expect that they are fairly robust. The issue of the disc's lifetime will be discussed in Section 5.

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

As the disc gets older and thinner, low-mass planets may eventually be able to open a gap in the disc (e.g. Crida et al. 2006). In our simulations, at the very late stages of the disc the disc's aspect ratio is about ~0.03 near 0.1 au. For a 10 $M_{\oplus}$ planet in the very inner regions of the disc 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 disc-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; Bouvier et al. 2007; Flock et al. 2017). In this context, the inner edge of the circumstellar disc 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 and 10 d (e.g. Bouvier 2013), which suggest that the corotation radius is 0.01–0.2 au. The orbital period distribution of the innermost Kepler planets is consistent with this. We have included this characteristic of discs in our simulations considering that fixing the inner edge of the disc at 0.1 au is reasonable. In all our simulations, the disc extends from 0.1 to 100 au. At the inner edge of the disc, the surface density is artificially changed to create a planet trap at about ~0.1 au. This is done by multiplying the surface density by the following rescaling factor:

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

## 2.2 Disc–planet interaction: type-I migration

Based on the underlying disc 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

$$ x = -\frac{\partial \ln \Sigma_{\text{gas}}}{\partial \ln r}, \quad \beta = -\frac{\partial \ln T}{\partial \ln r}. \qquad (5) $$

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

$$ \Gamma_{\text{tot}} = \Gamma_{\text{L}}\Delta_{\text{L}} + \Gamma_{\text{C}}\Delta_{\text{C}}, \qquad (6) $$

where $\Gamma_{\text{L}}$ is the Lindblad torque and $\Gamma_{\text{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, Raymond & Pierens 2013; Pierens, Cossou & Raymond 2013; Fendyke & Nelson 2014). To account for this, we calculate $\Gamma_{\text{tot}}$ with two rescaling functions to reduce the Lindblad and corotation torques according to the planet’s eccentricity and orbital inclination (Cresswell & Nelson 2008; Coleman & Nelson 2014). The reduction of the Lindblad torque can be expressed as

$$ \Delta_{\text{L}} = \left[ P_e + \frac{P_e}{|P_e|} \times \left\{ 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 \right\} \right]^{-1}, \qquad (7) $$

where

$$ P_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}. \qquad (8) $$

The reduction factor of the co-orbital torque $\Delta_{\text{C}}$ is simply given by

$$ \Delta_{\text{C}} = \exp\left(\frac{e}{e_f}\right) \left\{1 - \tanh\left(\frac{i}{h}\right)\right\}, \qquad (9) $$

where $e$ and $i$ are the planet orbital eccentricity and inclination, respectively. $e_f$ is defined in Fendyke & Nelson (2014) as

$$ e_f = 0.5h + 0.01. \qquad (10) $$

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

$$ \begin{aligned} \Gamma_{\text{C}} ={}& \Gamma_{\text{c,hs,baro}} F(p_{\text{v}})G(p_{\text{v}}) + (1 - K(p_{\text{v}}))\Gamma_{\text{c,lin,baro}} \\ & + \Gamma_{\text{c,hs,ent}} F(p_{\text{v}})F(p_{\chi})\sqrt{G(p_{\text{v}})G(p_{\chi})} \\ & + \sqrt{(1 - K(p_{\text{v}}))(1 - K(p_{\chi}))}\Gamma_{\text{c,lin,ent}}. \end{aligned} \qquad (11) $$

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

$$ \Gamma_{\text{L}} = (-2.5 - 1.5\beta + 0.1x) \frac{\Gamma_0}{\gamma_{\text{eff}}}, \qquad (12) $$

$$ \Gamma_{\text{c,hs,baro}} = 1.1 \left(\frac{3}{2} - x\right) \frac{\Gamma_0}{\gamma_{\text{eff}}}, \qquad (13) $$

$$ \Gamma_{\text{c,lin,baro}} = 0.7 \left(\frac{3}{2} - x\right) \frac{\Gamma_0}{\gamma_{\text{eff}}}, \qquad (14) $$

$$ \Gamma_{\text{c,hs,ent}} = 7.9\xi \frac{\Gamma_0}{\gamma_{\text{eff}}^2}, \qquad (15) $$

and

$$ \Gamma_{\text{c,lin,ent}} = \left(2.2 - \frac{1.4}{\gamma_{\text{eff}}}\right) \xi \frac{\Gamma_0}{\gamma_{\text{eff}}}, \qquad (16) $$

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 disc aspect ratio, $\Sigma_{\text{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 co-orbital torque. For example, the barotropic part of the co-orbital torque is not affected by thermal diffusion while the entropy related part is affected by both thermal and viscous diffusions. The parameter governing viscous saturation is defined by

$$ p_{\nu} = \frac{2}{3} \sqrt{\frac{r^2 \Omega_k}{2\pi\nu}} x_s^3, \qquad (17) $$

where $x_s$ is the non-dimensional half-width of the horseshoe region,

$$ x_s = \frac{1.1}{\gamma_{\text{eff}}^{1/4}} \sqrt{\frac{q}{h}}. \qquad (18) $$

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

$$ p_{\chi} = \frac{2}{3} \sqrt{\frac{r^2 \Omega_k}{2\pi \chi}} x_s^3, \qquad (19) $$

where $\chi$ is the thermal diffusion coefficient that reads as

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

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

$$ Q = \frac{2\chi}{3h^3r^2\Omega_k} \qquad (21) $$

to define the effective $\gamma$, used in equation (18), as

$$ \gamma_{\text{eff}} = \frac{2Q\gamma}{\gamma Q + \frac{1}{2}\sqrt{2(\gamma^2 Q^2 + 1)^2 - 16Q^2(\gamma - 1) + 2\gamma^2 Q^2 - 2}}. \qquad (22) $$
The functions F, G and K in equation (11) are given by

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

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

and

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

Note that the p takes the form of $p_v$ (equation 16) or $p_\chi$ (equation 18) as defined above.

Following Papaloizou & Larwood (2000) and Cresswell & Nelson (2008), we define the planet's migration time-scale as

$$t_{\text{in}} = -\frac{L}{\Gamma_{\text{tot}}}. \qquad (26)$$

With this definition of $t_m$, the respective time-scale for a planet on circular orbit to reach the star is $t_m/2$. In equation (26), the quantity L is a planet's orbital angular momentum and $\Gamma_{tot}$ is the type-I torque defined in equation (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 time-scales are given by $t_e$ and $t_i$, respectively. They are defined as

$$t_e = \frac{t_{\text{wave}}}{0.780} \left( 1 - 0.14 \left(\frac{e}{h}\right)^2 + 0.06 \left(\frac{e}{h}\right)^3 + 0.18 \left(\frac{e}{h}\right) \left(\frac{i}{h}\right)^2 \right), \qquad (27)$$

and

$$t_i = \frac{t_{\text{wave}}}{0.544} \left( 1 - 0.3 \left(\frac{i}{h}\right)^2 + 0.24 \left(\frac{i}{h}\right)^3 + 0.14 \left(\frac{e}{h}\right)^2 \left(\frac{i}{h}\right)^2 \right), \qquad (28)$$

where

$$t_{\text{wave}} = \left(\frac{M_{\odot}}{m_p}\right) \left(\frac{M_{\odot}}{\Sigma_{\text{gas}} a^2}\right) h^4 \Omega_k^{-1}, \qquad (29)$$

and $M_{\odot}$, $a_p$, $m_p$, i and e are the solar mass and the embryo semimajor axis, mass, orbital inclination and eccentricity, respectively.

Using the previously defined time-scales, 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

$$\mathbf{a}_m = -\frac{\mathbf{v}}{t_m}, \qquad (30)$$

$$\mathbf{a}_e = -2 \frac{(\mathbf{v} \cdot \mathbf{r}) \mathbf{r}}{r^2 t_e}, \qquad (31)$$

and

$$\mathbf{a}_i = -\frac{v_z}{t_i} \mathbf{k}, \qquad (32)$$

where k is the unit vector in the z-direction. Equations (30)–(32) are given in Papaloizou & Larwood (2000) and Cresswell & Nelson (2006, 2008).

<figure><img src="image_2.png" /><figcaption>**Figure 1.** Evolution of the migration map calculated in a disc with metallicity equal to 1 per cent and $\alpha = 0.0054$ (dimensionless viscosity). The grey lines in each panel correspond to zero-torque locations and they delimit outward migration regions. At $t_{disc} = 0$ yr two different regions where outward migration is possible are shown. The most prominent one extends from about 5 to 20 au for planets with masses from 10 to more than 40 $M_\oplus$. As the disc evolves these regions move inward, shrink and eventually merge. The yellow region at about 0.1 au corresponds to the planet trap set at the disc inner edge.</figcaption></figure>

## 2.3 A migration map

Combining our disc 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 semimajor axis and mass.

Fig. 1 presents an evolving migration map of our chosen disc. It shows the direction and relative speed of migration of planets on circular and coplanar orbits at different locations within the disc. The direction of migration is represented by the colour. A negative torque (reddish to black) implies inward migration while a positive one (orange to light yellow) represents outward 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 disc inner edge (Masset et al. 2006).
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Evolution of a characteristic simulation of our fiducial set during the 5 Myr gas disc phase. Each blue circle corresponds to one embryo and the size of the point scales as $M_p^{1/3}$, where $M_p$ is the planet mass. The vertical axis shows the mass and the horizontal one shows the planet's semimajor axis. The grey dashed line delimits regions of outward migration and the disc inner edge. See also animated figure online.</figcaption></figure>

## 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.

In equation (33), $a_i$ and $a_j$ are the semimajor 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 that planetary embryos start to evolve in the disc corresponds to $t_{disc}=0$ yr.

### 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 disc model (Bitsch et al. 2015). The embryos’ initial masses are randomly and uniformly selected in the range from 0.1 to 4.5 $M_\oplus$. The total mass in planetary embryos in each simulation is about 60 $M_\oplus$. Adjacent planetary embryos were initially spaced by ~5 mutual hill radii $R_{H,m}$ (Kokubo & Ida 2000), where

$$R_{H,m} = \frac{a_i + a_j}{2} \left( \frac{M_i + M_j}{3M_\odot} \right)^{1/3}. \qquad (33)$$

### 3.2 Fiducial model: dynamical evolution

Figs 2 and 3 show the dynamical evolution of two characteristic simulations during the 5 Myr gas phase. The grey lines in each panel delimit the outward migration regions shown in Fig. 1. The imposed planet trap at the disc inner edge at ~0.1 au is also evident. Embryos migrate inward and converge to form resonant chains. Resonant planet pairs often become unstable and collide. As planetary embryos collide and grow (and the disc 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 inward, 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
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Another example of the dynamical evolution of the planetary embryos during the 5 Myr gas disc phase. The size of each blue circle scales as $M_p^{1/3}$, where $M_p$ is the planet mass. The vertical axis shows the mass and the horizontal one shows the planet's semimajor axis. The grey dashed line delimits regions of outward migration and the disc inner edge. See also animated figure online.</figcaption></figure>

planets larger than 5 $M_⊕$ migrate inward to the inner regions of the disc. This deposits planets in long resonant chains at the inner edge of the disc. Small planetary embryos – typically smaller than 1 $M_⊕$ migrate very slowly and stay beyond 1 au. At the end of the gas disc phase planetary systems exhibits compact resonant configurations with 5–10 planets inside 0.5 au.

Figs 4 and 5 show the dynamical evolution during and after the gas disc phase for the simulations shown in Figs 2 and 3. Recall that the gas lasts 5.1 Myr and the entire simulations last 100 Myr. The simulations from Figs 4 and 5 behaved in a similar fashion during the gas disc phase, but their later evolutionary paths are contrasting examples. The system from Fig. 4 remained stable after the gas disc dissipated, surviving in a long resonant chain containing six 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 Fig. 4 are less than 0.05 and their orbital inclinations are smaller than $0.1°$. The final eccentricities of planets in Fig. 5 are as high as 0.1 and their inclinations are as high as $6°$. The most massive planets were somewhat larger in the unstable simulation as

well (~18 $M_⊕$ for the simulation from Fig. 5 versus 11 $M_⊕$ for the simulation from Fig. 4).

## 4 TURBULENT DISCS

Hydrodynamical instabilities in gaseous protoplanetary disc generate turbulence in the disc 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 disc mid-plane (Johansen, Henning & Klahr 2006) and the vertical shear instability (Nelson, Gressel & Umurhan 2013; Stoll & Kley 2014). Another potentially important source for the observed gas accretion rate on young stars is turbulence driven by the magnetorotational instability (MRI; Balbus & Hawley 1998). In a sufficiently ionized and magnetized disc, the MRI generates magnetohydrodynamic turbulence that leads to outward angular momentum transport (Brandenburg et al. 1995; Armitage 1998). MRI turbulence produces large-scale, axisymmetric and long-lived density and pressure perturbations in the disc (Hawley, Gammie & Balbus 1996). Here, we perform simulations
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Dynamical evolution of planets in one simulation during and after the gas disc dispersal. The panels show the temporal evolution of planets' semimajor axis, eccentricity, mass and orbital inclination. The same line colour is used to consistently represent an individual planet in all panels. The gas dissipates at 5.1 Myr and the system is numerically integrated up to 100 Myr. This planetary system is dynamically stable after the gas disc phase for at least 100 Myr. The grey dashed vertical line shows the time of the disc dissipation.</figcaption></figure>

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, Steinacker & Adams 2004; Baruteau & Lin 2010). However, one should also note that recent studies have proposed that a number of non-ideal effects can suppress magnetorotational turbulence in a large region of the disc (Turner et al. 2014).

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

To model the turbulent motion of the gas, which essentially corresponds to density fluctuations in the gaseous disc, 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

$$F = -\Gamma \nabla \Phi, \qquad (34)$$

where

$$\Gamma = \frac{10^3 \sum_{\text{gas}} r^2}{\pi^2 M_{\star}}. \qquad (35)$$

Note that our $\Gamma$ is larger than that in Laughlin et al. (2004) by a factor of ~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

$$\Phi = \gamma r^2 \Omega^2 \sum_{i=1}^{N} \Lambda_{c,m}, \qquad (36)$$

where a single oscillation mode is defined as

$$\Lambda_{c,m} = \xi e^{-\frac{(r-r_c)^2}{\sigma^2}} \cos(m\theta - \phi_c + \Omega_c t) \sin\left(\frac{\tilde{\tau}}{\Delta t}\right). \qquad (37)$$
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Dynamical evolution of planets in one simulation during and after the gas disc dispersal. The panels show the temporal evolution of planets' semimajor axis, eccentricity, mass and orbital inclination. The same line colour is used to consistently represent an individual planet in all panels. The gas dissipates at 5.1 Myr and the system is numerically integrated up to 100 Myr. This planet system presents a phase of dynamical instability after the gas disc phase which lead to collisions and, consequently, to a planetary system dynamically less compact but relatively excited. The grey vertical line shows the time of the disc dissipation.</figcaption></figure>

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

$$ \gamma = 0.085h\sqrt{\alpha}. \qquad (38) $$

Also, as noted in equation (36), a specific mode is determined by the azimuthal wavenumber *m*, which we randomly sort with a logarithmic distribution between 1 and 96, the centre of its initial radial location *r*<sub>c</sub>, and azimuthal phase *φ*<sub>c</sub>. Only modes with wavenumber *m* smaller than 6 are considered (Ogihara et al. 2007). To sample *r*<sub>c</sub>, we use a lognormal distribution to select values between *r*<sub>in</sub> and *r*<sub>out</sub>. In our simulations *r*<sub>in</sub> = 0.1 au and *r*<sub>out</sub> = 25 au. The outer edge of the turbulent region (*r*<sub>out</sub>) 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 disc than most of previous studies. We set *N* = 125, similar to Horn et al. (2012). *φ*<sub>c</sub> is sorted with a uniform distribution between 0 and 2π. The dimensionless parameter *ξ* is sorted according to a Gaussian distribution with a unitary standard deviation and mean-value zero. Each mode has a radial extent of *σ* = πr<sub>c</sub>/4m.

The planet coordinates are represented by the radial distance *r* and its azimuthal coordinate *θ*. *Ω*<sub>c</sub> is the Keplerian frequency calculated at *r*<sub>c</sub>. Still, turbulent fluctuations appear and disappear in the disc. To account for this phenomenon in equation (35), the wave mode lifetime is defined by Δ*t* = 0.2πr<sub>c</sub>/mc<sub>s</sub> (Baruteau & Lin 2010), where c<sub>s</sub> is the local sound speed. Thus, a given mode *m* may only exist from its birth time *t* = *t*<sub>0</sub> to *t̃* = *t* − *t*<sub>0</sub> = Δ*t*. If *t̃* > Δ*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

$$ F_{\text{turb},r} = \gamma \Gamma r \Omega^2 \sum_{i=1}^{N} \left( 1 + \frac{2r(r - r_c)}{\sigma^2} \right) \Lambda_{c,m}, \qquad (39) $$

$$ F_{\text{turb},\theta} = \gamma \Gamma r \Omega^2 \sum_{i=1}^{N} \Lambda_{c,m}, \qquad (40) $$

and

$$ F_{\text{turb},z} = 0. \qquad (41) $$
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Dynamical evolution of planets in one simulation including the turbulent effects during and after the gas disc dispersal. The panels show the temporal evolution of planets' semimajor axis, eccentricity, mass and orbital inclination. The same line colour is used to consistently represent each planet in all panels. The gas dissipates at 5.1 Myr and the system is numerically integrated up to 100 Myr. This planetary system is dynamically stable after the gas disc phase for at least 100 Myr. The grey vertical line shows the time of the disc dissipation.</figcaption></figure>

For simplicity, in our model we assume that there is no feedback from the stochastic density/pressure fluctuations in the disc on our disc structure model. Basically, we assume the same underlying disc 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 Fig. 4. The result of this simulation is shown in Fig. 6.

Fig. 6 shows that turbulence is important in the outermost parts of the disc, visible by the 'random walks' in semimajor axis of planets beyond ~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 weaken. When planets reach the inner edge of the disc 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 per cent for the inner regions of the gaseous disc and a given $m$, the mode lifetime scales as $\Delta t \sim r_c^{3/2}$. Thus, modes in the inner regions of the disc (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 disc, 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 disc dissipates, its effects become 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.
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Cumulative distribution of the time of the last collision in our simulations. Cumulative distributions are separately calculated for the gas disc phase and after gas disc dissipation phase. Left: computed from collisions happening between 0 and 5 Myr. Right: computed from collisions happening from 5 to 100 Myr.</figcaption></figure>

## 5 SIMULATION OUTCOMES

Our simulations follow a bifurcated evolutionary path. During the gas disc phase, planetary embryos grow and migrate inward to the inner edge of the disc. 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 disc with a reduced lifetime would not help in avoiding these configurations. Resonant chains are typically established in less than 1.5 Myr (e.g. Figs 4–6). After the gas disc 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 analyse more in detail the evolution and outcome of simulations. We present both the fiducial and turbulent sets of simulations.

Fig. 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 disc phase, most systems have their last collisions after 1 Myr (left-hand panel). This is expected since this corresponds to when most planets are approaching the disc inner edge (see for example Figs 5 and 6), and have already reached relatively more compact configurations. This generates dynamical instabilities and collisions. Fig. 7 also shows that during the gas disc phase the cumulative distributions of our fiducial and turbulent models are similar. The cumulative distributions of the last collision epoch grow 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 to take place and every forming planetary system exhibits at least a few collisions during the gas disc phase. Because the disc 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-hand 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 starts to drop. Also, the cumulative distributions do not reach 1. This indicates that not all planetary systems underwent instabilities

after the gas disc phase (or did not experience any collision). About 60 per cent of the fiducial simulations and 50 per cent of the turbulent simulations experienced dynamical instabilities and at least one collision after the gas disc phase. Comparatively, only 20 per cent 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 time-scales because of the very small time-step necessary to resolve the orbits of planets that reach the inner orbit of the disc and the consequent long CPU time demanded. However, we do not expect a linear growth with time of the number of unstable systems.

Fig. 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 disc (Fig. 8, left-hand 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 disc (Izidoro, Morbidelli & Raymond 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 (Fig. 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 disc 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 stable 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 Fig. 9.

We first focus on the results of our fiducial model. Observing the period ratio distribution of adjacent planet pairs in Fig. 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). Fig. 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. Fig. 9 also show how the masses of planets in stable and unstable systems compare to each other at 5 Myr. The mass cumulative distributions of Fig. 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.
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Representative sample of planetary systems produced in our fiducial and turbulent models at 5 and 100 Myr. The vertical axis shows the mass and the horizontal one represents the planets' semimajor axis. **Left:** planetary systems at 5 Myr. **Middle:** planetary systems at 100 Myr. For each model, we randomly selected 20 planetary systems that are plotted together. **Right:** selected observed planetary systems. The line connects planets in the same system.</figcaption></figure>

The results of Fig. 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 Fig. 9 are qualitatively similar in some key aspects. The most notable one is that, before gas dispersal, planet pairs are essentially found in compact resonant configurations. However, this is about to change for unstable systems. While stable systems configuration 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 5 Myr) 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.

Fig. 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 two to five 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 per cent more fiducial simulations were unstable at later times (Fig. 7). This difference probably comes from a subtle difference in some characteristics of the resonant chains (e.g. the number of planets in the chain, commensurabilities and masses of planet pairs (Matsumoto, Nagasawa & Ida 2012), amplitude of libration of resonant angles (Adams et al. 2008), etc. Perhaps part of this difference is also due to the small number of statistics.

Fig. 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 d and radii smaller than 4 R⊕. The motivation for choosing these cut-offs comes from the completeness of Kepler observations (Silburt, Gaidos & Wu 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 disc 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 Fig. 9). There is little difference between the fiducial and turbulent simulations. Instabilities break the resonant chains created during the gas disc phase, promote collisions and scattering events that reduce the number of planets in the system (Fig. 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 Figs 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 3.<fnref n="1" /> This suggests 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.)

Fig. 12 shows the cumulative distribution of semimajor axis, mass, orbital eccentricity and mutual inclination of simulated planets inside 0.5 au. Cumulative distributions of semimajor axis (left-hand upper panel of Fig. 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

<footnote n="1">¹ The tail of the period ratio distribution probably suffers from selection bias effects.</footnote>
<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 9.** Cumulative distributions of period ratio of adjacent planets (top), the number of planets in the resonant chain (middle) and masses of planets (bottom) at 5 Myr. Only planets inside 0.5 au are considered. The grey line in the top panel shows the period ratio distribution of adjacent planet pairs in Kepler systems.</figcaption></figure>

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 pile-up of planets at about 0.15 au, near the disc 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.

<figure><img src="image_13.png" /><figcaption>**Figure 10.** Cumulative distribution of the total number of planets inside 0.5 au produced in our fiducial and turbulent model. The cumulative distributions are shown at 5 and 100 Myr. Only unstable planetary systems are considered here.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 11.** Cumulative period ratio distributions of adjacent planets produced in our fiducial and turbulent model for dynamically unstable systems, before and after the dynamical instability. The grey line corresponds to the observed period ratio of adjacent planets in the Kepler data.</figcaption></figure>

The median planet mass after instability is about three times larger. This is simply due to the fact that the unstable systems underwent a late phase of collisions<fnref n="2" />. 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; Figueira et al. 2012; Kane et al. 2012; Tremaine & Dong 2012; Fabrycky et al. 2014; Plavchan, Bilinski & Currie 2014; Van Eylen & Albrecht 2015; Ballard & Johnson 2016). To represent the eccentricity distribution of observations derived from statistical studies we used a Rayleigh distribution with $\sigma_e = 0.1$ (e.g. Moorhead et al. 2011). For the inclination distribution we used a Rayleigh distribution with $\sigma_i = 1.5$ (e.g. Fang & Margot 2012).

<footnote n="2">² Note that, in the inner regions, collisions are much more common than ejections. It is an issue of Safronov number (Safronov 1972).</footnote>
<figure><img src="image_15.png" /><figcaption>**Figure 12.** Cumulative distribution of semimajor axis (top-left), mass (top-right), orbital eccentricity (bottom-left) and mutual inclination (bottom-right) of planets inside 0.5 au in our simulations. The eccentricity and orbital inclination distributions of observed Kepler planets derived from statistical analysis are shown for comparison in the respective lower panels. The eccentricity distribution shown by the grey line follows a Rayleigh distribution with $\sigma_e = 0.1$ (Moorhead et al. 2011), while the mutual inclination distribution (grey line) follows a Rayleigh distribution with $\sigma_i = 1.5$ (Fang & Margot 2012). The latter put more than 85 per cent of the planets with orbital inclination smaller than 3°.</figcaption></figure>

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 compares 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 disc 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

Fig. 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 Fig. 11 (see Fig. 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 modestly more spread out than observed Kepler systems.
<figure><img src="image_16.png" /><figcaption>**Figure 13.** The separation of planet pairs in unstable simulations as a function of the total mass. The low-mass planet pairs have $M_{\text{tot}} = M_1 + M_2 < 25 \, \text{M}_{\oplus}$, the high-mass have $M_{\text{tot}} > 34.5 \, \text{M}_{\oplus}$, and medium-mass pairs are in between. Left: period ratio distribution. Right: separation in units of mutual Hill radii. Masses for the Kepler systems were ributed using the probabilistic mass-radius relation of Wolfgang et al. (2016).</figcaption></figure>

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_{\text{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_{\text{tot}} < 25 \, \text{M}_{\oplus}$, the high-mass pairs have $M_{\text{tot}} > 34.5 \, \text{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.

Fig. 13 (right-hand 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-hand 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, Rogers & Ford (2016): $M = 2.7(R/R_{\oplus})^{1.3}$. This fits nicely with the results of Fang & Margot (2012), who inferred that Kepler systems are typically spaced by roughly 20 ± 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 indicates that they are signposts of stable systems. Matching observations thus requires that a fraction of resonant chains remain stable after dissipation of the gas disc.

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 cut-off of $4 \, \text{R}_{\oplus}$ in the Kepler data, the maximum mass inferred using the Wolfgang mass-radius relationship is $16.4 \, \text{M}_{\oplus}$ and the median planet mass in the sample is $6.3 \, \text{M}_{\oplus}$. The median mass of unstable planets at the end of our simulations is $14 \, \text{M}_{\oplus}$.$^3$

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 2011; 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.

Fig. 14 (left-hand 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 Figs 4 and 5) are more strongly affected by observational bias. Naturally, the orbital configuration of the

<footnote>$^3$It is interesting to note that the median mass of the planets trapped in resonant chains is $5.6 \, \text{M}_{\oplus}$. This is comparable to the inferred value for Kepler sample of $6.3 \, \text{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.</footnote>
<figure><img src="image_17.png" /><figcaption>**Figure 14.** The outcome of our simulated observations. Left-hand panel: the period ratio distribution of the observed Kepler systems (thick grey line), as well as our stable (solid red line) and unstable (solid blue line) simulations. The dashed lines show the planet pairs retrieved by simulated observations. Right-hand panel: *P* values from K-S tests comparing the Kepler systems with a sample of simulated planet pairs with *M* < 16.4 M<sub>⊕</sub>, after taking observational bias into account with simulated observations. The dashed/dotted line is at *p* =10 per cent/1 per cent.</figcaption></figure>

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-fitting 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 cut-off of 16.4 M<sub>⊕</sub>, which corresponds to the mass of a 4 R<sub>⊕</sub> planet with the Wolfgang et al. (2016) mass-radius relation. This is our upper size cut-off for Kepler planets. We also restricted this analysis to systems with P<sub>2</sub>/P<sub>1</sub> ≤ 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.

Fig. 14 (right-hand panel) shows that our simulated observations roughly match the Kepler sample if less than 25 per cent of planet pairs come from stable simulations (for a probability *p* ≥ 10 per cent 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 per cent (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 correlate 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 \text{ M}_\oplus, P_2/P_1 \le 3$), the best match to observations was for the smallest contribution from stable systems. However, up to nearly 20 per cent of stable systems were allowed while maintaining $p \ge 10$ per cent.

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 ~25 per cent 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 2016). 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° above the initial *i* = 0 plane to 30° 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, Fig. 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 per cent of detections were of a single planet in transit, and 66 per cent of detections had *N* ≥ 3. In contrast, for unstable systems 78 per cent of outcomes were single-planet detections and only 7 per cent of detections had three or more planets.
<figure><img src="image_18.png" /><figcaption>**Figure 15.** The number of planets detected in synthetic observations of our simulated planetary systems. The blue/red curves represent the unstable/stable simulations (combining the fiducial and turbulent sets). The grey curve is a 90–10 mixture of the unstable and stable simulated systems, respectively. The thick green curve is the Kepler sample, removing single giant planet systems but keeping systems with giants and super-Earths.</figcaption></figure>

If we combine stable and unstable systems in a 1-to-9 ratio (i.e. with 90 per cent unstable and 10 per cent 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 \ge 4$) are more often stable. The Kepler-223 (Mills et al. 2016) and TRAPPIST-1 (Gillon et al. 2017) multiresonant super-Earth systems appear to be good examples of high-$N$ stable systems.

This would suggest that the Kepler dichotomy is simply an observational artefact. 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 Fig. 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 varied 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 (2016) 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^{\circ}$ to $10^{\circ}$. The results of this experiment matched qualitatively those of Johansen et al. (2012). Fig. 16 (left-hand 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 (Fig. 16, right-hand 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 is 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 disc 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 discs 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 discs may vary from disc to disc but *in situ* accretion simulations also predict that the radial distribution of planets in multiplanet systems should be mass ranked. More specifically, they predict that planet mass should decrease with semimajor 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
<figure><img src="image_19.png" /><figcaption>**Figure 16.** Comparison of the distribution of the number of planets in our simulations and statistical studies. Left: distribution of the number of planets of the observed Kepler systems (thick grey line), as well as our unstable simulations (dashed line), and the distribution derived from different statistical distributions of mutual inclinations of planets. Right: inclination distributions in our stable systems (blue line), as well as in all our unstable systems (red line), the mix of unstable of unstable systems (90 per cent unstables plus 10 per cent stables; green line), Rayleigh distribution with $\sigma = 1.5$ (dashed line; see Fang & Margot 2012) and Rayleigh distribution with $\sigma = 4^\circ$ (dash-dotted line).</figcaption></figure>

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 \ge 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).

Finally, it is also important to note that Kepler planets in single and multiple planet systems have statistically different orbital eccentricities distributions. Single Kepler planets are prone to have much higher orbital eccentricities than planets in multiple planet systems (Xie et al., 2016). Our results are consistent with this finding.

## 7 COMMENTS ON DIFFERENT DISC 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 disc 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 disc, both at the end of the disc 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 disc 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 us 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 disc 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 disc 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 provide only 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 catalogue. Planets' eccentricities and orbital inclinations were set to zero. To mimic type-I migration and eccentricity damping, Rein (2012) used simple migration time-scales of about $10^3-10^4$ yr. 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 us first focus on the gas disc phase.

The first great difference between our model and Rein's is how synthetic forces were implemented. In the Rein (2012) 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 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
catalogue (see fig. 1 in Rein 2012). Basically, the approach in Rein (2012) strongly favours the production of dynamically relaxed planetary systems. While Rein's simulations had short integration time, our protoplanetary discs lived for about 5 Myr. Thus, planets in our simulations had the time to migrate to the inner edge of the disc and pile up in long chain of resonances before the gas dissipated (e.g. Fig. 6). In our simulations, the disc surface density decreases with increasing time according to the disc model derived from hydrodynamical simulations (Bitsch et al. 2015). Finally, the third reason that 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 disc 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 disc surface density, aspect ratio and the distance between the centre 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 equation 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 \ge 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 discs.

It is also interesting to note that our results seem to agree well with the conclusions from a very recent paper by Batygin & Adams (2017), where they derived an analytic criterion for turbulent disruption of mean-motion resonances. According to these authors – at the inner regions of the disc – 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 per cent (50 per cent) of our planet pairs have combined masses larger than $6 M_{\oplus} (10 M_{\oplus})$, this fits nicely with our results.

## 8 WHY DO SO MANY RESONANT CHAINS GO UNSTABLE?

Fig. 14 shows that our simulations can match the sample of Kepler planets if at least 75 per cent of resonant chains go unstable and at most 25 per cent 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$ per cent 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$ per cent of resonant chains remain stable. This yields an acceptable match to the Kepler period ratio distribution (Fig. 14).

Only 50–60 per cent of resonant chains were unstable in our fiducial and turbulent sets of simulations (Fig. 7). This is a far cry from the roughly 95 per cent 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 disc, 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 disc remains fixed at $\sim 0.1$ au during all disc lifetime. However, in reality, the inner edge of the disc should move as gas the disc and star evolve, due to the balance between the stellar magnetic pressure and the viscous torque of the disc (e.g. Koenigl 1991). Liu, Ormel & Lin (2017) showed that under certain conditions two planets can experience divergent migration as the disc inner edge moves outward. This mechanism could be also another important inductor of dynamical instabilities.

We have considered neither gas accretion on to planetary cores (e.g. Ginzburg, Schlichting & Sari 2016), nor mass-loss during planetary collisions (e.g. Inamdar & Schlichting 2016) and 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 disc. Leftover planetesimals in the disc 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 disc viscosity, the initial distribution of planetary embryos and photoevaporation in the disc. 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 one to two 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 disc model (Bitsch et al.
2015). We found that stochastic forcing from disc 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 disc and migrate inward due to torques from the gaseous disc (Figs 1–3). When the first embryo reaches the inner edge of the disc, its migration is stopped by the planet disc-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 disc dissipates, about 50–60 per cent 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-hand panel). Our simulations match the period ratio distribution of the observed Kepler planet pairs if 75–100 per cent 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 per cent. 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. AI, SNR, AM, AP and FH 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. AI thanks also partial financial support from FAPESP (Proc. 16/19556-7 and 16/12686-2) during the final preparation of this paper. BB thanks the Knut and Alice Wallenberg Foundation for their financial support.

## REFERENCES

Adams F., Laughlin G., Bloch A., 2008, ApJ, 683, 1117

Armitage P. J., 1998, ApJ, 501, L189

Artymowicz P., 1993, ApJ, 419, 155

Balbus S. A., Hawley J. F., 1998, Rev. Mod. Phys., 70, 1

Ballard S., Johnson J. A., 2016, ApJ, 816, 66

Baruteau C., Lin D. N. C., 2010, ApJ, 709, 759

Baruteau C., Masset F., 2008, ApJ, 672, 1054

Baruteau C., Papaloizou J. C. B., 2013, ApJ, 778, 7

Baruteau C. et al., 2014, in Beuther H., Klessen R. S., Dullemond C. P., Henning T., Protostars Planets VI. Univ. Arizona Press, Tucson, AZ, p. 4293

Batalha N. M. et al., 2013, ApJS, 204, 24

Batygin K., 2015, MNRAS, 451, 2589

Batygin K., Adams F. C., 2017, ApJ, 153, 120

Batygin K., Morbidelli A., 2013, AJ, 145, 1

Bell K. R., Lin D. N. C., 1994, ApJ, 427, 987

Bitsch B. et al., 2013, A&A, 549, A124

Bitsch B., Kley W., 2010, A&A, 523, A30

Bitsch B., Kley W., 2011, A&A, 530, A41

Bitsch B., Morbidelli A., Lega E., Crida A., 2014, A&A, 564, A135

Bitsch B., Johansen A., Lambrechts M., Morbidelli A., 2015, A&A, 575, A28

Boley A. C., Ford E. B., 2013, preprint (arXiv:1306.0566)

Boley A. C., Morris M. A., Ford E. B., 2014, ApJ, 792, L27

Borucki W. J. et al., 2010, Science, 327, 977

Borucki W. J. et al., 2011, ApJ, 736, 19

Bouvier J., 2013, in Hennebelle P., Charbonnel C., eds, EAS Publ. Ser. Vol. 62, Role and Mechanisms of Angular Momentum Transport During the Formation and Early Evolution of Stars, Evry Schatzman School 2012. p. 143

Bouvier J., Alencar S. H. P., Harries T. J., Johns-Krull C. M., Romanova M. M., 2007, in Reipurth B., Jewitt D., Keil K., eds, Protostars and Planets V. Univ. Arizona Press, Tucson, AZ, p. 479

Brandenburg A., Nordlund A., Stein R. F., Torkelsson U., 1995, ApJ, 446, 741

Chambers J. E., 1999, MNRAS, 304, 793

Charbonneau D., Brown T. M., Burrows A., Laughlin G., 2007, in Reipurth B., Jewitt D., Keil K., eds, Protostars and Planets V. Univ. Arizona Press, Tucson, AZ, p. 701

Chatterjee S., Ford E. B., 2015, ApJ, 803, 33

Chatterjee S., Tan J. C., 2014, ApJ, 780, 53

Chatterjee S., Tan J. C., 2015, ApJ, 798, L32

Chiang E., Laughlin G., 2013, MNRAS, 431, 3444

Coleman G. A. L., Nelson R. P., 2014, MNRAS, 445, 479

Coleman G. A. L., Nelson R. P., 2016, MNRAS, 457, 2480

Cossou C., Raymond S. N., Pierens A., 2013, A&A, 553, L2

Cossou C., Raymond S. N., Hersant F., Pierens A., 2014, A&A, 569, A56

Cresswell Nelson R. P., 2008, A&A, 482, 677

Cresswell P., Nelson R. P., 2006, A&A, 450, 833

Crida A., Morbidelli A., Masset F., 2006, Icarus, 181, 587

Delisle J.-B., Laskar J., 2014, A&A, 570, L7

Delisle J.-B., Laskar J., Correia A. C. M., Boué G., 2012, A&A, 546, A71

Delisle J.-B., Laskar J., Correia A. C. M., 2014, A&A, 566, A137

Dong S., Zhu Z., 2013, ApJ, 778, 53

Fabrycky D. C. et al., 2012, ApJ, 750, 114

Fabrycky D. C. et al., 2014, ApJ, 790, 146

Fang J., Margot J.-L., 2012, ApJ, 761, 92

Fendyke S. M., Nelson R. P., 2014, MNRAS, 437, 96

Figueira P. et al., 2012, A&A, 541, A139

Flock M., Fromang S., Turner N. J., Benisty M., 2017, ApJ, 835, 230

Fressin F. et al., 2013, ApJ, 766, 81

Fromang S., Papaloizou J., 2007, A&A, 476, 1113

Gillon M. et al., 2017, Nature, 542, 456

Ginzburg S., Schlichting H. E., Sari R., 2016, ApJ, 825, 29

Goldreich P., Sari R., 2003, ApJ, 585, 1024

Goldreich P., Schlichting H. E., 2014, AJ, 147, 32

Goldreich P., Tremaine S., 1979, ApJ, 233, 857

Goldreich P., Tremaine S., 1980, ApJ, 241, 425

Haghighipour N., 2013, Annu. Rev. Earth Planet. Sci., 41, 469

Hansen B. M. S., 2014, MNRAS, 440, 3545

Hansen B. M. S., Murray N., 2012, ApJ, 751, 158

Hansen B. M. S., Murray N., 2013, ApJ, 775, 53

Hartmann L., Calvet N., Gullbring E., D'Alessio P., 1998, ApJ, 495, 385

Hawley J. F., Gammie C. F., Balbus S. A., 1996, ApJ, 464, 690

Hellary P., Nelson R. P., 2012, MNRAS, 419, 2737

Horn B., Lyra W., Mac Low M.-M., Sándor Z., 2012, ApJ, 750, 34

Howard A. W., 2013, Science, 340, 572

Howard A. W. et al., 2012, ApJS, 201, 15

Hu X., Zhu Z., Tan J. C., Chatterjee S., 2016, ApJ, 816, 19

Ida S., Lin D. N. C., 2008, ApJ, 685, 584

Ida S., Lin D. N. C., 2010, ApJ, 719, 810

Inamdar N. K., Schlichting H. E., 2016, ApJ, 817, L13

Izidoro A., Morbidelli A., Raymond S. N., 2014, ApJ, 794, 11

Izidoro A., Raymond S. N., Morbidelli A., Hersant F., Pierens A., 2015, ApJ, 800, L22

Johansen A., Henning T., Klahr H., 2006, ApJ, 643, 1219

Johansen A., Davies M. B., Church R. P., Holmelin V., 2012, ApJ, 758, 39
Kane S. R., Ciardi D. R., Gelino D. M., von Braun K., 2012, MNRAS, 425, 757

Ketchum J. A., Adams F. C., Bloch A. M., 2011, ApJ, 726, 53

Klahr H. H., Bodenheimer P., 2003, ApJ, 582, 869

Kley W., Crida A., 2008, A&A, 487, L9

Kley W., Nelson R. P., 2012, ARA&A, 50, 211

Kley W., Bitsch B., Klahr H., 2009, A&A, 506, 971

Koenigl A., 1991, ApJ, 370, L39

Kokubo E., Ida S., 2000, Icarus, 143, 15

Laughlin G., Steinacker A., Adams F. C., 2004, ApJ, 608, 489

Lecoanet D., Adams F. C., Bloch A. M., 2009, ApJ, 692, 659

Leinhardt Z. M., Stewart S. T., 2012, ApJ, 745, 79

Lin D. N. C., Papaloizou J., 1979, MNRAS, 186, 799

Lin D. N. C., Papaloizou J., 1986, ApJ, 309, 846

Lissauer J. J. et al., 2011a, ApJS, 197, 8

Lissauer J. J. et al., 2011b, Nature, 470, 53

Lithwick Y., Wu Y., 2012, ApJ, 756, L11

Liu B., Zhang X., Lin D. N. C., Aarseth S. J., 2015, ApJ, 798, 62

Liu B., Zhang X., Lin D. N. C., 2016, ApJ, 823, 162

Liu B., Ormel C. W., Lin D. N. C., 2017, A&A, 601, A15

Lovelace R. V. E., Li H., Colgate S. A., Nelson A. F., 1999, ApJ, 513, 805

Lyra W., Paardekooper S.-J., Mac Low M.-M., 2010, ApJ, 715, L68

MacDonald M. G. et al., 2016, AJ, 152, 105

Masset F. S., Morbidelli A., Crida A., Ferreira J., 2006, ApJ, 642, 478

Mathis S., 2015, A&A, 580, L3

Matsumoto Y., Nagasawa M., Ida S., 2012, Icarus, 221, 624

Mayor M. et al., 2011, preprint (arXiv:1109.2497)

McNeil D. S., Nelson R. P., 2010, MNRAS, 401, 1691

Mills S. M., Fabrycky D. C., Migaszewski C., Ford E. B., Petigura E., Isaacson H., 2016, Nature, 533, 509

Moorhead A. V. et al., 2011, ApJS, 197, 1

Morbidelli A., Raymond S. N., 2016, J. Geophys. Res., 121, 1962

Morbidelli A., Lambrechts M., Jacobson S., Bitsch B., 2015, Icarus, 258, 418

Morbidelli A. et al., 2016, Icarus, 267, 368

Moriarty J., Ballard S., 2016, ApJ, 832, 34

Mullally F. et al., 2015, ApJS, 217, 31

Nelson R. P., 2005, A&A, 443, 1067

Nelson R. P., Gressel O., 2010, MNRAS, 409, 639

Nelson R. P., Papaloizou J. C. B., 2003, MNRAS, 339, 993

Nelson R. P., Gressel O., Umurhan O. M., 2013, MNRAS, 435, 2610

Ogihara M., Ida S., 2009, ApJ, 699, 824

Ogihara M., Kobayashi H., 2013, ApJ, 775, 34

Ogihara M., Ida S., Morbidelli A., 2007, Icarus, 188, 522

Ogihara M., Morbidelli A., Guillot T., 2015a, A&A, 578, A36

Ogihara M., Morbidelli A., Guillot T., 2015b, A&A, 584, L1

Paardekooper S.-J., Mellema G., 2006, A&A, 459, L17

Paardekooper S.-J., Mellema G., 2008, A&A, 478, 245

Paardekooper S.-J., Papaloizou J. C. B., 2008, A&A, 485, 877

Paardekooper S.-J., Baruteau C., Crida A., Kley W., 2010, MNRAS, 401, 1950

Paardekooper S.-J., Baruteau C., Kley W., 2011, MNRAS, 410, 293

Papaloizou J. C. B., 2011, Celest. Mech. Dyn. Astron., 111, 83

Papaloizou J. C. B., Larwood J. D., 2000, MNRAS, 315, 823

Papaloizou J., Lin D. N. C., 1984, ApJ, 285, 818

Papaloizou J. C. B., Nelson R. P., 2003, MNRAS, 339, 983

Papaloizou J. C. B., Terquem C., 2010, MNRAS, 405, 573

Petigura E. A., Howard A. W., Marcy G. W., 2013, Proc. Natl. Acad. Sci., 110, 19273

Pierens A., Baruteau C., Hersant F., 2011, A&A, 531, A5

Pierens A., Baruteau C., Hersant F., 2012, MNRAS, 427, 1562

Pierens A., Cossou C., Raymond S. N., 2013, A&A, 558, A105

Plavchan P., Bilinski C., Currie T., 2014, PASP, 126, 34

Pu B., Wu Y., 2015, ApJ, 807, 44

Raymond S. N., Cossou C., 2014, MNRAS, 440, L11

Raymond S. N., Morbidelli A., 2014, in Knežević Z., Lemaitre A., eds, Proc. IAU Symp. 310, Complex Planetary Systems. Cambridge Univ. Press, Cambridge, p. 194

Raymond S. N., Barnes R., Mandell A. M., 2008, MNRAS, 384, 663

Rein H., 2012, MNRAS, 427, L21

Rein H., Papaloizou J. C. B., 2009, A&A, 497, 595

Rein H., Payne M. J., Veras D., Ford E. B., 2012, MNRAS, 426, 187

Rivera E. J., Laughlin G., Butler R. P., Vogt S. S., Haghighipour N., Meschiari S., 2010, ApJ, 719, 890

Romanova M. M., Lovelace R. V. E., 2006, ApJ, 645, L73

Romanova M. M., Ustyugova G. V., Koldoba A. V., Wick J. V., Lovelace R. V. E., 2003, ApJ, 595, 1009

Rowe J. F. et al., 2014, ApJ, 784, 45

Safronov V. S., 1972, Evolution of the Protoplanetary Cloud and Formation of the Earth and Planets. Keter Publishing House, Israel

Sari R., Goldreich P., 2004, ApJ, 606, L77

Schlaufman K. C., 2014, ApJ, 790, 91

Schlichting H. E., 2014, ApJ, 795, L15

Schneider J., Dedieu C., Le Sidaner P., Savalle R., Zolotukhin I., 2011, A&A, 532, A79

Shakura N. I., Sunyaev R. A., 1973, A&A, 24

Silburt A., Gaidos E., Wu Y., 2015, ApJ, 799, 180

Spalding C., Batygin K., 2016, ApJ, 830, 5

Stewart S. T., Leinhardt Z. M., 2012, ApJ, 751, 32

Stoll M. H. R., Kley W., 2014, A&A, 572, A77

Suzuki T. K., Ogihara M., Morbidelli A., Crida A., Guillot T., 2016, A&A, 596, A74

Tanaka H., Ward W. R., 2004, ApJ, 602, 388

Tanaka H., Takeuchi T., Ward W. R., 2002, ApJ, 565, 1257

Terquem C., Papaloizou J. C. B., 2007, ApJ, 654, 1110

Tremaine S., Dong S., 2012, AJ, 143, 94

Turner N. J., Fromang S., Gammie C., Klahr H., Lesur G., Wardle M., Bai X.-N., 2014, in Beuther H., Klessen R. S., Dullemond C. P., Henning T., eds, Protostars and Planets VI. Univ. Arizona Press, Tucson, AZ, p. 411

Van Eylen V., Albrecht S., 2015, ApJ, 808, 126

Ward W. R., 1986, Icarus, 67, 164

Ward W., 1997, Icarus, 126, 261

Weiss L. M., Marcy G. W., 2014, ApJ, 783, L6

Williams J. P., Cieza L. A., 2011, ARA&A, 49, 67

Winn J. N., 2011, in Seager S., ed., Exoplanet Transits and Occultations. Univ. Arizona Press, Tucson, AZ

Winters W. F., Balbus S. A., Hawley J. F., 2003, ApJ, 589, 543

Wolfgang A., Rogers L. A., Ford E. B., 2016, ApJ, 825, 19

Wright J. T., Gaudi B. S., 2013, in Oswalt T. D., French L. M., Kalas P., eds, Planets, Stars and Stellar Systems. Springer-Verlag, Berlin, p. 489

Xie J. W. et al., 2016, Proc. Natl. Acad. Sci., 113, 11431

## SUPPORTING INFORMATION

Supplementary data are available at MNRAS online.

**Figure 2.** Evolution of a characteristic simulation of our fiducial set during the 5 Myr gas disc phase.

**Figure 3.** Another example of the dynamical evolution of the planetary embryos during the 5 Myr gas disc phase.

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the article.

This paper has been typeset from a TEX/LATEX file prepared by the author.