Formation of planetary systems by pebble accretion and migration: Hot super-Earth systems from breaking compact resonant chains

TL;DR

This study uses N-body simulations with pebble accretion and migration models to explain the formation and dynamical evolution of hot super-Earth systems, matching observed period ratios and instability rates.

Contribution

It introduces a comprehensive model combining pebble accretion and migration to reproduce the formation, resonant chain breaking, and observed characteristics of super-Earth systems.

Findings

Majority of resonant chains become unstable after gas dispersal

Simulations match observed period ratio distributions

Most systems are inherently multiple with high mutual inclinations

Abstract

At least 30\% of main sequence stars host planets with sizes of between 1 and 4 Earth radii and orbital periods of less than 100 days. We use N-body simulations including a model for gas-assisted pebble accretion and disk--planet tidal interaction to study the formation of super-Earth systems. We show that the integrated pebble mass reservoir creates a bifurcation between hot super-Earths or hot-Neptunes ($\lesssim15M_{\oplus}$) and super-massive planetary cores potentially able to become gas giant planets ($\gtrsim15M_{\oplus}$). Simulations with moderate pebble fluxes grow multiple super-Earth-mass planets that migrate inwards and pile up at the inner edge of the disk forming long resonant chains. We follow the long-term dynamical evolution of these systems and use the period ratio distribution of observed planet-pairs to constrain our model. Up to $\sim$95\% of resonant chains become dynamically unstable after the gas disk dispersal, leading to a phase of late collisions that breaks the original resonant configurations. Our simulations naturally match observations when they produce a dominant fraction ($\gtrsim95\%$) of unstable systems with a sprinkling ($\lesssim5\%$) of stable resonant chains (the Trappist-1 system represents one such example). Our results demonstrate that super-Earth systems are inherently multiple (${\rm N\geq2}$) and that the observed excess of single-planet transits is a consequence of the mutual inclinations excited by the planet--planet instability. In simulations in which planetary seeds are initially distributed in the inner and outer disk, close-in super-Earths (abridged).

Tables

Table 1. Table 1: Initial conditions of our simulations. From left to right, columns mark the simulation model, the region where planetary seeds are initially distributed in the disk, their initial mutual radial separation, and the disk age at the start of the simulation. The mutual separation is also given in units of mutual hill radii, Rhm = ( m 1 + m 2 3 ​ M ⊕ ) 1 / 3 ​ a 1 + a 2 2 Rhm superscript subscript m 1 subscript m 2 3 subscript M direct-sum 1 3 subscript a 1 subscript a 2 2 {\rm Rhm=\left(\frac{m_{1}+m_{2}}{3M_{\oplus}}\right)^{1/3}\frac{a_{1}+a_{2}}{2}} , where m 1 subscript m 1 {\rm m_{1}} and m 2 subscript m 2 {\rm m_{2}} are the masses and a 1 subscript a 1 {\rm a_{1}} and a 2 subscript a 2 {\rm a_{2}} are the semi-major axes of any two adjacent planetary embryos.

		Radial	Mutual	${\text{t}}_{\text{start}}$
		Distribution	separation
Model	I	0.7 to $\sim 20\mathrm{AU}$	0.25 AU ($\sim 5-100\mathrm{Rhm}$)	0.5 and 3 Myr
	II	0.7 to $\sim 60\mathrm{AU}$	$\sim 0.01-3\mathrm{AU}$ ($\sim 10-30\mathrm{Rhm}$)	0.5 and 3 Myr
	III	$0.2\left(\mathrm{or}0.3\right)$ to $\sim 2\mathrm{AU}$	0.025 AU ($\sim 15-55\mathrm{Rhm}$)	0.5 and 3 Myr

Table 2. Table 2: Selected setups to model the formation and long-term dynamical evolution of close-in super-Earths based on the results of Section 3 . From left to right, the columns are name of the model, disk age at the start of the simulation ( t start subscript 𝑡 start {t_{{\rm start}}} ), the scaled pebble flux ( S peb subscript 𝑆 peb {S_{{\rm peb}}} ), and the integrated pebble flux ( I peb subscript I peb {\rm I_{{\rm peb}}} ). In our nominal simulations, Rocky peb = 0.1 ​ cm subscript Rocky peb 0.1 cm {\rm Rocky_{{\rm peb}}=0.1\leavevmode\nobreak\ cm} but we also performed simulations with Rocky peb = 1 ​ cm subscript Rocky peb 1 cm {\rm Rocky_{{\rm peb}}=1\leavevmode\nobreak\ cm} .

Scenarios	$t_{\mathrm{start}}$	$S_{\mathrm{peb}}$	$I_{\mathrm{peb}}\left(M_{\oplus }\right)$
Model-I	0.5 Myr	1	$\sim 194$
Model-I	3.0 Myr	5	$\sim 150$
Model-II	3.0 Myr	5	$\sim 150$
Model-III	0.5 Myr	10	$\sim 1940$

Equations

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

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

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

$$\Re =\tanh \left(\frac{r-r_{\mathrm{in}}}{0.05r_{\mathrm{in}}}\right),$$

$${\dot{M}}_{\mathrm{core}}=\pi R_{\mathrm{acc}}^2{\rho }_{\mathrm{p},\mathrm{mid}}\bar{S}\delta v,$$

$${\rho }_{\mathrm{p},\mathrm{mid}}={\Sigma }_{\mathrm{peb}}/(\sqrt{2\pi }{H}_{\mathrm{peb}}),$$

$$\bar{S}=\frac{1}{\pi R_{\mathrm{acc}}^2}{\int }_{z_0-R_{\mathrm{acc}}}^{z_0+R_{\mathrm{acc}}}2\exp [-z^2/(2H_{\mathrm{peb}}^2)]\sqrt{R_{\mathrm{acc}}^2-(z-z_0{)}^2}\mathrm{dz},$$

$${\dot{M}}_{\mathrm{peb}}=2\pi r_{\mathrm{g}}\frac{\mathrm{d}{r}_{\mathrm{g}}}{\mathrm{d}t}Z{\Sigma }_{\mathrm{g}}(r_{\mathrm{g}}).$$

$$r_{\mathrm{g}}={\left(\frac{3}{16}\right)}^{1/3}(G{M}_{\star }{)}^{1/3}({ϵ}_{\mathrm{D}}Z{)}^{2/3}{t}^{2/3},$$

$$\frac{\mathrm{d}{r}_{\mathrm{g}}}{\mathrm{d}t}=\frac{2}{3}{\left(\frac{3}{16}\right)}^{1/3}(G{M}_{\star }{)}^{1/3}({ϵ}_{\mathrm{D}}Z{)}^{2/3}{t}^{-1/3},$$

$${\tau }_{\mathrm{f}}=\frac{\sqrt{3}}{8}\frac{{ϵ}_{\mathrm{P}}}{\eta }\frac{{\Sigma }_{\mathrm{peb}}(r_P)}{{\Sigma }_{\mathrm{g}}(r_P)},$$

$$\eta =-\frac{1}{2}{\left(\frac{H}{r}\right)}^2\frac{\partial \ln P}{\partial \ln r}.$$

$${\Sigma }_{\mathrm{peb}}=\sqrt{\frac{2S_{\text{peb}}{\dot{M}}_{\mathrm{peb}}{\Sigma }_{\mathrm{g}}}{\sqrt{3}\pi {ϵ}_{\mathrm{P}}{r}_{\mathrm{P}}{v}_{\mathrm{K}}}},$$

$${\tau }_{\mathrm{f},\mathrm{rock}}=\frac{{\rho }_{\mathrm{rock}}{\mathrm{Rocky}}_{\mathrm{peb}}}{{\Sigma }_{\mathrm{g}}},$$

$${\Sigma }_{\mathrm{peb},\mathrm{rock}}=\frac{{\dot{M}}_{\mathrm{peb},\mathrm{rock}}}{2\pi r{v}_{\mathrm{peb},\mathrm{rock}}},$$

$$v_{r_{\mathrm{peb},\mathrm{rock}}}=-\frac{2.0\eta v_{\mathrm{k}}}{{\tau }_{\mathrm{f},\mathrm{rock}}+{\tau }_{\mathrm{f},\mathrm{rock}}^{-1}}-\frac{\nu }{r},$$

$$f_{\mathrm{acc}}=\frac{{\dot{M}}_{\mathrm{core}}}{{\dot{M}}_{\mathrm{peb}}}.$$

$$M_{\mathrm{iso}}=25f_{\mathrm{fit}}{\mathrm{M}}_{\mathrm{E}}+\frac{{\Pi }_{\mathrm{crit}}}{\lambda }{\mathrm{M}}_{\mathrm{E}},$$

$$f_{\mathrm{fit}}={\left[\frac{H/r}{0.05}\right]}^3\left[0.34{\left(\frac{\log ({\alpha }_3)}{\log (\alpha )}\right)}^4+0.66\right]\left[1-\frac{\frac{\partial \ln P}{\partial \ln r}+2.5}{6}\right],$$

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

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

$$x=-\frac{\partial ln\text{leavevmode}{\Sigma }_{gas}}{\partial ln\text{leavevmode}r},\text{leavevmode}\text{leavevmode}\text{leavevmode}\beta =-\frac{\partial ln\text{leavevmode}T}{\partial ln\text{leavevmode}r}.$$

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

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

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

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

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

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

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

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

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

Full text

1 Introduction

Exoplanet systems present a diversity of architectures compared with the structure of our home planetary system. Planets with sizes between those of Earth and Neptune, that is, with between 1 and 4 Earth radii, have been found in compact multi-planet systems orbiting their host stars at orbital periods shorter than 100 days (Lissauer et al., 2011b; Marcy et al., 2014; Fabrycky et al., 2014). These systems are typically referred to as hot super-Earth systems and their high abundance (e.g., Mayor et al., 2011; Batalha et al., 2013; Howard, 2013; Fressin et al., 2013) is one of the greatest surprises in exoplanet science. Observations and planet occurrence studies suggest that at least 30% of the FGK-type stars in the galaxy host hot super-Earths with a period of less than 100 days (Mayor et al., 2011; Howard et al., 2012; Fressin et al., 2013; Petigura et al., 2013; Zhu et al., 2018; Mulders, 2018; Mulders et al., 2018). Hot super-Earths are found by inference to have orbits with low orbital eccentricities and mutual inclinations (Mayor et al., 2011; Lissauer et al., 2011b; Johansen et al., 2012; Fang & Margot, 2012; Xie et al., 2016; Zhu et al., 2018).

Our understanding of the origins of hot super-Earths remains incomplete. Many models have been proposed in the last decade or so, but these scenarios have been gradually refined by observational constraints and simulations and many have already been discarded (see discussions in Raymond et al., 2008; Raymond & Morbidelli, 2014; Raymond & Cossou, 2014; Schlichting, 2014; Morbidelli & Raymond, 2016; Ogihara et al., 2015a; Chatterjee & Tan, 2015; Izidoro et al., 2017). We now briefly discuss three scenarios: (a) in-situ accretion ; (b) drift-then-assembly; and (c) migration.

The in-situ scenario for the formation of hot super-Earths proposes that hot super-Earths formed more or less where they are seen today. This scenario requires (i) large amounts of mass in solids in the inner protoplanetary disk in order to form massive planets and (ii) the late formation of super-Earths in gas-poor or gas-free environments in order to avoid gas-driven migration (Raymond et al., 2008; Chiang & Laughlin, 2013; Raymond & Cossou, 2014; Schlichting, 2014; Schlaufman, 2014; Dawson et al., 2015, 2016). The required high density of solids in the protoplanetary disk is problematic if one imposes the canonical dust/gas ratio, because it would imply that the disk is gravitationally unstable. Thus, in-situ formation models often assume an enhanced dust/gas ratio (e.g., Hansen & Murray, 2012, 2013; Hansen, 2014; Dawson et al., 2015, 2016; Lee & Chiang, 2016), but without modeling how such an enhancement occurred or its properties (e.g., its radial extent). The late formation of super-Earths is sometimes justified by invoking that, as long as there is a significant amount of gas in the disk, gas dynamical friction damps the eccentricities of growing planets and prevents mutual orbit crossing. However, numerical simulations including all dynamical effects show that, whenever the density of planetesimals is large, planet growth is inexorably fast (Ogihara et al., 2015a); thus a “late formation” seems implausible. An alternative possibility is that planetesimal formation itself occurred late, but, given the assumed large dust/gas ratio, this also appears implausible because the streaming instability becomes effective as soon as the metallicity is larger than a few percent (Yang et al., 2017). If super-Earths form too early in the gas-rich inner parts of the disk, fast inward type-I migration leads to very compact systems of hot super-Earths that fail to match observations (Ogihara et al., 2015a). Effects of disk winds have been invoked to avoid large-scale fast inward type-I migration (Ogihara et al., 2015b). Inward type-I migration is only completely suppressed if the the corotation torque is assumed to be fully unsaturated, which may not be realistic for Earth-mass planets (Ogihara et al., 2018). This suggests that some level of inward gas-driven migration is difficult to avoid.

The drift-then-assembly model proposes that small particles such as pebbles or planetesimals drift inwards by gas drag and pile up at the pressure bump that may form at the transition between the magnetorotational instability-active inner regions and the exterior dead zone (Chatterjee & Tan, 2014; Boley & Ford, 2013; Boley et al., 2014; Chatterjee & Tan, 2015; Hu et al., 2016). This collection of particles forms a ring of solid material that eventually becomes gravitationally unstable and collapses to form a planet. The newly created planet induces another pressure bump outside its orbits and the process repeats. Although this idea is interesting, the model remains to be further developed. For instance, it has been shown that Hall shear instability may dominate near the disk inner edge when the magnetic field is aligned with the spin vector of the disk. This instability may prevent the local formation of close-in super-Earths via the processes envisioned in the drift-then-assembly model (Mohanty et al., 2018). Ueda et al. (2019) also showed that the accumulation of dust at the disk inner edge occurs only if the dust grains are large enough and the viscosity in the dead-zone is sufficiently low.

The migration model proposes that super-Earths or their constituent planetary embryos formed in gas-rich disks and migrated inward from outside their current orbits by planet–disk gravitational interaction (Terquem & Papaloizou, 2007; Ida & Lin, 2008, 2010; McNeil & Nelson, 2010; Hellary & Nelson, 2012; Cossou et al., 2014; Coleman & Nelson, 2014, 2016; Izidoro et al., 2017; Ogihara et al., 2018; Raymond et al., 2018a; Carrera et al., 2018). Simulations modeling planet–disk interaction predict that hot super-Earths typically migrate inwards and pile up at the disk inner edge forming long chains of first-order mean motion resonances (Terquem & Papaloizou, 2007; Raymond et al., 2008; McNeil & Nelson, 2010; Rein, 2012; Rein et al., 2012; Horn et al., 2012; Ogihara & Kobayashi, 2013; Cossou et al., 2014; Raymond & Cossou, 2014; Ogihara et al., 2015a; Liu et al., 2015, 2016; Izidoro et al., 2017; Ormel et al., 2017; Unterborn et al., 2018; Ogihara et al., 2018). During the gas disk phase, the orbital eccentricities and inclinations of super-Earths are tidally damped by the gaseous disk (Papaloizou & Larwood, 2000; Goldreich & Sari, 2003; Tanaka & Ward, 2004; Cresswell & Nelson, 2008; Bitsch & Kley, 2010, 2011a; Teyssandier & Terquem, 2014). As the disk evolves and loses mass, eccentricity and inclination damping become less efficient. Once the gas dissipates, these effects vanish. If eccentricities and orbital inclinations of planets in the chain grow because of mutual interactions, the orbits of the planets may eventually cross each other leading to collisions and scattering events (e.g., Kominami & Ida, 2004; Iwasaki & Ohtsuki, 2006; Matsumoto et al., 2012; Morbidelli, 2018). This evolution typically breaks the resonant configurations established during the gas disk phase, leading to a phase of giant impacts (for a detailed discussion on the onset of dynamical instabilities in resonant systems, see Meyer & Wisdom (2008, and references therein); Goldreich & Schlichting (2014, and references therein); Deck & Batygin (2015, and references therein); Xu et al. (2018, and references therein); Pichierri et al. (2018, and references therein); Pichierri & Morbidelli (2020, and references therein). The final (post-instability) configuration of such systems is nonresonant. Thus, the fact that most super-Earths are not found in resonant systems (Lissauer et al., 2011c; Fabrycky et al., 2014) should not be used as an argument against the migration model. In fact, the current distributions of super-Earths are consistent with all systems emerging from resonant chains. Systems like Kepler-223 (Mills et al., 2016) and TRAPPIST-1 (Gillon et al., 2017; Luger et al., 2017) have multiple-planet resonant chains that are naturally produced by migration. These resonant chains represent the small fraction of systems that did not become unstable after gas dispersal (Cossou et al., 2014; Izidoro et al., 2017; Ogihara et al., 2018).

The migration model can match the period ratio distribution of Kepler planets by combining a fraction of unstable and stable systems, typically $\lesssim 10\%$ of stable and $\gtrsim 90\%$ of unstable systems (Izidoro et al., 2017). The model also suggests that the large number of Kepler systems with single transiting planets versus multiple transiting planets – known as the Kepler dichotomy (Johansen et al., 2012) – is a consequence of the dispersion of orbital inclinations of super-Earths rather than a true dichotomy in planetary multiplicity. Dynamical instabilities excite the orbital inclinations of planets such that observations are likely to miss transits of mutually inclined planets (Winn, 2010), raising the number of single-transiting systems. Although some studies claim that the Kepler dichotomy is real  (e.g., Fang & Margot, 2012; Johansen et al., 2012; Moriarty & Ballard, 2016), the migration-instability model is arguably the simplest explanation for an apparent dichotomy, predicting an insignificant number of real single-planet systems in the Kepler sample.

To date, a downside of the migration model has been that simulations have assumed from the beginning that several Earth-mass planetary embryos formed in different parts of the disk. However, it remains unclear as to whether such distributions of embryos could really have arisen naturally, or as to how the details of initial conditions affect the final systems. It is crucial to evaluate the legitimacy of these assumptions and more importantly to assess whether the migration model remains viable when a more self-consistent approach is used.

The goal of this paper is to revisit the migration model (Cossou et al., 2014; Izidoro et al., 2017) and to build a comprehensive scenario for the origins of super-Earths that is consistent with a broad picture of planet formation. This is the main upgrade of this study relative to Izidoro et al. (2017). A key new ingredient in our scenario is pebble accretion. Pebble accretion plays a role after the formation of planetesimals, which are themselves thought to form by clumping of drifting pebbles via the so-called streaming instability  (Youdin & Goodman, 2005; Johansen et al., 2009; Carrera et al., 2015; Simon et al., 2016; Carrera et al., 2017). Planetesimals then grow by mutual collisions, and once they approach lunar mass, they start to accrete pebbles efficiently  (Johansen & Lacerda, 2010; Ormel & Klahr, 2010; Lambrechts & Johansen, 2012; Johansen et al., 2015; Xu et al., 2017). Pebble accretion can explain the rapid growth of the building blocks of terrestrial planets, super-Earths, and ice giants (e.g., Lambrechts & Johansen, 2012; Levison et al., 2015, 2015; Bitsch et al., 2015; Bitsch & Johansen, 2016; Ndugu et al., 2018; Chambers, 2016; Johansen et al., 2015; Johansen & Lambrechts, 2017; Lambrechts et al., 2019). However, what sets the destiny of planets in becoming either hot super-Earths or a different class of planet is not entirely clear.

This paper is part of a trilogy that develops a unified model to explain the formation of rocky Earth-like planets, hot super-Earths, and giant planets from pebble accretion and migration. This paper is dedicated to the formation pathways, dynamical evolution, and compositions of hot super-Earths. The other two companion papers of this trilogy focus on (1) the formation of terrestrial planets and super-Earths inside the snow line, highlighting the role of the pebble flux (Lambrechts et al., 2019, hereafter referred to as Paper I), and (2) understanding the conditions required for gas giant planet formation in the face of orbital migration (Bitsch et al., 2019, a companion paper of this series, hereafter referred to as Paper III)

Paper I models planetary growth exclusively inside the snow line. It shows that sufficiently low pebble fluxes – integrated pebble fluxes of $114{\rm M_{\oplus}}$ – lead to the slow growth of protoplanetary embryos that do not migrate substantially during the gas disk lifetime. These embryos are typically the same mass as Mars at the end of the gas disk phase. An increased pebble flux by a simple factor of two bifurcates the evolution of these systems inducing the formation of more massive rocky planetary embryos that migrate inwards and pile up at the inner edge of the disk. These different growth histories separate the formation of truly Earth-like planets from that of rocky super-Earths. The long-term dynamical evolution of these systems reveals that dynamical instabilities after gas dispersal finally set the architecture of these systems. Dynamical instabilities among small Mars-mass planetary embryos result in collisions that lead to the formation of Earth-like planets of no more than 4 Earth masses. Instabilities among large rocky planetary embryos near the disk inner edge assemble rocky hot super-Earth systems. An extensive analysis of the formation of hot super-Earths, also accounting for their possible origins beyond the snow line, is not performed in Paper I, but is shown here in the present work. Finally, Paper III shows that if the pebble flux is large enough (e.g., integrated pebble flux of $350-480M_{\oplus}$) super-Earths turn into gas giant planets. Paper III presents a self-consistent model of the growth and migration of gas giant planets. It also highlights the role of migration in the formation of gas giants, an aspect that is typically ignored in simulations modeling the formation of the Solar System from pebble accretion (Levison et al., 2015; Chambers, 2016).

This trilogy of papers is designed to provide a comprehensive view of planet formation and evolution, revealing the possible broad diversity of planetary systems produced from pebble accretion, disk evolution, migration, and long-term dynamical evolution of planetary systems.

The current paper is structured as follows. Section 2 describes the methods, namely our gas disk model, and pebble accretion prescription. Section 3 presents simulations designed to understand the role of the pebble flux in the formation and early evolution of the system. We also tested the role of the initial distribution of protoplanetary embryos and pebble sizes inside the snow line. In Section 4 we discuss our results in light of those of Paper I. In Section 5 we present the results of our simulations dedicated to modeling the growth and long-term dynamical evolution of hot super-Earth systems. In Section 6 we lay out observational constraints that we use to evaluate the success of different models in matching the true super-Earth distributions. We focus on the period ratio distribution, the Kepler “dichotomy” and the compositions (rocky vs. icy) of super-Earths. In Section 7 we place the Solar System in the context of our model. In Section 8 we present our conclusions. In Appendix A we detail our prescription for gas-driven migration.

2 Method

Our simulations are performed with FLINTSTONE, our new N-body code built on MERCURY (Chambers, 1999). In addition to the standard orbital integration, we include modules to consistently compute the following:

• •

the evolution of an underlying, gas-dominated protoplanetary disk;

• •

the gas-assisted accretion of pebbles onto growing planetary cores; and

• •

the orbital migration and tidal damping of cores.

Other aspects of the MERCURY code were left intact. For instance, collisions are modeled as inelastic merging events that conserve linear momentum. We also do not model volatile loss during giant impacts. Thus, the final water and ice content of planets in our simulations should be interpreted as upper limits. We also do not perform modeling of planetary interior. In this work we use the term “ice” indiscriminately as a proxy for all physical states of water in our planets. In all our simulations, planetary objects are ejected from the system if they reach heliocentric distances of greater than 100 AU.

We compared FLINTSTONE with the code used in our companion paper (Paper I) by Lambrechts et al. (2019) which is built on Symba (Duncan et al., 1998) and found similar results for test problems regarding planet migration, pebble accretion, and damping of eccentricity and inclination.

Below we describe our disk model and our prescriptions for pebble accretion, gas-driven migration, inclination, and eccentricity gas-tidal damping.

2.1 Gas disk model

Our underlying disk is represented by 1D radial profiles derived from 3D hydrodynamical simulations modeling gas disk evolution (Bitsch et al., 2015).111We note that the gas disk model considered in this work is more sophisticated than the one from in Paper I The hydrodynamical model accounts for the effects of viscous heating, stellar irradiation, and radial diffusion.

In the standard alpha-disk paradigm, the gas accretion rate on the young star is written as

[TABLE]

where ${h=H/r}$ is the disk aspect ratio with H being the disk scale height, $\alpha$ the dimensionless viscosity parameter (Shakura & Sunyaev, 1973), r the distance to the star, ${\Omega_{k}=\sqrt{GM_{\odot}/r^{3}}}$ the orbital the keplerian frequency, and $\Sigma_{gas}$ the gas surface density.

Following Hartmann et al. (1998) and Bitsch et al. (2015), the relationship between the disk and star age and the gas accretion rate onto the star is given by

[TABLE]

Finally, the hydrostatic equilibrium yields

[TABLE]

where ${\rm t_{disk}}$ is the disk age, $T$ is the temperature at the gas disk midplane, ${\mu}$ is the gas mean molecular weight set to ${\rm 2.3\leavevmode\nobreak\ gmol^{-1}}$, $G$ is the gravitational constant, ${M_{\odot}}$ is the mass of the star set equal to 1 solar mass, and $\mathcal{R}$ is the ideal gas constant. We use ${t_{{\rm start}}}$ to represent the disk age at the starting time of our simulations.

From Eqs. 1, 2, and 3 we obtain the gas surface density ($\Sigma_{{\rm gas}}$) of the disk by using the fits of the disk temperature at the midplane provided in Bitsch et al. (2015). The disk metallicity and ${\rm\alpha}$-viscosity parameter are set to 1% and to $\alpha=5.4\times 10^{-3}$, respectively. The main source of viscosity in protoplanetary disks is relatively unknown (e.g., Turner et al., 2014). Recent ALMA observations (e.g., Pinte et al., 2016; Dong et al., 2017; Zhang et al., 2018) have suggested disk viscosities that can be one order of magnitude lower than the value assumed in our simulations (however, see Dullemond et al. (2018) for a substantially higher radial diffusion). Lower disk viscosities could increase the efficiency of pebble accretion because the pebble scale height depends on viscosity. At the same time, in sufficiently low-viscosity disks lower mass planets would be able to open up gaps in the disk. It would be interesting to test the impact of the disk viscosity on our model but this remains a subject for future investigations.

As the disk evolves, we recalculate the disk structure every 500 years rather than at every time-step to save computation time. This does not significantly impact the quality of our approach because the disk structure only significantly changes on long timescales ($\gg$500 years).

Protoplanetary disks are also expected to have inner cavities in their gas density distribution created due to the magnetic star–disk interaction (Koenigl, 1991). The dipolar magnetosphere of a rotating young star may disrupt the very inner parts of the protoplanetary disk dictating the gas accretion flow onto the stellar surface (Romanova et al., 2003; Bouvier et al., 2007; Flock et al., 2017). The truncation radius is probably at a few stellar radii, inside the co-rotation radius with the star and where the magnetic field pressure balances the pressure of the accreting disk. In standard disks with typical magnetized young stars the disk truncation radius is expected to be around $\sim$0.03-0.2 AU (e.g., Bouvier, 2013).

The inner cavity of a protoplanetary disk is expected to have an important impact on planet formation because it is likely to act as an efficient planet trap, preventing inwardly migrating planets from simply falling onto the star (Masset et al., 2006; Romanova & Lovelace, 2006; Romanova et al., 2019). The innermost planets in several Kepler systems indeed have orbital periods corresponding to the expected truncation radius of disks, corroborating the existence of disk inner edges (e.g., Mulders et al., 2018).

To account for an inner cavity, our disk inner edge is fixed at $\sim 0.1$ AU in our nominal simulations unless stated otherwise. In order to represent the drop in surface density at this location we multiply the gas density by the following factor

[TABLE]

where $r$ is the heliocentric distance and $r_{{\rm in}}$ is the location of the disk inner edge. This approach has also been used in previous works (Cossou et al., 2014; Izidoro et al., 2017).

In Appendix A we describe how planet migration is modeled in our simulations. We emphasize that the migration prescription considered in this work is more sophisticated than that of Paper I. Unlike Paper I, here we take into account corotation (entropy and vortensity driven) effects in order to obtain a more quantitatively realistic model.

2.2 Model setup

Our simulations start with a distribution of small planetary embryos with masses randomly chosen between 0.005 and 0.015 Earth masses. Each simulation starts with a slightly different distribution of planetary seed masses. For this mass regime, pebble accretion typically dominates over planetesimal accretion (Johansen & Lambrechts, 2017). In all our simulations, the initial orbital inclination of planetary embryos is randomly selected between 0 and 0.5 degrees. The orbital eccentricities of all planets are set initially to $10^{-4}$. Other orbital angles were randomly selected between 0 and 360 degrees.

The birth location of the first planetesimals is unconstrained by observations. Some models suggest that planetesimal formation is likely to initiate just outside the snow line (Dra̧żkowska & Dullemond, 2014; Armitage et al., 2016; Dra̧żkowska & Alibert, 2017; Carrera et al., 2017; Schoonenberg & Ormel, 2017; Dra̧żkowska & Dullemond, 2018). Other models suggest planetesimal formation takes place just inside the snow line (Ida & Guillot, 2016) or even around 1 AU (Dr\każkowska et al., 2016). In our simulations we test different distributions of seeds where they naturally account for planetesimal formation: only throughout the inner disk (inside the snow line), only throughout the outer disk (outside the snow line), and both inside and outside the snow line (throughout the disk). In our disk model the snow line moves inward as the disk evolves; at 0.5 Myr, it is around 3 AU but by 3 Myr is already around 0.7 AU. The details of our different initial distributions of protoplanetary embryos are presented in Table LABEL:tab:model. In simulations of Models I and II, planetary seeds are initially distributed past 0.7 AU. Because at ${\rm t_{disk}=3\leavevmode\nobreak\ Myr}$ the snow line is at 0.7 AU, simulations with ${\rm t_{start}=3\leavevmode\nobreak\ Myr}$ correspond to scenarios where planetary seeds formed only outside the snow line.

The formation time of the first planetesimals is also poorly constrained. It has been proposed that in our inner Solar System at least two distinct generations of planetesimals were born: one forming early at about 0.5 Myr after the so-called calcium–aluminium-rich inclusions (CAIs; Villeneuve et al., 2009; Kruijer et al., 2012) and others late at about 3 Myr after CAIs. However, it is not clear whether this scenario is a generic outcome of planet formation. Given our limited understanding of the timing of planetesimal formation, we explore in our simulations two end-member scenarios. For simplicity, we consider a single generation of planetesimals for all our simulations. Our first scenario corresponds to the case where planetesimals form very early. In this case our simulations start with ${\rm t_{start}=0.5\leavevmode\nobreak\ Myr}$. In the second scenario, planetesimals are assumed to form late and our simulations start with ${\rm t_{start}=3.0\leavevmode\nobreak\ Myr}$ (see also Paper III). We note that different ${t_{{\rm start}}}$ imply different disk structures at the beginning of our simulations and this has an important impact on the system evolution both in terms of planet migration and pebble accretion (Bitsch et al., 2015).

In our simulations, planetary seeds starting initially inside the snow line are assumed to be rocky while those outside are considered to have 50% of their mass as ice. Rocky and icy planetary seeds have bulk densities of 5.5${\rm\leavevmode\nobreak\ g/cm^{3}}$ and 2${\rm\leavevmode\nobreak\ g/cm^{3}}$, respectively.

2.3 Pebble accretion

As in Paper I we do not model drifting pebbles as individual particles because of the high computational cost (but see Kretke & Levison, 2014; Levison et al., 2015). Instead, pebbles in the disk are modeled as a background field, which depends on the flux and Stokes number of the pebbles, which evolve over time. Our modeling of the pebble flux is more sophisticated than that from Paper I, where the number of pebbles in the disk decays exponentially during the disk lifetime and the Stokes number of the pebbles is fixed. Here, the pebble flux and Stokes number are computed in the context of dust coagulation models (Birnstiel et al., 2012; Lambrechts & Johansen, 2012, 2014), as summarized below. We follow the prescription of pebble accretion by protoplanets from Johansen et al. (2015), which can also account for reduced accretion rates for eccentric and inclined bodies.

The accretion rate onto the planetary core is given as

[TABLE]

where $\rho_{\rm p,mid}$ is the midplane density of pebbles, related to the pebble surface density layer $\Sigma_{\rm peb}$ via

[TABLE]

where $H_{\rm peb}=H_{\rm gas}\sqrt{\alpha_{\rm set}/\tau_{\rm f}}$ is the pebble scale height (Youdin & Lithwick, 2007), $\tau_{\rm f}$ is the Stokes number, and $\alpha_{set}$ is the dimensionless viscosity parameter representing disk midplane turbulence which determines the vertical settling of pebbles. In our nominal simulations, $\alpha_{\rm set}=\alpha$, although some studies argue for $\alpha_{\rm set}\ll\alpha$ (Zhu & Stone, 2014). $R_{\rm acc}$ denotes the accretion radius of the planet, $\delta v=v_{\rm rel}+\Omega R_{\rm acc}$ with $v_{\rm rel}$ being the relative velocity between pebbles and planets. To calculate the mass accretion rate, we first need the accretion radius and the pebble density averaged over the accretion radius. The stratification integral $\bar{S}$ is defined as the mean pebble density normalized by the pebble density in the midplane, namely:

[TABLE]

where $z_{0}$ is the height over the disk midplane where the protoplanet is located along its orbit. This equation is derived by considering layers of constant $z$ and consequently constant pebble density, but there is no analytical solution to this integral and the solution has to be obtained numerically. However, Johansen et al. (2015) used a square approximation that integrates the pebble density over a square instead of a circle rendering the integral analytically solvable, which we use here. The exact solution is shown in the Appendix of Johansen et al. (2015).

The pebble flux is calculated as

[TABLE]

Here, $r_{\rm g}$ represents the heliocentric distance at which dust particles have grown to pebble sizes and start drifting inwards by gas drag, and $\Sigma_{\rm g}(r_{\rm g})$ is the gas surface density at the pebble production line location $r_{\rm g}$ (for more details see Lambrechts & Johansen (2014)). The quantity $Z$ represents the dust to gas ratio in the disk that can be converted into pebbles at the pebble production line $r_{\rm g}$ at time $t$. In our nominal model, we take $Z=1\%$ . Lambrechts & Johansen (2014) derived the time-dependent radial location of the pebble production line as

[TABLE]

and thus

[TABLE]

where $M_{\star}$ is the stellar mass, which we set to $1M_{\odot}$, $G$ is the gravitational constant, and $\epsilon_{\rm D}=0.05$ is associated with the logarithmic growth range from dust grains to pebble sizes. In our model, at 0.5 Myr, 3 Myr, and 5 Myr the pebble production line is at ${\sim 77{\rm\leavevmode\nobreak\ AU}}$, ${\sim 250{\rm\leavevmode\nobreak\ AU}}$, and ${\sim 350{\rm\leavevmode\nobreak\ AU}}$, respectively222Our disk is probably on the large side of the spectrum of typical radial disk sizes ($\sim$100 AU), but some disks are as large as 500-1000 AU (e.g., Hughes et al., 2008). We note that in the drift-limited pebble-growth model, the pebble flux depends on the gas surface density at the pebble production line (Eq. 8). As the pebble production line moves beyond $\sim$50 AU, the gas disk and drift-limited pebble-growth models assumed in this work can strongly underestimate the pebble column density (Bitsch et al., 2018a) compared to observations (Wilner et al., 2005; Carrasco-González et al., 2016). This is a consequence of the low gas surface density in the remote regions of our evolving protoplanetary disk. In our nominal disk with ${S_{{\rm peb}}=1}$, the total mass in pebbles between 50 and 90 AU is only about 10${{\rm\leavevmode\nobreak\ M_{\oplus}}}$ at 1 Myr. Observations suggest that the total mass in dust and pebbles in the outer parts (50-90au) of the HL Tau disk, which is about 1 Myr old, is between 80 and 280 ${{\rm\leavevmode\nobreak\ M_{\oplus}}}$ (Carrasco-González et al., 2016). On the other hand, some observed disks seem not to be massive enough to form the observed exoplanet population (Manara et al., 2018; Mulders, 2018). To account for a possible range of disk masses, we rescale the gas surface density at the pebble production line by a factor ${S_{{\rm peb}}}$in order to increase the pebble surface density (or decrease in some extreme cases; see eq. 13). We note that $S_{{\rm peb}}$ is treated as a free-parameter in our model.

Particles in the gas disk grow and drift, so the local Stokes numbers of the particles are limited by drift (Birnstiel et al., 2012; Lambrechts & Johansen, 2012). By equating the pebble growth timescale with the drift timescale (i.e., valid for the drift-limited dominant pebble size; see section 2.4 in Lambrechts & Johansen (2014)) the stokes number may be written as

[TABLE]

where $\Sigma_{\rm g}(r_{\rm P})$ and $\Sigma_{\rm peb}(r_{\rm P})$ are the gas and pebble surface densities at the location of the planets $r_{\rm P}$. We represent the radial pressure support of the disk through the dimensionless parameter

[TABLE]

The pebble flux decreases in time as the disk evolves (Birnstiel et al., 2012; Lambrechts et al., 2014; Bitsch et al., 2018a). The (evolving) pebble surface density $\Sigma_{\rm peb}$ can be calculated using the pebble flux and Stokes number (Eqs. 8 and 11) as

[TABLE]

where $r_{\rm P}$ denotes the semi-major axis of the planet, and $v_{\rm K}=\Omega_{\rm K}r_{\rm P}$. ${S_{{\rm peb}}}$ is a nondimensional linear scaling factor of the pebble flux (see Eq. 8). In Eq. 13, $\epsilon_{\rm P}=0.5$ represents the pebble sticking efficiency under the assumption of near-perfect sticking (Lambrechts & Johansen, 2014). ${S_{{\rm peb}}=1}$ corresponds to an integrated pebble flux ${\rm I_{{\rm peb}}\approx 194\leavevmode\nobreak\ M_{\oplus}}$ beyond the snow line (see Figure 1). It remains possible that a disk with a higher(lower) pebble flux could be simply associated with a disk with higher(lower) metallicity (see Eq. 11 of Lambrechts & Johansen (2014)). For simplicity, we assume a unique disk metallicity set equal to 1% during the entire disk lifetime to model planet migration in all simulations. The same approach was taken in our companion paper by Bitsch et al. (2019).

We assume that the water component of pebbles crossing the snow line sublimates, releasing the rocky or silicate counterpart in the form of smaller dust grains. As the snow line moves inwards, the boundary between big (icy) and small (rocky) pebbles moves with the snow line. The Stokes number of icy pebbles is given by Eq. 11 in the drift-limited solution (dominant pebble size). We do not use Eq. 11 to calculate the Stokes number of silicate pebbles in the inner disk. Instead, we set the size of silicate pebbles to sizes of 1 mm which correspond to the typical chondrule sizes of ordinary chondrites (Friedrich et al., 2015) in our nominal simulations. We then calculate the stokes number of 1mm silicate pebbles in the inner disk via the classical definition of Stokes number (instead of using Eq. 11) as

[TABLE]

(e.g., Lambrechts & Johansen, 2012) where ${\rm Rocky_{{\rm peb}}}$ and $\rho_{\rm rock}$ define the size and bulk density of silicate pebbles, respectively. The bulk density of silicate pebbles is assumed to be 5.5${\rm g/cm^{3}}$ in all simulations. Inside the snow line, the pebble flux is reduced by a factor of two to account for the sublimation of the water component of the pebbles ($\dot{M}_{\rm peb,rock}=0.5\dot{M}_{\rm peb}$). We calculate the pebble surface density inside the snow line using the following equation.

[TABLE]

where the radial velocity of rocky pebbles is given by (e.g., Ormel & Klahr, 2010; Bitsch et al., 2018a)

[TABLE]

where ${\rm\nu=\alpha H^{2}\Omega_{k}}$ is the gas kinematic viscosity.

In our nominal simulations, we set ${\rm Rocky_{{\rm peb}}}$=1 mm which corresponds to the size of chondrules, but we also test the effects of larger silicate pebbles in our models. These pebble sizes were already used in Morbidelli et al. (2015) to reproduce the different growth speeds of the terrestrial planets compared to the gas giants in the Solar System. On the other hand, although laboratory experiments (e.g., Windmark et al., 2012; Blum, 2018) and numerical models (Birnstiel et al., 2010; Banzatti et al., 2015) suggest that the growth of silicate pebbles beyond milimeter size is challenging, centimeter-sized chondrule clusters are found in unequilibrated ordinary chondrites (Metzler, 2012). Silicate centimeter(cm)-sized pebbles probably also existed in our early inner Solar System at least to some (small) level. Thus, in our simulations we also analyze the effects of considering 1cm pebbles inside the snow line. We have verified that in our disk model, at the very early stages of the disk (${\rm t_{disk}\approx 0.0\leavevmode\nobreak\ Myr}$), 1cm (1mm) pebbles inside 0.5 AU would be in the Stokes (Epstein) regime of gas–particle coupling (Epstein, 1924). However, as our simulations start with ${\rm t_{start}=t_{disk}=0.5\leavevmode\nobreak\ Myr}$ or 3 Myr, both 1mm and 1cm silicate pebbles beyond $\sim$0.2 AU (the starting location of our innermost seeds) are in the Epstein regime of gas-particle coupling. The top panel of Figure 2 shows the gas and pebble surface densities as a function of orbital distance as the disk evolves. The bottom panel of Figure 2 shows the Stokes number and sizes of pebbles in our disk with nominal pebble fluxes and sizes at different disk ages.

To account for the sublimation of the water component of pebbles crossing the ice line we assume a reduction of the pebble mass fluxes to half. This is consistent with the assumed composition of seeds forming inside and outside the snow line. In our disk model, the water ice line moves inwards as the disk dissipates, and reaches 1 AU at 2 Myr.

A planet only accretes a fraction $f_{\rm acc}$ of the flux of pebbles $\dot{M}_{\rm peb}$ drifting past its orbit:

[TABLE]

The pebble flux arriving at interior planets is thus reduced by this fraction $f_{\rm acc}$, reducing the accretion rates onto the interior planets.

A sufficiently massive planet opens a partial gap and creates an inversion in the radial pressure gradient of the disk, halting the inward drift of pebbles (Paardekooper & Mellema, 2006a; Morbidelli & Nesvorny, 2012; Lambrechts et al., 2014). The mass at which this takes place is called the pebble isolation mass. The pebble isolation mass in itself is a function of the local properties of the protoplanetary disks, namely the viscosity, aspect ratio, and radial pressure gradient, and of the Stokes number of the particles, which can diffuse through the partial gap of the planet (Bitsch et al., 2018b; Weber et al., 2018). Here, we follow the exact description of Bitsch et al. (2018b), who give the pebble isolation mass including diffusion as

[TABLE]

with $\lambda\approx 0.00476/f_{\rm fit}$, $\Pi_{\rm crit}=\frac{\alpha}{2\tau_{\rm f}}$, and

[TABLE]

where $\alpha_{3}=0.001$.

After planets reach the pebble isolation mass, they can start to accrete gas from the disk. That process is modeled in our companion paper by Bitsch et al. (2019). Here we assume that the contraction of the gaseous envelope (Piso & Youdin, 2014) is sufficiently slow333A sufficiently high envelope opacity may prevent fast contraction (Piso & Youdin, 2014). that our seeds would not transition into runaway gas accretion and stay at super-Earth mass. Three-dimensional hydrodynamical simulations show that only planets more massive than about 15 ${\rm M_{\oplus}}$ and in disks with opacities lower than $\sim$0.01 $cm^{2}/g$ can transition to runaway gas accretion (Lambrechts & Lega, 2017). In this paper, we neglect the effects of gas accretion in all simulations but we flag a growing planet as a giant planet core when it reaches a mass larger than 15 ${\rm M_{\oplus}}$ during the gas disk phase.

3 The role of the pebble flux

Our simulations were conducted in two phases. First, we ran simulations considering a wide range of pebble fluxes (${S_{{\rm peb}}=0.2}$, 0.4, 1, 2.5, 5, and 10; see Table LABEL:tab:model). We used the outcome of these simulations to inspect which pebble fluxes could lead to the formation of planets in the super-Earth/Neptune mass range – with masses smaller than $\sim$15 ${\rm M_{\oplus}}$ – during the gas disk phase. The second phase of long-term integrations beyond the gas disk phase is discussed in Section 5. We remind the reader that our goal is to model the formation of hot super-Earth systems; the formation of rocky terrestrial planets and rocky super-Earths as well as giant planets are modeled in companion papers of this trilogy (Lambrechts et al., 2019; Bitsch et al., 2019). In Paper I we presented our findings that an integrated pebble flux of ${\rm 114\leavevmode\nobreak\ M_{\oplus}}$ leads to the formation of classical terrestrial planets while simulations with integrated pebble fluxes of ${\rm 190\leavevmode\nobreak\ M_{\oplus}}$ (or more) form super-Earth systems. As discussed before, one should not expect that by considering the same integrated pebble flux as Paper I our simulations will produce the same types of planets (super-Earths or terrestrial planets). This is because our planets may grow outside the snow line by accreting larger pebbles. We discuss this issue again below.

We performed ten simulations for each value of the pebble flux scaling ${S_{{\rm peb}}}$. We did not model the long-term dynamical evolution of these systems. We stop our simulations at the end of the gas disk phase, namely at 5 Myr. Due to the large number of simulations to be conducted in this section and in order to save central processing unit (CPU) time, we also set the disk inner edge in these simulations to about ${\rm r_{\rm in}\simeq 0.3\leavevmode\nobreak\ AU}$. The location of the disk inner edge in simulations modeling the formation of hot super-Earth system have been typically assumed to be around 0.1 AU (Cossou et al., 2014; Izidoro et al., 2017; Ogihara et al., 2018; Lambrechts et al., 2019). We likewise set the disk inner edge to ${\rm r_{in}\simeq 0.1\leavevmode\nobreak\ AU}$ in our simulations of Section 5. However, here we set ${\rm r_{in}\simeq 0.3\leavevmode\nobreak\ AU}$ and use a larger integration time-step to conduct this large batch of simulations without degrading the quality of our results. In this section, we infer the final compositions of the planets by tracking the source of the accreted material in terms of icy and silicate pebbles.

3.1 Model I

Figure 3 shows the growth and dynamical evolution of protoplanetary embryos in one simulation of Model I (see Table LABEL:tab:model for the definition of our models) with ${\rm t_{start}=0.5\leavevmode\nobreak\ Myr}$ and ${S_{{\rm peb}}=1}$ during the gas disk phase. In this simulation, planetary embryos grow by pebble accretion more quickly outside the snow line than inside because pebbles are typically larger in the cold regions of the disk in our model. At 1 Myr, the largest planetary embryo outside the snow line is about ${\rm 1\leavevmode\nobreak\ M_{\oplus}}$. At 2 Myr, the mass of the largest planetary embryo is about ${\rm\sim 7\leavevmode\nobreak\ M_{\oplus}}$. As the disk evolves, the pebble isolation mass (black dashed line) decreases across the entire disk because the disk cools down and gets thinner (Bitsch et al., 2015). The gray curves give the boundary of the (a,M) region where migration is outwards. We note that the pebble isolation mass is within the range of the region of outward migration in some parts of the disk.

Some planetary embryos in the simulation from Fig. 3 grew massive enough to enter the outward migration region, yet they did not always migrate outwards. This is because the mutual gravitational interaction with other growing planetary embryos acts to excite eccentricities and to reduce the contribution of the co-orbital torque responsible for driving outward migration (see Bitsch & Kley, 2010; Cossou et al., 2013). Outer nearby embryos may also act as a dynamical barrier for embryos in the outward-migration region. As the disk evolves further, the outward migration region quickly shrinks. Typically, the outermost embryo in the outward migration zone is eventually caught by the inward-moving outer edge of the outward-migration region. We do not see a significant level of outward migration in these simulations (see Figure 4). Instead, planetary embryos inside the outward migration region typically migrate very slowly inwards in a chain of mean-motion resonances. In this configuration, the outermost embryo sits at the edge of the outward migration region (see panel corresponding to 2 Myr in Figure 4).

Dynamical instabilities and collisions take place during the migration phase. As embryos in the outer disk grow fast and migrate inwards, they rapidly join the chain of embryos trapped at the outer edge of the outward migration region. This strong convergent migration tends to destabilize the resonant chain and promote additional collisions between embryos and further growth beyond the pebble isolation mass (Wimarsson et al., 2020). As the outward-migration region shrinks further, the more massive planetary embryos eventually get out of the region, becoming free to quickly migrate inwards.

At 3 Myr, several embryos have reached the inner edge of the disk to form a long resonant chain in mutual first-order mean motion resonances. At this time, some planets within 0.5 AU lie inside the outward-migration region. These planets have been pushed inwards by planets quickly migrating inwards with masses that fall above the outward-migration region. Figure 3 shows that at the end of the gas disk phase, planetary embryos in the resonant chain at the inner edge of the disk have masses lower than ${\rm\sim 10\leavevmode\nobreak\ M_{\oplus}}$.

The final compositions of all planets in the simulation from Fig. 3 are dominated by ice-rich material. Although this simulation starts with small planetary embryos inside the snow line (see snapshot corresponding to ${\rm t=0.5\leavevmode\nobreak\ Myr}$ in Figure 3) the growth of planetary embryos beyond the snow line is much more efficient. Planetary embryos growing beyond the snow line quickly reach pebble isolation mass blocking the pebble flux to inner regions, and consequently starving the innermost planetary embryos, in particular the rocky ones. Moreover, as larger planetary embryos migrate inwards they either collide with or scatter away small rocky planetary embryos. Consequently, the final system of close-in planets has a water-rich composition.

Figure 4 shows the mass and orbital evolution of the simulation from Fig.  3. In Fig. 4, the evolutions of the planets with final orbital semi-major axes within ${\rm 0.7\leavevmode\nobreak\ AU}$ are shown in color. We highlight the innermost objects of the system because we are interested in the formation of close-in super-Earths. The Kepler sample is almost complete for transiting planets larger than Earth and orbital periods smaller than 200 days. Thus, we compare our results with Kepler observations taking into account only planets with orbital periods shorter than about 200 days. As our simulations consider a solar-mass central star, this corresponds to planets with an orbital radius smaller than about 0.7 AU. In Figure 4, embryos that underwent collisions or were ejected from the system are shown in bold black. The gray color is used to show planetary objects on final orbits outside 0.7 AU and also leftover planetary embryos.

Figure 4 shows how the orbital eccentricities of embryos on average increase from the start of the simulation to about 1 Myr as they grow by pebble accretion and consequently their mutual gravitational interaction becomes stronger. The gas tidal effects damp the orbital eccentricity and inclination and counter-balance the effects of mutual gravitational stirring. Orbits of larger protoplanetary embryos are more efficiently damped by the gas. Figure 4 shows that only after ${\rm 1\leavevmode\nobreak\ Myr}$ planetary embryos have reached masses large enough to enter in a regime where type-I migration is reasonably fast. Migration leads to resonant shepherding and scattering events among planetary embryos. As consequence of close encounters, leftover planetary embryos were typically scattered by the largest migrating bodies when the latter moved towards the disk inner edge. Scattered bodies tend to reach dynamically excited orbits which may not be efficiently damped during the gas lifetime to allow residual growth by pebble accretion and migration (Levison et al., 2015). At the end of the gas disk phase, this simulation produced six planets with masses larger than ${\rm 2\leavevmode\nobreak\ M_{\oplus}}$ inside 0.7 AU. As shown in the panel illustrating mass evolution, the first phase of growth of these final planets is characterized by pebble accretion. However, collisions with other embryos, which also take place during early times (e.g., ${\rm\sim 1-1.5\leavevmode\nobreak\ Myr}$), become more common when the planets approach the inner edge of the disk. Of course, this is a consequence of convergent migration but also an effect of short dynamical timescales in the inner parts of the disk. At 3 Myr, multiple planetary embryos have reached the inner edge of the disk forming a long resonant chain.

Figure 5 shows the evolution of a simulation like the one in Figure 4 but for a case where the pebble flux is 2.5 times higher. Overall, the dynamical behavior of protoplanetary embryos in this simulation is similar to that shown in Figure 4. However, as expected, a higher pebble flux promotes a much faster growth of protoplanetary embryos by pebble accretion. In this simulation, during the first ${\rm\sim 1.5\leavevmode\nobreak\ Myr}$ protoplanetary embryos have already reached masses of greater than ${\rm\sim 6\leavevmode\nobreak\ M_{\oplus}}$ by pure accretion of pebbles. We note that, broadly speaking, this corresponds to the final masses of planetary objects at the end of the gas disk phase in the simulation of Figure 4. As these larger planetary objects migrate inwards, they collide with each other and grow even further. Most of these collisions happen inside the disk inner cavity as this region becomes overcrowded because of the successive arrival of planetary embryos migrating from more distant regions of the disk. At the end of the gas disk phase, the most massive planetary embryo in this simulation is ${\rm\sim 50\leavevmode\nobreak\ M_{\oplus}}$. Thus, the final masses of planetary objects in the simulations of Figures 4 and 5 are drastically different.

In our simulations we neglect gas accretion onto protoplanetary embryos. The three most massive final planets produced in the simulation of Figure 5 – with masses larger than ${\rm\sim 20\leavevmode\nobreak\ M_{\oplus}}$ – represent very good candidates for accretion of massive gas atmospheres to become gas giants. We do not model the formation of gas giant planets in this paper, but we dedicated a companion paper (Bitsch et al., 2019) to address this issue more carefully. In this work, whenever a planetary embryo becomes more massive than ${\rm\sim 15\leavevmode\nobreak\ M_{\oplus}}$ during the gas disk phase we flag it as a giant planet core rather than a super-Earth. We use the simulations of this section only to identify the pebble fluxes that lead to the formation of super-Earths rather than giant planet cores. Once we know which pebble fluxes produce planets of ${\rm\leq 15\leavevmode\nobreak\ M_{\oplus}}$ during the gas disk phase, in Section 5 we use this information to conduct a larger set of simulations dedicated to modeling the formation and long-term dynamical evolution of super-Earths. However, we note that our assumed threshold mass of ${\rm\sim 15\leavevmode\nobreak\ M_{\oplus}}$ for transition to runaway gas accretion (our flag of giant planet core) depends on dust opacity in the envelope, which is uncertain (Pollack et al., 1996; Ormel et al., 2015; Lee & Chiang, 2015; Lambrechts & Lega, 2017).

Figure 6 shows the final masses of planetary embryos in all simulations of Model I with scaling pebble fluxes varying from ${S_{{\rm peb}}=0.2}$ to 2.5. Each panel of Figure 6 shows the outcome of ten simulations in a diagram depicting semi-major axis versus mass. Figure 6 shows a clear trend: the final planet masses are higher for higher pebble fluxes (larger ${S_{{\rm peb}}}$). Planetary embryos growing in low-pebble-flux environments do not grow massive enough to migrate to the disk inner edge. We recall that in these particular simulations, planetary embryos are initially distributed from 0.7 to 20 AU. For ${S_{{\rm peb}}=0.2}$ and ${S_{{\rm peb}}=0.4}$ the amount of radial migration of planetary embryos is typically modest. In these two cases, most planetary embryos remain at sub-Earth mass and beyond 0.5 AU (see also Paper I for simulations showing the long-term dynamical evolution of a similar population of rocky embryos). Earth-mass or more-massive planets at 1-2 AU produced for ${S_{{\rm peb}}=0.2}$ and 0.4 (for Model I) are planets that started forming farther out and migrated down to their final position. We compare our results with those of Paper I in the following section. Before discussing the details of these results we also recall that the integrated pebble flux is the flux of pebbles during the entire course of the simulation. A simulation with ${S_{{\rm peb}}=0.2}$ and ${\rm t_{start}=0.5\leavevmode\nobreak\ Myr}$, for example, features an integrated pebble flux of ${\rm\sim 0.2\times 194\leavevmode\nobreak\ M_{\oplus}=39\leavevmode\nobreak\ M_{\oplus}}$. A simulation with ${S_{{\rm peb}}=1}$ and ${\rm t_{start}=0.5\leavevmode\nobreak\ Myr}$ results in an integrated pebble flux of ${\rm\sim 1\times 194\leavevmode\nobreak\ M_{\oplus}}$ (see Figure 1). In both cases, only a fraction of the integrated pebble flux is used to build planets. In Figure 6 these numbers are 17%, 25%, 32%, and 30% for $S_{\rm peb}=0.2$, 0.4, 1, and 2.5, respectively.

Figure 6 also shows that our nominal pebble flux ${S_{{\rm peb}}=1}$ results in the formation of planets that do indeed migrate to the inner edge of the disk (assumed at ${\rm\sim 0.3\leavevmode\nobreak\ AU}$ in these simulations). Planets inside ${\rm 0.7\leavevmode\nobreak\ AU}$ in simulations with ${S_{{\rm peb}}=1}$ have masses lower than ${\rm\sim 10-15\leavevmode\nobreak\ M_{\oplus}}$ (see also Figure 4). A further factor of 2.5 increase in the pebble flux results in the formation of multiple planets with masses as large as ${\rm\sim 40-50\leavevmode\nobreak\ M_{\oplus}}$ (see also Fig. 5). We note that planets do not necessarily stay beyond the disk inner edge. Planets migrating to the inner edge of the disk pile up in long resonant chains. In some cases, the innermost planets anchored at the inner edge of the disk are pushed inside the disk cavity (Cossou et al., 2014; Izidoro et al., 2017; Brasser et al., 2018; Carrera et al., 2019). Some planets are pushed to distances as close as 0.06 AU from the star.

Our results also show that a simple difference of a factor of approximately two in pebble flux (from ${S_{{\rm peb}}=1}$ to 2.5) is enough to bifurcate the evolution of our planetary systems from one leading to a system of super-Earth mass planets to one hosting several massive protoplanetary cores which are very likely to become gas giants. For example, several planets forming in our simulations with ${S_{{\rm peb}}=2.5}$ are more massive than the Solar System ice giants by a factor of a few. Although not shown in Figure 6, for completeness, we also inspected the results of simulations considering higher scaling pebble fluxes of ${S_{{\rm peb}}=5}$ and 10. In these simulations, the final planets reaching the inner edge of the disk are as massive as ${\rm 100\leavevmode\nobreak\ M_{\oplus}}$ due to convergent migration towards the disk inner edge and successive collisions. Bitsch et al. (2019) present the results of simulations with higher pebble fluxes where gas accretion onto planetary cores is consistently modeled.

3.2 Model

II

Figure 7 shows the distribution of planets produced in simulations of Model II, which is similar to Model I but the seeds are assumed to appear only at ${\rm t_{start}=3.0\leavevmode\nobreak\ Myr}$, i.e., towards the end of the lifetime of the disk (5 Myr). For the nominal pebble flux scaling ${S_{{\rm peb}}=1}$, the integrated drifted pebble flux for a disk with ${\rm t_{start}=0.5\leavevmode\nobreak\ Myr}$ is about ${\rm\sim 194\leavevmode\nobreak\ M_{\oplus}}$. In a disk with ${\rm t_{start}=3.0\leavevmode\nobreak\ Myr}$ and ${S_{{\rm peb}}=1}$ the total integrated pebble flux is about ${\rm\sim 30\leavevmode\nobreak\ M_{\oplus}}$ (see Figure 1). This is the total reservoir of mass available for planetary growth.

We note that the results from Models I and II shown in Figures 6 and 7, respectively, are not dramatically different: with ${S_{{\rm peb}}=0.2}$ and ${\rm t_{start}=0.5\leavevmode\nobreak\ Myr}$ versus ${S_{{\rm peb}}=1}$ and ${\rm t_{start}=3.0\leavevmode\nobreak\ Myr}$. For example, the largest planets in the top left-hand panel of Figure 7 are around 1 AU and have masses of between 1 and ${\rm 2\leavevmode\nobreak\ M_{\oplus}}$. Planets around 1 AU in the top left-hand panel of Figure 6 have also masses in this range. The reason for this is that the total drifted pebble flux is about the same in these two setups (${\rm\sim 39\leavevmode\nobreak\ M_{\oplus}}$ for Model-I with ${S_{{\rm peb}}=0.2}$ and ${\rm\sim 30\leavevmode\nobreak\ M_{\oplus}}$ for Model II with ${S_{{\rm peb}}=1.0}$; note that they have different ${\rm t_{start}}$). However, even though the integrated pebble fluxes are relatively similar, the Stokes numbers of the pebbles (see Eq. 11) are different because of the different pebble and gas surface densities at the planetary location, which can lead to different growth. This effect becomes more clear in the outer regions of the disk when comparing, for example, the top right panel of Figure 6 (${S_{{\rm peb}}=0.4}$) with the top right panel of Fig. 7 (${S_{{\rm peb}}=2.5}$). We note that these two cases have also comparable integrated pebble fluxes (${\rm\sim 75\leavevmode\nobreak\ M_{\oplus}}$). In the top right-hand panel of Figure 6, the typical final mass of planetary embryos around 10 AU is sub-Mars mass while those in top right-hand panel of Figure 7 have masses of about ${\rm\sim 2\leavevmode\nobreak\ M_{\oplus}}$, at least a factor of 20 larger.

Figure 7 shows that increasing pebble flux from ${S_{{\rm peb}}=2.5}$ to 5 is enough to promote the delivery of planetary embryos to the inner edge of the disk (see bottom left-hand panel of Figure 7). The final planets anchored at the inner edge of the disk in this case have typical masses of a few Earth masses to Neptune mass. A further increase in the pebble flux with ${S_{{\rm peb}}=10}$ (see bottom left-hand panel of Figure 7) promotes the formation of planetary embryos with masses greater than ${\rm\sim 20\leavevmode\nobreak\ M_{\oplus}}$. This reinforces our previous finding of Figure 6 where a simple increase in the pebble flux by a factor of two is enough to bifurcate the growth pattern of planetary embryos from one leading to super-Earths or hot-Neptunes to another producing massive planetary cores which would be very likely to become gas giants (Lambrechts & Lega, 2017).

3.3 Model

III

In Model III we test an extreme scenario in which planetary seeds are assumed to have only formed in the innermost rocky regions of the disk (see Table LABEL:tab:model). Our goal is not only to inspect the final masses of the planets as a function of the pebble flux but also to study how the initial distribution of planetary embryos may impact the final planet compositions. We note that all our simulations of Models I and II failed to grow and deliver rocky planetary embryos to the disk inner edge. We naively expect this to be more easily achieved if the initial distribution of planetary embryos is restricted to the region inside the snow line.

Figures 8 and 9 show the planets produced by simulations of Model III with pebbles of different sizes inside the snow line. In Figure 8 the size of the silicate pebble is ${\rm Rocky_{{\rm peb}}=0.1\leavevmode\nobreak\ cm}$ while in Figure 8 this is ${\rm Rocky_{{\rm peb}}=1\leavevmode\nobreak\ cm}$. In both cases ${\rm t_{start}=0.5\leavevmode\nobreak\ Myr}$.

In simulations of Figure 8, small planetary embryos are initially distributed between $\sim 0.3$ and 2 AU. The disk snow line is at about 3 AU at 0.5 Myr. Consequently, all planetary embryos grow at first by the accretion of silicate pebbles. However, as the disk evolves, the snow line moves inwards eventually sweeping the outermost seeds from outside in. Thus, the outermost seeds (at about 2 AU) eventually switch to accreting larger icy pebbles until they reach pebble isolation mass, and thus block the flux of pebbles from the outer disk, starving the inner disk. We note that the final ice mass fraction of a planetary seed initially inside the snow line is primarily regulated by the total mass in icy dust particles available in the outer disk at the time when the snow line crosses the seed’s orbit (Ida et al., 2019). Nevertheless, the pebble flux filtering by outer seeds (if they exist) also plays a crucial role. Figure 8 shows that protoplanetary seeds swept up by the snow line may reach masses large enough to move by type-I migration to the inner edge of the disk before the gas is dispersed. However, at the end of the gas disk phase the very final shape of the outward migration region is such that it favors ${\rm\sim 2\leavevmode\nobreak\ M_{\oplus}}$ protoplanetary embryos becoming stranded between 1 and 2 AU (see e.g., top left-hand panel of Fig. 8 and the shape of the migration zone in Fig. 3).

Figure 8 shows that simulations with ${S_{{\rm peb}}=1}$ failed to form multiple planets with masses larger than a few ${\rm M_{\oplus}}$ of any composition anchored at the inner edge of the disk. Increased pebble fluxes (${S_{{\rm peb}}=2.5}$, 5, and 10) successfully promote the formation and delivery of icy planetary embryos with masses of ${\rm\sim 5-10\leavevmode\nobreak\ M_{\oplus}}$ to the inner edge of the disk at 0.3 AU. However, only simulations with ${S_{{\rm peb}}=10}$ (botton right-hand panel of Fig. 8) produce final rocky protoplanetary embryos with masses larger than ${\rm\sim 1\leavevmode\nobreak\ M_{\oplus}}$, i.e., in the super-Earth mass range. Higher pebble fluxes also tend to produce a larger number of closer-in icy planetary embryos compared to lower pebble fluxes because initially more distant seeds (swept by the snow line) can grow faster and to larger masses, consequently leaving the outward migration region and migrating inwards.

Figure 9 shows the planets produced in our Model III simulations with larger silicate pebbles (${\rm Rocky_{{\rm peb}}=1\leavevmode\nobreak\ cm}$) for different pebble fluxes. As expected, comparing the results of Figure 8 and Figure 9 it is clear that the growth of rocky planetary embryos is far more efficient when ${\rm Rocky_{{\rm peb}}=1\leavevmode\nobreak\ cm}$ for all pebble fluxes. Interestingly, only simulations with ${S_{{\rm peb}}\gtrsim 2.5}$ produced concomitantly and systematically rocky and icy super-Earths of similar mass. Low pebble fluxes tend to favor the formation of large icy super-Earths as opposed to rocky ones. Simulations with ${S_{{\rm peb}}=10}$ produce at least a few planetary embryos as large as ${\rm\sim 20\leavevmode\nobreak\ M_{\oplus}}$. Results presented in Figs. 8 and 9 clearly show that avoiding the formation of icy super-Earths is a difficult task even in the scenario where seeds only form well inside the snow line. Figures 8 and 9 show that Model-III only produces systems dominated by rocky super-Earths if seeds starting inside 1 AU grow to Earth-mass or larger sufficiently quickly, allowing them to grow above pebble-isolation mass. In order to ensure fast growth of these seeds we invoked the existence of cm-sized silicate pebbles in the inner disk (${\rm Rocky_{{\rm peb}}=1\leavevmode\nobreak\ cm}$). However, faster growth could be equally achieved with our nominal mm-sized silicate pebbles if we reduce the level of turbulent vertical stirring of mm-sized silicate pebbles. We performed two simple simulations to confirm this claim. We find that the growth of a ${\rm 0.01\leavevmode\nobreak\ M_{\oplus}}$ planetary seed at 1 AU growing in a 0.5 Myr-old disk where ${S_{{\rm peb}}=5}$ and ${\rm Rocky_{{\rm peb}}=1\leavevmode\nobreak\ cm}$ is very similar to that of an identical seed in a simulation with ${S_{{\rm peb}}=5}$, ${\rm Rocky_{{\rm peb}}=0.1\leavevmode\nobreak\ cm,}$ and a disk where the scale height of the pebble layer is smaller by a factor of approximately three. This scale height corresponds to a disk where $\alpha_{set}\simeq\alpha/10$ (see Section 2.3). We discuss this issue again in Section 6.

The main difference between the results of Model III and Models I and II is in the efficiency of pebble accretion; in other words, the fraction of the integrated pebble flux converted into planets. In Models I and II – depending on ${S_{{\rm peb}}}$ and ${t_{{\rm start}}}$ – from 17% to 49% of the pebble flux is converted into planets while in Model III this quantity varies between 1.4% and 5.6%. This represents a significant reduction. We now list a few possible reasons for this difference but probably not all. First, in Model III (both Figures 8 and 9), seeds starts only inside the snow line where pebbles are typically smaller than beyond the snow line (even when ${\rm Rocky_{{\rm peb}}=1.0\leavevmode\nobreak\ cm}$) and, consequently, the efficiency of pebble accretion is reduced. Second, in simulations of Model-III with ${\rm Rocky_{{\rm peb}}=0.1\leavevmode\nobreak\ cm}$ the outermost seed is typically swept by the moving snow line and grows really fast to pebble isolation mass suppressing planetary growth in the inner regions. Finally, the pebble isolation masses in the inner regions of the disk are small thus planets stop growing at lower masses.

Although the parameters ${S_{{\rm peb}}}=5\leavevmode\nobreak\ {\rm and}\leavevmode\nobreak\ 10$ result in a large mass in pebbles available for planetary growth (if ${t_{{\rm start}}}=0.5\leavevmode\nobreak\ {\rm Myr}$ they imply ${\rm\sim 970\leavevmode\nobreak\ M_{\oplus}}$ and ${\rm\sim 1940\leavevmode\nobreak\ M_{\oplus}}$, respectively), we stress that the actual amount of mass in pebbles required to grow our planetary systems is always much smaller. Most of the planets in our simulation are already fully formed after a relatively short time, especially the ones migrating to the inner system. Thus, to form most of these planets it is only required to have a sufficiently large pebble flux in the beginning of our simulations and not necessarily during the whole gas disk lifetime.

The pebble accretion efficiency in our simulations varies between 1 and 50%. We note that in simulations with single planets, the reported pebble accretion efficiencies are typically lower, only a few percent per planet (e.g., Lin et al., 2018; Liu & Ormel, 2018). Our simulations start with many protoplanetary seeds, so pebble accretion becomes far more efficient simply because multiple planets can accrete pebbles.

4 Comparison with the results of Paper I

Paper I defined super-Earths as planets that reached a mass large enough for significant migration during the gas disk lifetime, and showed that rocky super-Earths tend to be clustered close to the inner edge of the disk, particularly at the end of the gas-disk phase (some spreading occurring during instabilities after gas removal; see section 4 below). The rocky planets at ${\rm\sim 1\leavevmode\nobreak\ AU}$ are typically terrestrial planets, that is, planets that grew mostly after gas removal from a disk of planetary embryos too small to migrate (less than 3 Mars masses). The results presented here in Figs. 6-9 show that super-Earths (in the sense of Paper I) can also be at 1 AU or even significantly beyond this region. Their existence is also suggested by micro-lensing observations (Beaulieu et al., 2006; Muraki et al., 2011; Bennett et al., 2007; Kubas et al., 2012; Sumi et al., 2010, 2016; Koshimoto et al., 2017). However, in our model these planets are never rocky. They either start their formation beyond the snow line and migrate inwards towards their final position, or are swept by the snow line (in the case of model III), so that they accrete most of their mass from icy pebbles and experience a temporary phase of outward migration. In this paper, rock-dominated planets are seen only in model III and their distribution mirrors the results of Paper I: only small rocky embryos remain near 1 AU (e.g., Fig. 8 upper-left panel) while rocky super-Earths migrate towards the inner edge of the disk (here placed farther out than in Paper I for the reasons explained in section 2).

The sensitive dependency of the final planet mass on the pebble flux discussed in paper I is recovered here. However, a quantitative comparison with Paper I shows some apparent differences. For instance, in Paper I a total pebble flux of 200 Earth masses (5/3$\times$nominal) leads to the formation of a super-Earth close to the disk’s inner edge. Here, a pebble flux of 194 Earth masses (Fig. 8) leads only to small (typically submartian) rocky planetary embryos. However, the pebble flux reported in Paper I is the flux within the snow line. The flux reported here is beyond the snow line, and should be divided by two inside of the snow line because of ice sublimation. Moreover, pebble size, pebble flux, pebble scale height, and the gas disk model differ from the model presented in Paper I. Even when the pebble fluxes are truly equivalent, the final masses of planets in our simulations and those in Paper I will probably differ. In Paper I, the disk snow line location is fixed throughout the entire lifetime of the disk. Thus, because in our simulations the disk snow line moves inside 1 AU at late stages, the pebble sizes around 1-2 AU may become much larger than the fixed pebble size considered in Paper I, resulting in different growth modes. Moreover, in our simulations, a significant fraction of the pebble flux is filtered by the growing planets in the outer disk. Hence, the rocky pebble flux seen by the embryos in the inner disk in the simulation of Fig. 8 (upper-left panel), for example, corresponds to a subnominal flux in Paper I. There is therefore no inconsistency on the final masses of rocky bodies.

The main difference between this paper and Paper I is that we simulate the concurrent growth of rocky and ice-rich super-Earths —as well as dynamical perturbations between them— whereas the existence of icy planets was neglected in Paper I. Moreover, we show that, if seeds are present throughout the outer disk, the growth of rocky planets is precluded because almost all of the pebble flux is intercepted or blocked by the outer, icy embryos. Only if the seed distribution ends near the snow line (model -III) can the growth of rocky planets be significant. If the rocky pebbles are smaller than the icy pebbles, the ice-rich planets and the smaller, rocky-dominated planets are radially mixed (Fig. 8). If the rocky pebbles are as big as icy pebbles, rocky-dominated planets are typically the innermost ones and ice-rich planets are those farther out (Fig. 9). Although both cases can be consistent with the observations of extra-solar planets (see Sections 5 and 6 below), they are inconsistent with the Solar System, where rocky planets are much smaller than the ice-rich planets (Uranus, Neptune, and the cores of Jupiter and Saturn) but the two types of planets are not radially mixed. A discussion of the Solar System case is deferred to Section 7. Paper I can therefore be seen as a subcase where some process, for example the formation of giant planets (see also Paper III), impedes the migration of ice-rich planets into the inner system.

5 Formation of close-in super-Earths

We now aim to model the formation of super-Earth systems like those observed. As shown in the previous section, different combinations of parameters (${t_{{\rm start}}}$, ${S_{{\rm peb}}}$, and ${\rm Rocky_{{\rm peb}}}$) produce dramatically different planetary systems. To conduct our new simulations we have purposely selected pebble fluxes that can successfully lead to the formation of hot super-Earth systems rather than systems of terrestrial planets (see Lambrechts et al., 2019) or gas giants (Bitsch et al., 2019). We are interested in setups where at the time of the gas disk dispersal most planets anchored at the disk inner edge have masses of between ${\rm\sim 1\leavevmode\nobreak\ M_{\oplus}}$ and ${\rm\sim 15\leavevmode\nobreak\ M_{\oplus}}$. This makes them reasonably consistent with the expected masses for most close-in super-Earths observed by Kepler (e.g., Wolfgang et al., 2016).

To systematically evaluate the performance of the migration model we performed 250 simulations considering four different selected scenarios. These are shown in Table LABEL:tab:2, which also presents the integrated pebble flux of each setup, which depends on ${S_{{\rm peb}}}$ and ${t_{{\rm start}}}$.

To remain consistent with previous simulations of the formation of close-in super-Earth systems (e.g., Izidoro et al., 2017) we set the disk inner edge at $\rm r_{in}$=0.1 AU rather than at $\sim$0.3 AU. This also matches the location of the actual inner edge of the disk according to Mulders (2018). This is essentially the only difference between the setup of the simulations presented in this section and those presented above. However, to have better statistics when analyzing the long-term dynamical evolution of planetary systems, for each scenario of Table LABEL:tab:2 we performed 50 simulations with slightly different initial conditions. As before, each model is represented by an initial distribution of planetary embryos as described in Table LABEL:tab:model. After the gas disk dispersal, simulations are continued for an additional $\sim$50 Myr. Some particularly interesting cases were integrated up to 300 Myr.

5.1 Long-term dynamical evolution

Figure 10 shows the evolution of a simulation of Model II where ${\rm t_{start}=3.0\leavevmode\nobreak\ Myr}$, ${S_{{\rm peb}}=5}$, and ${\rm Rocky_{{\rm peb}}=0.1\leavevmode\nobreak\ cm}$. The dynamical evolution during the gas disk phase of the growing planetary embryos is qualitatively similar to those in Figures 4 and 5. Essentially, planetary embryos grow and migrate inwards establishing a long chain of planets mutually captured in first-order mean motion resonances. At the end of the gas disk phase, planetary embryos anchored at the inner edge of the disk have masses lower than ${\rm\sim 7\leavevmode\nobreak\ M_{\oplus}}$. After gas dispersal, the simulation is continued for an additional 45 Myr in a gas-free scenario but the resonant chain remains dynamically stable. The right-hand panels of Figure 10 show the orbital eccentricities and inclinations of all planetary embryos in this simulation. As shown, planetary embryos inside 0.7 AU at the end of the simulation have orbits with very low eccentricities and inclinations.

Figure 11 shows the evolution of another simulation of Model II in which $t_{{\rm start}}={\rm 3.0\leavevmode\nobreak\ Myr}$, $S_{{\rm peb}}={\rm 1}$, and ${\rm Rocky_{{\rm peb}}}$=0.1 cm. Again, planetary embryos grow and migrate to the disk inner edge forming a long resonant chain. The long-term dynamical evolution of planets in this simulation is in great contrast with that of Figure 10. After the gas disk phase (see vertical gray line in the figure), the resonant chain of planets anchored at the inner edge becomes dynamically unstable at about 6 Myr. The dynamical instability promotes collisions and further planetary growth. The instability phase leads to a final planetary system which is dynamically more excited. Planets’ orbital eccentricities reach values of a few percent while orbital inclinations increase by a few degrees.

The dynamical evolutions presented in Figures 10 and 11 are representative of all our simulations. Some planetary systems present a phase of dynamical instability after gas dispersal but some do not. Following Izidoro et al. (2017), we refer to these two classes of planetary systems as “stable” and “unstable”. The fraction of stable and unstable systems varies in our simulations. We stress that once a system becomes dynamically unstable, the duration of the instability phase is typically short and the system rapidly evolves to a less compact but typically long-term stable configuration.

Figure 12 shows selected stable (left panel) and unstable (middle panel) planetary systems produced in our different models. The right-hand panel shows selected observed planetary systems. Planets in stable systems have masses lower than $\sim$10${\rm M_{\oplus}}$. There is no clear radial mass ranking in planetary systems produced in our simulations (in contrast, see Figure 3 of Ogihara et al. (2015a)). In the following section we take a closer look at the orbital architecture of these systems. Finally, we note from comparing the left and middle panels of Figure 12 that unstable systems are relatively more spread and have fewer but larger planets. The results of Figure 12 are qualitatively similar to those in Izidoro et al. (2017) where the formation of close-in super-Earth systems is modeled assuming ad hoc initial distributions of Earth-mass planetary embryos in the outer parts of the disk.

5.2 Fraction of unstable systems and timing of the instability

Figure 13 presents a cumulative distribution of the epoch of the last collision after gas disk dispersal for different sets of simulations. Dynamical instabilities start to occur as soon as the gas goes away (at 5 Myr). In all our simulations, the fraction of dynamically unstable systems is always smaller than 1. However, some scenarios of Table LABEL:tab:2 produce a much higher rate of instabilities than others. For example, red and purple dashed lines representing Model I show that more than 90% of planetary systems anchored at the inner edge of the disk become dynamically unstable after gas dispersal. On the other hand, the fractions of unstable systems in simulations of Models II and III drop to 70% and 45%, respectively.

In the simulations of Izidoro et al. (2017), only $\sim$50% of planetary systems became dynamically unstable after gas dispersal. However, Izidoro et al. (2017) also showed that, in order to match the period ratio distribution of observations, more than 75% of the planetary systems are required to become dynamically unstable. Why so many systems become dynamically unstable after the gas dispersal remained unsolved in Izidoro et al. (2017). Figure 13 shows that in some of our systems this fraction may be as high as $\sim 95\%$. A discussion about why, in some of our models, a higher fraction of the systems become dynamically unstable compared to that of Izidoro et al. (2017) is presented in Section 5.3. In later sections we also evaluate how these simulations match other observational constraints.

5.3 Why do some systems become dynamically unstable?

A critical question remains as to what determines whether or not a given close-in planetary system will become dynamically unstable. To answer this question, we analyze the orbital architecture of our protoplanetary systems at the end of the gas disk dispersal and prior to the timing of the instability, at 5 Myr. Some planetary systems started the instability phase at the very end of the disk dispersal (e.g., at $\sim$4.9 Myr) and these systems were discarded from this analysis. Of course, because the rate of dynamically unstable systems is very high in some of our models we end up with only a few stable planetary systems. We do not include in our analysis cases that could suffer dramatically from small number statistics (e.g., stable systems of Model I).

The left-hand panel of Figure 14 shows the period ratio of unstable and stable systems at 5 Myr (before the instability). The higher rate of dynamical instability was observed for simulations of Model I (purple and red dashed lines in Figure 13). The purple dashed line of Figure 14 shows that Model I with ${\rm t_{start}=3.0\leavevmode\nobreak\ Myr}$ and ${S_{{\rm peb}}=5}$ produces the most compact planetary systems at the end of the gas disk phase. The middle panel of Figure 14 shows that these same planetary systems correspond to those with the larger number of planets anchored at the disk inner edge.

The lowest fraction of dynamically unstable systems belongs to Model III with ${\rm t_{start}=0.5\leavevmode\nobreak\ Myr}$, ${S_{{\rm peb}}=10}$, and ${\rm Rocky_{{\rm peb}}=0.1\leavevmode\nobreak\ cm}$, in which fewer than half of the planetary systems become dynamically unstable after gas dispersal (yellow dashed line in Figure 13). The period ratio distribution of planet pairs for this set of simulations shows that they are the least dynamically compact systems and also the least crowded ones, although stable and unstable systems within Model III do not show significant differences. Thus, the fate of a planetary system in becoming dynamically unstable or not after gas dispersal depends on the number of planets in the chain, compactness, and planet masses. Indeed, Matsumoto et al. (2012) showed numerically that, for a given resonant chain, the less massive the planets are, the more numerous they need to be to become dynamically unstable or, equivalently, for a given mass, there is a maximal number N of planets that can be stable. Pichierri & Morbidelli (2020) performed a detailed analytical analysis of the dynamics of resonant chains and found a theoretical justification for this observation.

6 Matching observations

We now evaluate how well our simulations match observations. In Section 6.1 we first lay out five observational constraints related to the observed super-Earth population, mainly based on observations with the NASA Kepler space telescope (Borucki et al., 2010). Our simulations already match two constraints by consistently forming super-Earth systems that have similar masses to the real ones. Next we discuss three constraints in detail and perform synthetic observations to quantitatively compare our simulations with observations. We first explore the period ratio distribution (Section 6.2) and then the Kepler dichotomy (Section 6.3). Finally, we discuss the rocky versus icy nature of super-Earths (Section 6.4) and perform several additional sets of simulations to explore the conditions required to form rocky super-Earths.

6.1 Observational constraints

To be considered successful, any super-Earth formation model must match the available observational constraints. Yet all constraints are not equally important. We now briefly discuss five constraints that we, quite subjectively, order by relative strength.

• •

The large abundance of super-Earths. At least 30%, and perhaps up to 90%, of main sequence stars host close-in super-Earths (Mayor et al., 2011; Howard et al., 2012; Fressin et al., 2013; Petigura et al., 2013; Zhu et al., 2018; Mulders, 2018; Mulders et al., 2018). Any model must explain why such planets form so readily, and are found on a wide range of stellar metallicity, with more massive super-Earths being at most slightly more abundant around high-metalicity stars (e.g., Buchhave et al., 2012; Petigura et al., 2018; Kutra & Wu, 2020).

• •

The super-Earth period ratio distribution. In multiple-planet systems, the relative spacing of planetary orbits is a measure of the dynamical state of the system. This is measured via the period ratio of adjacent planets, which has been well characterized by the Kepler mission (modulo selection effects such as missing certain planets Lissauer et al., 2011c; Fabrycky et al., 2014).

• •

Approximate masses and sizes of super-Earths. While super-Earths are typically between 1 and 4 Earth radii, determining their masses has required great investment in radial velocity (e.g., Marcy et al., 2014) and transit timing variation (e.g., Lithwick et al., 2012; Mazeh et al., 2013) measurements. This has led to the derivation of mass–radius relationships for close-in small planets (Weiss et al., 2013; Weiss & Marcy, 2014; Wolfgang et al., 2016), which suggests that most super-Earths (and mini-Neptunes) are a few to ten Earth-masses. This is a mass at which migration is very efficient.

• •

The super-Earth multiplicity distribution, or ‘Kepler dichotomy’. While a large fraction of stars are found to host super-Earths, most systems seem to host only a single super-Earth whereas a small fraction of the systems host many planets. This has been referred to as the Kepler dichotomy (Johansen et al., 2012; Ballard & Johnson, 2016a). The debate continues as to whether or not systems with a single super-Earth are truly single, and as to whether or not these systems have different origins from systems with multiple super-Earths. In our previous paper we proposed that most single super-Earths are not single and that the dichotomy is simply a consequence of the relatively broad distribution of mutual inclinations between super-Earths which reduces the probability of multiple-transiting systems (Izidoro et al., 2017). From that perspective, the dichotomy reflects two kinds of planetary systems in nature: those that underwent dynamical instabilities and inclination excitation after gas dispersal and those that avoided this. However, we note that the rate of false positive single super-Earth systems is likely to be higher than for multiple systems, potentially affecting this constraint at a quantitative level.

• •

Compositional trends among super-Earths. Atmospheric photoevaporation models suggest that many super-Earths are likely to be purely rocky. On very close-in orbits, planets are unable to retain thick envelopes in the face of strong UV-driven photoevaporation (e.g., Hubbard et al., 2007; Owen & Wu, 2013). Fulton et al. (2017) showed that there is a dip in the size distribution of close-in planets between roughly 1.5 and 2${\rm\,R_{\oplus}}$. Interior modeling of planets with inferred masses and radii suggests that the photo-evaporation valley favors mainly a rocky composition of super-Earth cores over icy cores (Owen & Wu, 2017; Jin & Mordasini, 2018; Van Eylen et al., 2019). However, at present, there is no consensus as to whether or not all super-Earths are rocky. A recent analysis of super-Earth compositions by Zeng et al. (2019) suggests that super-Earths with sizes larger than 2-3${\rm\,R_{\oplus}}$ are more likely to be water-rich worlds.

While each of these constraints is important, we consider the compositional constraints on close-in super-Earths to be the weakest. This is also because (a) the division between rocky and ice-rich is poorly defined in terms of the actual water mass fraction; and (b) there are uncertainties in both measured planetary radii and masses and in the equation of state of planetary constituents at high pressure (e.g., Valencia et al., 2007; Rogers & Seager, 2010; Swift et al., 2012; Lopez & Fortney, 2014; Duffy et al., 2015; Zeng et al., 2016; Dorn et al., 2017; Mills & Mazeh, 2017; Berger et al., 2018; Wicks et al., 2018; Smith et al., 2018). Nonetheless, we discuss the composition of super-Earths in Section 6.4.

6.2 The period ratio distribution

In this section we first examine the period ratio distribution from our simulated systems and then perform synthetic observations of those systems to statistically compare them with Kepler data.

6.2.1 Dynamical architecture of unstable and stable simulated systems

We first divide our simulations within each scenario of Table LABEL:tab:2 into groups of stable and unstable planetary systems. As above, unstable systems are those that have undergone dynamical instabilities after the gas dispersal and stable systems are those that did not present instabilities from the end of the gas disk phase (5 Myr) to the end of our simulations (typically at 50 Myr, but in some cases at 300 Myr). In order to compare our results with the sample of planet candidates from the Kepler mission we only consider in our analysis planets with semi-major axes smaller than 0.7 AU (${\rm P\lesssim 200\leavevmode\nobreak\ days}$) and with masses larger than ${\rm 1\leavevmode\nobreak\ M_{\oplus}}$ at the end of the simulation (50 Myr). These cut-offs are applied because the Kepler sample is almost complete for transiting planets larger than Earth and orbital periods smaller than 200 days (Petigura et al., 2013; Silburt et al., 2015). In our analysis, we used observational data downloaded from NASA Exoplanet archive. We selected planets with sizes between 1 and ${\rm 4\leavevmode\nobreak\ R_{\oplus}}$, stars with effective temperature between 3660 and 7600 K and, finally, we removed potential false positives, that is, planet candidates in the dataset with score parameter smaller than 0.5. In this section, we compare the results of our simulations directly with the Kepler sample. However, as observational data are expected to suffer from observational bias, in the following section we perform a more systematic analysis where we attempt to quantify the effects of observational data when comparing our synthetic planetary systems with Kepler observations.

Figure 15 shows the cumulative distributions of (a) the period ratio of adjacent planet pairs, (b) planet masses, (c) number of planets, (d) orbital eccentricities, and (e) orbital inclinations of all four different setups of Table LABEL:tab:2. Planets or planet pairs belonging to stable systems are shown with solid lines and unstable ones with dashed lines. The period ratio distribution of adjacent planet pairs shows that stable systems are dynamically very compact at the end of the gas disk phase. Most adjacent planet pairs belonging to dynamically stable systems exhibit period ratios smaller than two. stable adjacent planet pairs are typically locked in first-order mean motion resonances. Systems becoming dynamically unstable at the very end or after the gas dispersal break resonances and become dynamically much less compact. Planet pairs in unstable systems are typically not locked in first-order mean motion resonances. These results are consistent with those of Izidoro et al. (2017).

The gray line in the top left-hand panel of Figure 15 shows the period ratio distribution of adjacent Kepler planets (${\rm R<4\leavevmode\nobreak\ R_{\oplus}}$). While the period ratio distribution of stable systems is drastically different from the Kepler sample for all our models, the period ratio distributions of adjacent planets in unstable systems in the various models span from slightly more compact than observed to even more spread out (see e.g., the purple and red dashed lines in the top-left panel of Figure 15).

Figure 15 (middle-top panel) shows that stable systems have lower mass planets than unstable ones, which is expected. Dynamical instabilities promote a late phase of accretion. Indeed, at the inner edge of the disk, dynamical instabilities rarely result in ejection of planetary bodies from the system, but instead result in accretion. Typical stable systems have planets with masses lower than ${\rm\sim 10\leavevmode\nobreak\ M_{\oplus}}$ and a median value of ${\rm\sim 3-4\leavevmode\nobreak\ M_{\oplus}}$ for all scenarios of Table LABEL:tab:2. However, the mass distributions of unstable systems are dramatically different.

The blue dashed line representing unstable systems of Model II with ${\rm t_{start}=3.0\leavevmode\nobreak\ Myr}$, ${S_{{\rm peb}}=5,}$ and ${\rm Rocky_{{\rm peb}}=0.1\leavevmode\nobreak\ cm}$ produced final planets with a median mass of ${\rm\sim 7.5\leavevmode\nobreak\ M_{\oplus}}$. The red dashed line representing unstable systems of Model I with ${\rm t_{start}=0.5\leavevmode\nobreak\ Myr}$, ${S_{{\rm peb}}=1,}$ and ${\rm Rocky_{{\rm peb}}=0.1\leavevmode\nobreak\ cm}$ shows a median planet mass of about ${\rm\sim 15\leavevmode\nobreak\ M_{\oplus}}$, a factor of two larger. This result becomes easier to understand by revisiting Figure 14. Although the mass distributions of Figure 15 are similar for stable systems, the differences among the mass distributions of unstable systems are remarkable because the typical number of planets in the system in each of these models is different (see middle panel of Figure 14). Thus, even if planets produced in different scenarios have similar masses before gas dispersal, very compact systems with a larger number of planets in the resonant chain generally produce more massive planets after instabilities.

It is not straightforward to compare the masses of simulated planets with those in the Kepler sample. Kepler observations typically provide planet radii rather than masses. The masses of Kepler planets have been estimated in several studies, for example by fitting empirical mass–radius relations from well-characterized planets or by exploring probabilistic aspects of the mass–radius relation (Lissauer et al., 2011a; Fang & Margot, 2012; Weiss & Marcy, 2014; Wolfgang et al., 2016). Here we use the probabilistic mass–radius relation of Wolfgang et al. (2016) which yields ${\rm M=2.6\left(R/R_{\oplus}\right)^{1.3}}$ with a standard deviation ${\rm\sigma_{M}=\sqrt{4.41+1.5\left(R/R_{\oplus}-1\right)}}$. Because we impose a cutoff of ${\rm R<4\leavevmode\nobreak\ R_{\oplus}}$ in the Kepler sample, the maximum mass of Kepler planets inferred from Wolfgang’s mass–radius relation is ${\rm\sim 18.9\leavevmode\nobreak\ M_{\oplus}}$. With this relation, the median planet mass in our Kepler sample is ${\rm 6.4\leavevmode\nobreak\ M_{\oplus}}$. The gray lines in the top-middle panel of Figure 15 show the expected mass (${\rm M}$; thick line) and mass dispersion (${\rm M\pm\sigma_{M}}$; thin lines) distributions of Kepler planets.

The planet-mass distribution of unstable systems of model I (${\rm t_{start}=0.5\leavevmode\nobreak\ Myr}$, ${S_{{\rm peb}}=1,}$ and ${\rm Rocky_{{\rm peb}}=0.1\leavevmode\nobreak\ cm}$; red dashed line in Figure 15) shows that these planets are too massive compared to the expected masses of Kepler planets. About 20% of planets in these systems have masses larger than ${\rm 18.9\leavevmode\nobreak\ M_{\oplus}}$. Planet masses in unstable systems of model I with ${\rm t_{start}=3.0\leavevmode\nobreak\ Myr}$ and ${S_{{\rm peb}}=5}$ are marginally consistent with those expected for Kepler planets. Masses of planets in unstable systems of Models II and III are typically lower than ${\rm 18.9\leavevmode\nobreak\ M_{\oplus}}$, in good agreement with the expected masses of Kepler planets. Figure 15 also shows a clear trend. Planet pairs in Model I with ${\rm t_{start}=0.5\leavevmode\nobreak\ Myr}$ typically have larger period ratios than planet pairs composed of lower mass planets (Model I with ${\rm t_{start}=3.0\leavevmode\nobreak\ Myr}$). This result is consistent with those of Izidoro et al. (2017).

The orbital eccentricities and inclinations of planets in stable and unstable systems are also dramatically different. Planets in stable systems have low eccentricity and obits with low orbital inclination, while planets in unstable systems are dynamically excited. Statistical analyses have inferred the orbital eccentricity and inclination distributions of Kepler planets (Lissauer et al., 2011c; Kane et al., 2012; Tremaine & Dong, 2012; Figueira et al., 2012; Fabrycky et al., 2014; Plavchan et al., 2014; Ballard & Johnson, 2016b; Van Eylen & Albrecht, 2015; Xie et al., 2016; Van Eylen et al., 2019). Following Izidoro et al. (2017), these distributions are represented by Rayleigh distributions with $\sigma_{e}=0.1$ (e.g., Moorhead et al. (2011)) and $\sigma_{i}=1.5^{\circ}$ (e.g., Fang & Margot (2012)) for eccentricity and inclination, respectively. These distributions are represented by the gray lines in the top right-hand and bottom left-hand panels of Figure 15. The median orbital eccentricities of planets in unstable systems range from 0.05 to 0.1. unstable planetary systems produced in Model I are clearly the most excited ones. This result is easy to understand by inspecting the corresponding mass distributions, which are the most skewed towards massive planets. The red dashed line representing Model I in Fig. 15 also corresponds to the most massive planets. Higher mass planets have more violent dynamical instabilities and consequently their final orbital architectures are more excited (Pu & Wu, 2015; Izidoro et al., 2017). Finally, the bottom middle panel of Figure 15 shows the planet multiplicity distributions in our systems and in observations. As already discussed, stable systems have a larger number of planets than unstable systems. None of our systems match the high number of Kepler systems with a single transiting planet but we revisit this issue below when we account for the observational effects in our simulations.

6.2.2 Synthetic observations and a quantitative comparison with Kepler data

Figure 15 shows that, for some of our models, the period ratio distributions of unstable systems constitute a reasonable match to observations by themselves. Although encouraging, a more effective comparison requires an attempt to quantify the effects of observational bias. We start by first comparing the period distribution of planets in our simulations with the debiased observed period distribution from Fressin et al. (2013). Later in this section, we describe simulations of observations of our planetary systems to compare with other observed distributions (e.g., period ratio and planet multiplicity) of the Kepler sample.

Figure 16 shows the period distribution of planetary systems in each of our models (see legend of Figure 15 for the description of each model). The Kepler period distribution (thick gray line) and the distribution corrected (thin gray line) for transit probability and detection biases (Fressin et al., 2013) are also shown. While the Kepler sample shows a distribution that is shifted towards low orbital periods, the period distribution of our simulations match the debiased observations fairly well. Our Model (purple lines) provides a good match to the debiased sample, in particular for P$\gtrapprox$10 days (which corresponds to orbital distances greater than $\sim$0.09 AU for a solar mass star). It is important to keep in mind that the inner edge of our disk is at about $\sim$0.1 AU (P$\approx$11 days), so it is natural to expect that planet occurrence in our simulation will be lower inside 0.1 AU. However, the innermost planets may still be pushed inside the disk inner cavity by the outermost migrating planets, as one can see for Model I (solid red line; see also Carrera et al. (2019)). Finally, the location of the disk inner edge is a free parameter in our model and we do not include the effects of tides in this work. A disk inner edge closer to the star and stellar tides may have a significant impact on the occurrence of planets at very short orbital periods (Lee & Chiang, 2017; Carrera et al., 2019).

To compare our model with the observed distribution of period ratios between adjacent planets (and planet multiplicity distribution), we need to apply observational biases to our model. To understand why biases may be important, imagine for instance a model system of three planets with adjacent period ratios P2/P1 and P3/P2. If, for some reason, the probability of observing the middle planet is small compared that of observing each of the other two, then the observed period ratio of what appears to be adjacent planets would be P3/P1. Unlike the period distribution, it is very hard to debias the observed data on period ratios, and therefore applying the observational biases to the model (planetary systems produced in our simulations) is a preferable approach.

The transit probability of a planet in a circular orbit is ${R_{\star}/a}$, where ${R_{\star}}$ and ${a}$ are the stellar radius and planet semi-major axis, respectively. More importantly, transits are only detectable if the planet’s orbital plane is sufficiently near to the line of sight between the observer and the star. Following Izidoro et al. (2017) we simulated transit observations of the planetary systems produced in our N-body simulations. Each planetary system coming from our N-body simulations is observed from a large number of lines of sight evenly spaced by 0.1 degree from angles spanning from 30 to -30 degrees in relation to the arbitrary plane i=0. Azimuthal viewing angles are evenly spaced by 1 degree and span from 0 to 360 degrees.

For a given line of sight, if at least one planet is detected444In our synthetic observation model, a detection is characterized by a planet transiting in front of the star. We do not attempt to model the signal-to-noise ratio of the light curve, for example, and so our biases are simply due to the geometry of the system relative to the view direction. we store the orbital details of each detect planet and from the detected planet(s) we create a new observed planetary system. The top left-hand panel of Figure LABEL:fig:simulated_observations shows the period ratio distribution of stable and unstable detected systems of Model I with ${\rm t_{start}=0.5\leavevmode\nobreak\ Myr}$. In the whole of our analysis, we only consider adjacent planet pairs with ${\rm P_{out}/P_{in}}<10$. The synthetic observed distributions (black dashed/solid lines) are more compact than the respective original ones. This is because if two adjacent planets are widely spaced radially we are more likely to observe only one of them (typically the inner one) and therefore miss the information on their wide separation. Instead, planets close to each other are more likely to be observed (i.e., both of them) and therefore the datum on their narrow separation is unlikely to be missed.

This effect is more pronounced for unstable systems where adjacent planet pairs are typically not resonant but mutually inclined. Only planet pairs with sufficiently small mutual orbital inclinations are detected in unstable systems. Planet pairs with small mutual orbital inclinations are generally the more compact planet pairs as well (see Figure 15). Thus, the period ratio distribution of detected unstable systems tends to be significantly more compact than the real one. We anticipate that this result is qualitatively different from that found in Izidoro et al. (2017), where synthetic observations of unstable systems produced an even greater spread of planetary systems for the following reason.

In our simplistic synthetic simulations of transit observations, a single planetary system produced in our N-body simulations is observed from several different lines of sight. Thus, synthetic observations of a single N-body system result in several observed systems, which can contribute to the observed period ratio distribution with thousands of identical planet pairs observed from different lines of sight. One of the challenges in performing this analysis is to decide how to weight such contribution when calculating the period ratio distribution. The total number of retrieved planet pairs after combining observations from all lines of sight may vary drastically from one N-body system to another. For example, a planetary system with planets in almost coplanar orbits is probably successfully observed for several viewing angles aligned with the planets’ orbital plane, thus producing a very large number of planet pairs through the survey simulator. On the other hand, a planetary system with planets on inclined orbits is probably more rarely observed, producing a smaller number of planet pairs. Izidoro et al. (2017) computed the period ratio distribution of observed planet pairs through the survey simulator, weighting the cumulative distribution by N-body system. In other words, even if a given N-body system was observed $n$ times by the survey simulator, its characteristics were counted in the resulting distribution only once. In this work, we decide to normalize the distributions by “detected system” rather than by N-body system. Thus, the characteristics of an N-body system detected $n$ times are counted $n$ times in the resulting distributions. We believe our new approach is more appropriate than that used in Izidoro et al. (2017).

The difference between the observed and real distributions is more marginal for all other models (see middle and bottom left-hand panels in Fig. LABEL:fig:simulated_observations) because the orbital separation between adjacent planets and mutual inclinations are significantly smaller than in the case of Model I with ${\rm t_{start}=0.5\leavevmode\nobreak\ Myr}$.

Matching Kepler observations with a mix of stable and unstable systems

Figure LABEL:fig:simulated_observations shows that the detected period ratio distributions of our unstable systems can in some case match the Kepler distribution reasonably well. Nevertheless, here we attempt to statistically match Kepler observations by mixing a fraction of stable and unstable systems. As discussed in Izidoro et al. (2017), though most Kepler planet pairs are not resonant, a few known planetary systems do have planets locked in long resonant chains. This is the case of Kepler-223 (Mills et al., 2016) and TRAPPIST-1 (Gillon et al., 2017; Luger et al., 2017), for instance. Detected unstable planet pairs (from simulations) alone do not account for these long resonant chain planetary systems. In order to truly match observations, we therefore need to invoke a mixture of stable and unstable systems (Izidoro et al., 2017). We recall that we have from our N-body simulations the ratio of stable and unstable systems in each of our models (see Figure 13). Now, to compare the outcome of our simulated detections with the Kepler data we take that ratio of real stable and real unstable systems to be a free parameter whose value is determined by fitting the observations. To perform this fit, we proceed as follows.

For each assumed fraction of detected stable planet pairs we randomly select from our extensive database of stable and unstable simulated observations a number of stable and unstable planet pairs that, when added together, yields a number of planet pairs equivalent to that in the Kepler sample. For each assumed fraction of detected stable pairs we repeat this procedure 1000 times calculating the KS p-value of each subsample and Kepler observation. The p-value associated with a given fraction of detected stable planet pairs is the mean p-value computed from these 1000 subsamples. As we are fundamentally interested in assessing the real fraction of stable systems, we transform the fraction of detected stable planet pairs into an estimate of the real fraction of stable systems. We do this by dividing the fraction of detected stable planet pairs by the ratio between the mean number of planet pairs detected in stable systems and the mean number of planet pairs detected in unstable systems. This transformation is important because synthetic observations of stable systems tend to retrieve a much larger number of planets pairs than those of unstable systems.

The right-hand panels of Figure LABEL:fig:simulated_observations show the p-values for samples mixing different real fractions of stable and unstable systems. We consider that whenever the p-value of the KS-test is higher than 5% the model distribution and the observed distribution are in good statistical agreement, that is, the assumed fraction of stable systems is acceptable.

The top-right panel of Figure LABEL:fig:simulated_observations shows that our simulated observations match the Kepler sample if less than ${\rm\sim 1-5\%}$ of systems remain stable. This also holds true for simulated observations of Model I with ${\rm t_{start}=3.0\leavevmode\nobreak\ Myr}$ and Model III (right-second-from-top and bottom panels of Figure LABEL:fig:simulated_observations). However, our KS-tests show that the period ratio distribution of observed planets in Model II do not match observations for any real fraction of stable systems because its unstable systems are not sufficiently spread.

Izidoro et al. (2017) obtained a larger upper bound to the fraction of possible stable systems ($\lesssim 25\%$) than the one obtained here (${\rm\lesssim 5-10\%}$) because the p-values calculated in Izidoro et al. (2017) take a very reduced effective sample size. Moreover, we note that. in their analysis. Izidoro et al. (2017) only considered planet pairs with ${\rm P_{out}/P_{in}<3}$ and with masses lower than ${\rm 18.9\leavevmode\nobreak\ M_{\oplus}}$. We have relaxed these cutoffs in our analysis because our simulations produced lower mass planets than those of Izidoro et al. (2017) (compare Figure 15 with Figure 12 of Izidoro et al. (2017)). More importantly, we recall that the construction of the observed distributions in this work is different from that of Izidoro et al. (2017). In the new method, stable systems —which are very coplanar and consequently commonly observed— have a huge weight in the distribution. Therefore, only a few of them can be accommodated when combining stable and unstable systems to match observations.

Although the results of our study suggest that $>$95% of the systems have to become unstable after gas dispersal, only the simulations of Model I achieve this result (Figure 13). These simulations match the period ratio distribution of observations naturally, with their intrinsic fraction of unstable and stable systems. Therefore, when we consider the mixing ratio of unstable and stable systems as a free parameter we are just inferring the range of mixing ratios that would still be consistent with observations. This is different for Model II and Model III. Both Model II and Model III give a significantly smaller fraction of unstable systems, and so their period ratio distribution produced via simulated observations does not naturally match Kepler observations. This implies that for these models, some additional mechanism is needed to trigger instabilities and rise the fraction of unstable systems (e.g., that discussed in Spalding & Batygin (2016)).

The high fraction of unstable systems produced in Model I is nevertheless a significant improvement relative to Izidoro et al. (2017), where only 50% of the systems happened to become unstable after gas removal. On the other hand, we reiterate the fact that not all our models are equally successful in matching other observational constraints. For instance, model I with ${\rm t_{start}=0.5\leavevmode\nobreak\ Myr}$ and ${S_{{\rm peb}}=1}$ and model II with ${\rm t_{start}=3.0\leavevmode\nobreak\ Myr}$ and ${S_{{\rm peb}}=5}$ are the least favored models overall. The former model produces planets that are systematically too massive (top-middle panel of Figure 15) while the latter produces systems that are too dynamically compact, and do not match the Kepler period ratio distribution for any real fraction of stable systems. Therefore, our favored models are model I with ${\rm t_{start}=3.0\leavevmode\nobreak\ Myr}$ and ${S_{{\rm peb}}=5}$ and model III with ${\rm t_{start}=0.5\leavevmode\nobreak\ Myr}$ and ${S_{{\rm peb}}=10}$.

Each of our favored models has its own caveats. The masses of planets in model I with ${\rm t_{start}=3.0\leavevmode\nobreak\ Myr}$ and ${S_{{\rm peb}}=5}$ are dangerously high compared to the estimated masses of Kepler planets, but our masses should be taken as upper limits because we do not model fragmentation and volatile loss in giant impacts (Marcus et al., 2010; Stewart & Leinhardt, 2012). Also, perhaps we are simply missing some ingredient to make the fraction of unstable systems in model III higher. A possibility is that our simulations, covering just 50 Myr, are too short. On longer timescales, more systems would become dynamically unstable. Another possibility is that we are missing some physics, such as for example the interaction of the planets with planetesimals remaining in the system (Chatterjee & Ford, 2015), tidal interactions with the star (Bolmont & Mathis, 2016), or even spin-orbit misalignments effects (Spalding & Batygin, 2016). Adams et al. (2008) suggested that turbulence in the disk may prevent deep locking in resonance, increasing the post-gas instability fraction. However, Izidoro et al. (2017) found no difference between systems produced in simulations with or without turbulence. Thus, we did not attempt turbulent simulations here.

Next, we evaluate how these simulations match other observational constraints.

6.3 The Kepler dichotomy

The Kepler sample has a much larger number of planetary systems exhibiting single transiting planets than multi-transiting ones (Lissauer et al., 2011c; Fabrycky et al., 2014). Almost 80% of the Kepler candidates are in single transiting systems. This has been referred to as the Kepler dichotomy (Johansen et al., 2012; Fang & Margot, 2012; Ballard & Johnson, 2016b; Moriarty & Ballard, 2016). Different scenarios have been proposed as an attempt to explain this dichotomy. On one hand, it has been suggested that this remarkable excess of single transiting planets is indeed real, that is, single transiting planets are truly singles (e.g., Johansen et al., 2012). On the other hand, it has been suggested that the dichotomy is only apparent and that the excess of single transiting planets arises from observing higher multiplicity systems where only one planet transits (e.g., Izidoro et al., 2017).

Synthetic transit observations of a mix of stable and unstable systems produced in Izidoro et al. (2017) suggest that the Kepler dichotomy arises from observing multi-planet systems with planets in mutually inclined orbits. Systems with low mutual inclination (mostly stable systems) contribute mostly to observed high-N planet systems and large mutual inclinations of unstable systems naturally produce a peak in the observed planet-multiplicity distribution at ${\rm N_{det}=1}$. Izidoro et al. (2017) matched the Kepler dichotomy assuming that Kepler planets comprise about $\lesssim 10\%$ of stable and $\gtrsim$90% of unstable systems. If this is correct, the Kepler dichotomy consists of a dichotomy in the inclination distribution rather than in planet multiplicity. Therefore, here, we also investigate how our simulations match the Kepler dichotomy and how our results compare to Izidoro et al. (2017).

Figure 17 shows the multiplicity distributions of our synthetically detected systems compared with that of the Kepler sample. Here we make use of the simulated detections produced earlier in this section. We note that more than 90% of the unstable planetary systems of Model I with ${\rm t_{start}=0.5\leavevmode\nobreak\ Myr}$ have two or more planets inside 0.7 AU (red dashed line in the left-hand panel of Figure 17). The two stable systems of Model I with ${\rm t_{start}=0.5\leavevmode\nobreak\ Myr}$ have six and eight planets in their chains.

Simulated detections of unstable systems of Model I with ${\rm t_{start}=0.5\leavevmode\nobreak\ Myr}$ provide a very good match to the Kepler multiplicity. Simulated observations of unstable systems rise the fraction of single planets from $\sim$10%555in reality, these unstable systems have additional planets but they are beyond 0.7 AU and thus do not account in our analysis. We also stress that these systems formed with multiple planets inside 0.7 AU but dynamical instabilities after gas dispersal resulted in multiple scattering events and collisions resulting in a single planet inside 0.7 AU. to about $\sim$75%. This result is consistent with that of Izidoro et al. (2017). Simulated observations of stable planetary systems result in a very spread and almost flat multiplicity distribution.

We show that our simulated observations match the period ratio distribution of the Kepler sample if one mixes typically $\lesssim 5\%$ of stable systems with $\gtrsim 95\%$ of unstable ones. Figure 17 shows that mixing 1% stable and 99% unstable systems provides an almost perfect match to the Kepler dichotomy in Model I with ${\rm t_{start}=0.5\leavevmode\nobreak\ Myr}$. The middle and right-hand panels of Figure 17 show that Model I with ${\rm t_{start}=3.0\leavevmode\nobreak\ Myr}$ and Model III with ${\rm t_{start}=0.5\leavevmode\nobreak\ Myr}$ do not match the Kepler dichotomy equally well, even assuming that only 1% of the systems are stable. The detected multiplicity distribution in the middle and right-hand panels of Figure 17 show a deficit at ${\rm N_{det}=1}$ (dashed green line) when compared to the Kepler population. Model II provides a poorer match to the observed multiplicity distribution compared to other models (see bottom-left panel of Figure 17). We also reiterate the fact that Model II fails to match Kepler observations period ratio distribution (Figure LABEL:fig:simulated_observations).

To better understand how planetary systems with different (true-)numbers of planets contribute to the origin of the observed dichotomy, we selected simulations of Model I to conduct a deeper analysis. Figure 18 shows again (as in Figure 17) the observed planet multiplicity distribution in synthetic transit observations of a mix of 1% stable and 99% unstable systems (dashed green histograms). The colors filling each ${\rm N_{det}}$ dashed green histogram bin show the fractional contributions as a function of true system multiplicity. In the left panel, where ${\rm t_{start}=0.5\leavevmode\nobreak\ Myr}$, 80% of all the observed systems show a single transit. However, only $\sim$5% of these systems have only one planet inside 0.7 AU, $\sim$39% of systems have two planets, $\sim$28% have three planets, and $\sim$8% have four or more planets. These results show that, in this case, the excess of single detected planets comes mainly from systems of intrinsically 2-3 super-Earths. In the right panel, where ${\rm t_{start}=3.0\leavevmode\nobreak\ Myr}$, the peak at ${\rm N_{det}=1}$ constitutes $\sim$65% of the systems, where $\sim$11% of the systems have two planets, $\sim$14% have three planets, $\sim$21% have four planets, and 18% of the systems have five or more planets. In this latter case, the relative contributions in terms of the systems’ true multiplicity are not dramatically different and are within a factor of two.

Although Model I with ${\rm t_{start}=3.0\leavevmode\nobreak\ Myr}$ and ${S_{{\rm peb}}=5}$ and Model III with ${\rm t_{start}=0.5\leavevmode\nobreak\ Myr}$ and ${S_{{\rm peb}}=10}$ do not perfectly match the dichotomy, there may be several false positives in the ${\rm N_{det}=1}$ planetary systems of Kepler. The rate of false positives for single Kepler planets is estimated to be of 20%, at least two times higher than the false-positive rate in multi-planet systems (Morton & Johnson, 2011; Fressin et al., 2013; Coughlin et al., 2014; Désert et al., 2015; Morton et al., 2016). Thus, the peak in the Kepler multiplicity distribution at ${\rm N_{det}=1}$ could be shorter in reality. Additionally, some single-planet systems in the Kepler sample may be truly singles formed by different mechanisms (Izidoro et al., 2015, 2017). Therefore, Model I with ${\rm t_{start}=3.0\leavevmode\nobreak\ Myr}$ and Model III with ${\rm t_{start}=0.5\leavevmode\nobreak\ Myr}$ might be good after all but with an extremely high instability fraction.

6.4 Super-Earths: Rocky or Icy?

As discussed above, the migration model has been shown to produce mostly ice-rich super-Earths (see Raymond et al., 2018a). This is in conflict with the results of atmospheric photoevaporation models, which suggest that the super-Earths closest to the parent star – and perhaps the majority of all super-Earths – are mostly rocky (Owen & Jackson, 2012; Lopez & Fortney, 2013; Owen & Wu, 2017; Lopez, 2017; Jin & Mordasini, 2018). The water-rich composition of super-Earths in our nominal simulations is probably more consistent with interior modeling of the composition of super-Earths, which suggests that super-Earths larger than 2-3${\rm\,R_{\oplus}}$ are water-rich worlds (Zeng et al., 2019).

In this section we first examine the compositions of close-in super-Earths from the simulations presented above. Then, inspired by the current debate on the true nature of super-Earths, we perform additional sets of simulations with different assumptions to find the conditions needed to form rocky super-Earths.

6.4.1 Compositions of super-Earths in Models I, II, and III

The super-Earths that grew in our Model I and Model II simulations (see Table LABEL:tab:2) grew mainly from pebbles originating from beyond the snow line. Thus, their final ice-mass fractions are as high as 50% (e.g., Figures 6, 7). Also, as suggested by Figure 8, only Model III simulations with large refractory pebbles (${\rm Rocky_{{\rm peb}}=1\leavevmode\nobreak\ cm}$) or extremely large pebble fluxes produced multiple inner planets with low water-ice contents.

Figure 19 shows the final planetary systems produced in 50 simulations of Model III with ${\rm t_{start}=0.5\leavevmode\nobreak\ Myr}$, ${S_{{\rm peb}}=10,}$ and ${\rm Rocky_{{\rm peb}}=0.1\leavevmode\nobreak\ cm}$. The results of these simulations are separated in two batches. The left-hand panel shows planetary systems that remained stable after the gas disk dispersal. The right-hand panel shows the final planetary systems that became dynamically unstable after the gas dispersal. It is clear that the former are much more compact than the latter for reasons explained above. More importantly, Figure 19 shows a diversity of planetary system compositions. The left-hand panel shows that in dynamically stable systems, the innermost planets are typically rocky. This is a consequence of resonant shepherding. As the snow line moves sufficiently inwards, it first sweeps the outermost planetary embryos in the disk; these objects grow faster and migrate inwards, and on their way inwards they encounter in resonance growing rocky planetary objects and shepherd them towards the inner edge of the disk. This result supports those of Raymond et al. (2018a) who proposed that the migration model is consistent with a variety of super-Earth compositions.

Some planetary systems in the left-hand panel of Figure 19 are particularly interesting, such as for example the second planetary system from the bottom. This planetary system shows five planets with masses ranging between 1.7 and ${\rm\sim 9\leavevmode\nobreak\ M_{\oplus}}$ with bulk compositions that alternate radially from rocky to icy. This shows that adjacent planets may have drastically different feeding zones. Even after dynamical instabilities (right-hand panel of Figure 19) the innermost planets in several systems are typically rocky. In some cases, rocky super-Earths end up on orbits that are immediately adjacent to ice-rich super-Earths that are reminiscent of the Kepler-36 system (Carter et al., 2012) and consistent with the simulations of Raymond et al. (2018a). It is also important to point out that in Model III with ${\rm Rocky_{{\rm peb}}=0.1\leavevmode\nobreak\ cm}$ the largest rocky close-in super-Earths in stable systems have typical masses of ${\rm 2M_{\oplus}}$ or lower. Observations of photoevaporated planets show rocky super-Earths up to at least 5${\rm M_{\oplus}}$. unstable systems of Model III with ${\rm Rocky_{{\rm peb}}=0.1\leavevmode\nobreak\ cm}$ do produce some ${\rm\sim 5M_{\oplus}}$ rocky super-Earths but we can produce them far more easily in Model III with ${\rm Rocky_{{\rm peb}}=1\leavevmode\nobreak\ cm}$.

We also note that in our unstable systems of Figure 19, planets farther out are typically icy, resulting in an overall population of mixed compositions with most planets being icy. These icy planets are beyond 0.1 AU, and as they are typically larger than ${\rm\sim 1-2\leavevmode\nobreak\ M_{\oplus}}$ they are very unlikely to lose their entire water content by photo-evaporation (Kurosaki et al., 2014). The co-existence of water-rich super-Earths and relatively less massive rocky super-Earths in the same system, in Figure 19, is at least qualitatively consistent with the findings of Zeng et al. (2019), which suggest that super-Earths larger than 2${\rm\,R_{\oplus}}$ are water-rich worlds. We note that Zeng et al. (2019) also suggest that the innermost (hottest) super-Earths (shown as redish dots in their Figure 1) are all rocky. This is roughly consistent with the results of Figure 19. Our results in Figure 19 are also consistent with a few planetary systems tested against the photoevaporation model (Owen & Campos Estrada, 2020, e.g., Kepler-100, Kepler-142, and Kepler-36) which suggest that adjacent planets in a system may have very distinct water contents.

We conclude this section by emphasizing that the absence of icy super-Earth systems in our companion paper (Paper I; Lambrechts et al., 2019) is consequence of assuming an initial distribution of seeds only inside the snow line and without taking into account the snow line’s evolution. The results of Paper I and this work agree on the fact that the formation of systems of rocky super-Earths requires special conditions. In the following section, we test if a different disk model or a putative lack of gas-driven planet migration can result in the formation of rocky super-Earths.

6.4.2 How can we produce rocky super-Earths?

In this section, we attempt to decipher whether or not we can produce close-in planetary systems dominated by rocky super-Earths by ignoring type-I migration, considering a nonevolving gas disk, or a combination of these two scenarios.

We performed three additional sets of simulations to answer this question. In all of these, seeds are initially distributed from 0.5 AU to $\sim$20 AU. The initial mutual radial separation of seeds are set by a geometric progression with common ratio equal to 1.06. Other orbital elements of our initial seed distribution were sampled as described in Section 2.2. We name this scenario Model IV. We stop the numerical integration of all simulations of this section at the end of the gas disk phase.

Figure 20 shows the final masses of planets growing by pebble accretion in a nonevolving disk (both gas surface density and temperature are kept constant) with different pebble fluxes. The gas disk structure is set by the starting time of the simulation (${\rm t_{start}=0.5\leavevmode\nobreak\ Myr}$). We keep the initial disk structure throughout the entire simulation, and therefore the snow line is kept fixed at about ${\rm\sim 3\leavevmode\nobreak\ AU}$. This allows us to test the effects of the movement of the snow line on our results. We recall that, in our nominal simulations, at the end of the gas disk phase the snow line has shifted into the inner solar system down to ${\rm\sim 1\leavevmode\nobreak\ AU}$.

Figure 20 shows that this scenario also fails to produce rocky super-Earths for all considered pebble fluxes. Close-in Earth-mass planets are predominately icy. Seeds from the outer disk grow faster and regulate the pebble flux to the inner region, impairing the growth of rocky seeds. Also, icy super-Earths eventually migrate into the inner disk, scattering or colliding with small rocky seeds.

Figure 21 shows the final masses of planets in simulations where type-I migration is neglected but the disk evolves as in our nominal simulations. The results of simulations with different pebble fluxes are shown. Unlike the previous scenario, seeds from the outer disk do not enter into the inner system but they still grow quickly, starving the seeds in the inner disk. Our results show that rocky seeds barely grow from their initial size because there is too much filtering of pebbles from outer, growing planets and pebbles are too small in the inner system, making accretion of rocky seeds very inefficient. For completeness, we also performed simulations where type-I migration is neglected and the gas disk is not evolving. The results of this scenario are qualitatively equivalent to those of Figure 21.

We conclude this section by stressing that we did not succeed in producing rocky super-Earth systems in any of our simulations where the initial distribution of seeds extends up to distances reasonably beyond the snow line (e.g., ${\rm\sim 10\leavevmode\nobreak\ AU}$). If super-Earths are mostly rocky, our best match to this constraint —while still far from ideal— comes from the results of the simulations of Model III with ${\rm Rocky_{{\rm peb}}=0.1\leavevmode\nobreak\ cm}$ and ${S_{{\rm peb}}=10}$. In the following section we revisit Model III, invoking cm-sized silicate pebbles in the inner disk as a path for more efficient planetary growth in the inner disk.

6.4.3 Making rocky super-Earths systems: A more successful case

Our goal for this section is to modify the parameters of our simulations in order to produce predominantly rocky super-Earths as well as a high fraction of unstable systems. Figure 14 shows that dynamical instabilities after gas dispersal are more likely in systems that are very dynamically compact at the end of the gas disk phase and more importantly with a sufficiently larger number of planets anchored at the disk inner edge (see also Matsumoto et al. (2012)). The simulations of Figure 9 were the only ones that produced predominantly rocky super-Earths inside of 0.7 AU. As a follow-up to that, we perform an additional set of simulations of Model III considering cm-sized silicate pebbles in the inner regions (${\rm Rocky_{{\rm peb}}=1\leavevmode\nobreak\ cm}$) and ${S_{{\rm peb}}=5}$. We do not re-use the results of Figure 9 because there the disk inner edge was far from the star (at about 0.3 AU). Therefore, in order to compare our results with observations and for consistency with our previous analyses, we redo these simulations considering $r_{{\rm in}}=0.1\leavevmode\nobreak\ {\rm AU}$. We note that we stop the simulations of Figure 9 at 5 Myr. We performed 50 new simulations. Because in this setup, most planets are fully formed during the first 2 Myr or less, and the disk lifetime in these simulations is set equal to 2.5 Myr, instead of the 5 Myr used for the simulations presented in Figure 9 (this also allows us to save CPU time because these simulations are very computationally expensive). We follow the long-term dynamical evolution of all these systems up to about 300 Myr after gas dispersal. Figure 22 shows our results.

The upper-panels of Figure 22 show planet mass, planet orbital inclination, and the last collision epoch distributions of stable and unstable systems. As expected, stable systems have lower mass planets and less mutually inclined planet pairs compared to unstable systems. The mass distribution of planets in unstable systems agrees quite well with the inferred masses of Kepler planets from mass–radius relationship models (Wolfgang et al., 2016). Moreover, the left-hand panel of Figure 22 shows that 88% of these systems became dynamically unstable after gas dispersal. This is different from simulations of Model III with ${\rm Rocky_{{\rm peb}}=0.1\leavevmode\nobreak\ cm}$ where only ${\sim 50\%}$ of the systems became unstable666We also integrated the stable systems of Model III with ${\rm Rocky_{{\rm peb}}=0.1\leavevmode\nobreak\ cm}$ up to 300 Myr but the fraction of unstable systems did not increase significantly, remaining around 50%.. This results in a much better though still imperfect match to observations (see bottom-right panel of 22).

The bottom-left panel of Figure 22 shows the period ratio distribution of planet pairs of stable and unstable systems and also their respective simulated detections. The bottom-middle panel shows statistics from a KS-test comparing Kepler observations with simulated detections of samples mixing different real fractions of stable and unstable systems. The bottom-left panel shows the planet multiplicity distributions from stable and unstable systems and also from our simulated detections. Overall, the results presented in this set of panels are very similar to those of Model III with ${\rm t_{start}=0.5\leavevmode\nobreak\ Myr}$, ${\rm Rocky_{{\rm peb}}=0.1\leavevmode\nobreak\ cm}$, and ${S_{{\rm peb}}=10}$ (see bottom panels of Figure LABEL:fig:simulated_observations and left-panel of 17). The period ratio distribution of our simulated planet pairs matches observations when one mixes about $<$2% stable planet pairs with $>98$% unstable planet pairs. Our simulated detections do not perfectly match the Kepler dichotomy but there is nevertheless a prominent peak at ${\rm N_{det}=1}$ qualitatively similar to that seen in observations. The fraction of systems with single detected planets in our simulated detections is about 60% compared to 80% of observations (but see section 6.3 for a discussion about the rate of false positives among single transiting planets in the Kepler data). Curiously, we also find that these simulations can provide a better match to the Kepler dichotomy if we rescale outwards the semi-major axis of planets in our simulations by a factor of 2-3. This is because close-in planets are more likely to be detected than those farther out. In our simulations, the disk inner edge is set at ${\rm\sim 0.1\leavevmode\nobreak\ AU}$ and this essentially sets the typical location of the innermost planets in our simulations. If we had considered the disk inner edge to be slightly further out, our simulations would probably better match observations. This remains an interesting issue for future studies.

Finally, Figure 23 shows all unstable systems produced in simulations of Model III with ${\rm Rocky_{{\rm peb}}=1\leavevmode\nobreak\ cm}$ at the very end of our simulations (300 Myr). Most planets produced in these simulations are rocky instead of icy.

These simulations confirm the expectations from Figure 9 that to form rocky super-Earths and achieve a final high instability fraction, growth inside the snow line has to be fast. Only in this case can planets migrate faster than the snow line and avoid accreting icy pebbles at all times. To achieve such fast growth, we invoked large silicate pebbles in the inner disk. However, a similar growth mode could be equally achieved by invoking less stirring of the silicate pebble layer with mm-sized pebbles. Nevertheless, we stress that any of these scenarios can only succeed in forming rocky super-Earths in our model if the initial distribution of seeds is restricted to regions well inside the snow line.

7 The Solar System in context

The Solar System is unusual among known exoplanet systems (see recent review by Raymond et al., 2018b). It has been estimated that only $\sim$8% of the planetary systems have the innermost planet (${\rm>1\leavevmode\nobreak\ R_{\oplus}}$) at an orbital period longer than that of Mercury (Mulders et al., 2018). The structure of the Solar System is also clearly segregated with low-mass rocky terrestrial planets residing in the inner regions and gaseous and/or icy giant planets in the outer region. The origin of this dichotomy in mass has been interpreted as a consequence of the process of pebble accretion and more precisely the presence of small silicate pebbles in the inner Solar System and larger icy pebbles in the outer Solar System (Morbidelli et al., 2015).

Jupiter and Saturn are the great architects of the Solar System. The cores of Jupiter and Saturn probably formed early and regulated the pebble flux to the terrestrial region. They first intercepted and consumed part of the pebbles drifting inwards but eventually reached pebble isolation disconnecting the inner and outer Solar System (Kruijer et al., 2017; Lambrechts et al., 2019; Weber et al., 2018). As a consequence of this process, terrestrial protoplanetary embryos were starved and only grew to about Mars-mass —at most— and did not migrate to the inner edge of the disk to become close-in super-Earths (see Paper I and Morbidelli et al. (2015)). As Jupiter and Saturn grew, they probably migrated into a resonant configuration —also because of the interaction with sun’s natal disk (e.g., Ward, 1986; Lin & Papaloizou, 1986). The exact gas-driven migration history of Jupiter and Saturn is not constrained, but their specific mass ratio prevented their migration into the Earth’s zone (Masset & Snellgrove, 2001; Morbidelli & Crida, 2007; Pierens & Raymond, 2011; Pierens et al., 2014). Jupiter and Saturn may have also blocked the inward gas-driven migration of Uranus and Neptune (and their precursors) to the inner Solar System, thus preventing the migration of icy super-Earths into the inner system (Izidoro et al., 2015, 2015).

In light of our current view of the formation mechanism of the Solar System, none of our simulations comes close to matching our own system. Perhaps our closest approximation is produced in Model III with ${\rm Rocky_{{\rm peb}}=0.1\leavevmode\nobreak\ cm}$ and ${S_{{\rm peb}}=1}$ (see top-left of Figure 8). Figure 8 suggests that if just a few seeds form near the snow line, they can grow much more than the inner planetary seeds. In some of these simulations, the more massive cores have a few Earth-masses and sit around 2-3 AU while the innermost planetary embryos all have masses below that of Mars. In reality, we need the two innermost cores to grow more than in this simulation in order to become gas-giant planets before migrating too much in the type-I regime (once giant planets, their mutual interactions can prevent their type-II migration as in Masset & Snellgrove (2001); but see also Paper III). We can see some hope of this taking place in the narrow corner of parameter space explored in this paper but we are far from coming to a firm conclusion. Certainly, this issue requires further investigation.

8 Conclusions

We used N-body numerical simulations to model the growth and migration of planetary embryos in gaseous protoplanetary disks simultaneously. Our simulations start with a distribution of roughly subMoon-mass planetary seeds that grow predominantly by pebble accretion. Our results show that the integrated pebble flux primarily sets the final planet masses. Planetary embryos growing in low-pebble-flux environments —as in simulations where the integrated pebble flux is $\sim 39\leavevmode\nobreak\ M_{\oplus}$— have final typical masses of $\sim 2\leavevmode\nobreak\ M_{\oplus}$ or lower. Simulations where the total pebble reservoir is of $\sim 194\leavevmode\nobreak\ M_{\oplus}$ produce multiple super-Earth-mass planets. Pebble accretion stops when planetary embryos reach pebble isolations mass. As the disk evolves, radial migration promotes mutual collisions and the delivery of multiple (super) Earth-mass planets to the disk inner edge. Our simulations also show that a simple increase by a factor of two in the total pebble flux from $\sim 194\leavevmode\nobreak\ M_{\oplus}$ to about $\sim 485\leavevmode\nobreak\ M_{\oplus}$ is enough to bifurcate the evolution of our planetary systems from one leading to typical super-Earth-mass planets ($\lesssim 15\leavevmode\nobreak\ M_{\oplus}$) to another achieving systems of super-massive planetary embryos which are very likely to become gas giants. We dedicated two companion papers to modeling the formation of truly terrestrial planets as opposed to rocky super-Earths (Lambrechts et al., 2019) and gas giant planets (Bitsch et al., 2019). The focus of this paper is to model the formation and long-term dynamical evolution of close-in super-Earths systems. Although some observed super-Earth systems host external gas giant planets, it is not clear whether this is predominant among such systems (Barbato et al., 2018; Zhu & Wu, 2018; Bryan et al., 2019). In this work, we do not model the effects of very massive external perturbers on our close-in systems. This issue was recently addressed by Bitsch et al. (2020).

In our simulations we tested the effects of considering different silicate pebble sizes and also different initial radial distributions of seeds. Some of our models are more successful than others in matching observations. In some models, up to $\sim 95\%$ of the resonant chains become dynamically unstable after gas dispersal. This fraction of naturally unstable systems after gas dispersal is significantly higher than that found by Izidoro et al. (2017), providing a better match to observations. Overall, our simulations match the period ratio distribution of the Kepler sample if one combines a fraction of stable and unstable planetary systems, typically $\lesssim$2% stable with $\gtrsim$98% unstable. Supporting the results of Izidoro et al. (2017), the results of our simulations also suggest that the excess of detected single-planet transiting systems compared to multiplanet transiting systems arises from a dichotomy in orbital inclination rather than in planet multiplicity.

The Kepler period distribution has been also used to constrain formation models. Lee & Chiang (2017) and Choksi & Chiang (2020) argue that super-Earths do not undergo significant orbital migration during their last doubling in mass which is consistent with our results. In the context of the migration model, the period distribution of planets with orbital periods shorter than 100 days has also been found to be relatively consistent with gas-driven migration and dynamical instabilities (Carrera et al., 2019). The migration model and dynamical instabilities can also account for the origin of nearly coplanar, high-planet-multiplicity, and nonresonant systems of close-in super-Earths, such as the Kepler-11 system (Esteves et al. (2020); see also McNally et al. (2019)). Future studies —ideally including interior modeling and effects of gas accretion on planets, stellar tides, fragmentation, and atmospheric loss generated by impacts— should also aim to test how the migration model matches the inferred ‘peas in a pod pattern” found in Kepler’s high-multiplicity systems, which suggests that planets in a given system have very similar sizes (e.g., Weiss & Petigura, 2020).

Our simulations where planetary seeds are initially distributed in the inner and outer disk (inside and outside the snow line) produce systematically close-in icy super-Earths. This is in contrasts with the results of atmospheric photoevaporation models that suggest that super-Earths are mostly rocky (Owen et al., 2012; Jin & Mordasini, 2018; Van Eylen et al., 2018). Our finding is probably more aligned with the results of Zeng et al. (2019), which suggest that super-Earths larger than 2${\rm\,R_{\oplus}}$ are water-rich planets. Finally, we also show that the formation of close-in systems dominated by rocky super-Earths requires special conditions such as the formation of planetary seeds only inside the snow line (e.g., $<$1-2 AU) and vigorous growth by pebble accretion in the inner disk (suggesting the existence of large rocky pebbles or of a low-scale-height pebble layer in the inner disk). These same systems typically show that the innermost rocky super-Earths tend to be rocky whereas the outermost ones tend to be icy. This is consistent with the results of Zeng et al. (2019) which suggest that the super-Earths exposed to the greatest stellar fluxes are mostly rocky.

Overall, the results presented here are in agreement with Izidoro et al. (2017) and suggest that pebble accretion, migration, and dynamical instabilities is a powerful combination of mechanisms to explain the bulk of the orbital and mass properties of super-Earth systems. However, if super-Earths are all rocky, very special conditions would be required to promote their formation in our models. The required conditions seem to diverge from our current understanding of planetesimal formation, which indicates the snow line as the most favorable place to produce planetary seeds (Armitage et al., 2016; Dra̧żkowska & Alibert, 2017), and are also in contrast with the assumption that larger icy —rather than rocky— pebbles are necessary to explain the structure of the Solar System. Instead, our results suggest that at least some super-Earths should be water-rich planets.

Acknowledgements

A. I. thanks FAPESP for support via grants 16/19556-7 and 16/12686-2, and CNPq via process 313998/2018-3. Most simulations of this work were performed on the computer cluster at UNESP/FEG acquired with resources from FAPESP. S. N. R and A. M. thank the Agence Nationale pour la Recherche for support via grant ANR-13-BS05-0003- 01 (project MOJO). B.B., thanks the European Research Council (ERC Starting Grant 757448-PAMDORA) for their financial support. A.J. was supported by the European Research Council under ERC Consolidator Grant agreement 724687-PLANETESYS, the Swedish Research Council (grant 2014-5775), and the Knut and Alice Wallenberg Foundation (grants 2012.0150, 2014.0017, and 2014.0048). Computer time for this study was also 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. This research has made use of the NASA Exoplanet Archive, which is operated by the California Institute of Technology, under contract with the National Aeronautics and Space Administration under the Exoplanet Exploration Program.

Appendix A Planet migration prescription

Growing planets interact gravitationally with the protoplanetary gas disk, which exerts a torque on the planet, and can thus migrate through the disk (see Baruteau et al. 2014 for a review). Low-mass planets migrate by type-I migration, whereas massive planets migrate by type-II migration. The torques responsible for type-I migration are the Lindblad torques, exerted from the spiral waves of the planet (Ward, 1986, 1997), and the corotation torques, originating from the horseshoe region around the planet (Goldreich & Tremaine, 1979; Ward, 1992). The Lindblad torque results in inward migration for most disk parameters, where the migration timescale is much shorter than the disk lifetime for bodies more massive than the Earth (Tanaka et al., 2002).

However, the corotation torque can, depending on the local disk properties, overcompensate the negative Lindblad torque, resulting in outward migration (Paardekooper & Mellema, 2006b; Paardekooper & Mellema, 2008; Baruteau & Masset, 2008; Paardekooper & Papaloizou, 2008; Kley et al., 2009; Paardekooper et al., 2011). This outward migration depends strongly on the radial profile of the gas surface density and temperature and thus entropy in the disk (Baruteau & Masset, 2008; Bitsch & Kley, 2011b). So-called regions of outward migration, where the total torque is positive and can be sustained by the entropy-driven corotation torque in the disk, are associated with shadowed regions, where the disk’s aspect ratio H/r decreases with orbital distance (Bitsch et al., 2013b, 2014, 2015; Baillié et al., 2015). Additionally, outward migration can exist if the radial gas surface density gradient increases with orbital distance (Masset et al., 2006), as at the inner edge of the disk.

For the migration of planets, we follow the torque formula given by Paardekooper et al. (2011), which has been extensively tested against 3D hydrodynamical simulations (Bitsch & Kley, 2011b; Lega et al., 2014, 2015).

In the torque formula by Paardekooper et al. (2011), the total torque $\Gamma_{\rm tot}$ is the sum of the Lindblad torque $\Gamma_{\rm L}$ and the corotation torque $\Gamma_{\rm C}$. The total type-I torque also depends on the planet’s orbital eccentricity and inclination. In order to account for the eccentricity and inclination effects, the total torque formula of Paardekooper et al. (2011) is rewritten as

[TABLE]

where $\Delta_{\text{L}}$ and $\Delta_{\text{C}}$ are rescaling factors accounting for torque reduction due to the orbital eccentricity and inclination of the planet. The Lindblad torque is

[TABLE]

where $\beta$ and x are the negative of the gas surface density and temperature gradients at the location of the planet (e.g., Lyra et al., 2010):

[TABLE]

The Lindblad torque reduction $\Delta_{\text{L}}$ (Cresswell & Nelson, 2008) is

[TABLE]

where

[TABLE]

where e and i are the planet’s orbital eccentricity and inclination, respectively, and h is the gas disk aspect ratio.

The co-orbital torque is written as

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

In the previous equations, ${\xi=\beta-(\gamma-1)x}$ is the negative of the entropy gradient and $\gamma=1.4$ is the adiabatic index. The scaling torque ${\Gamma_{0}=(q/h)^{2}\Sigma_{\text{gas}}r^{4}\Omega_{k}^{2}}$ is calculated at the planet’s location. Here, q is the planet–star mass ratio, and ${\Omega_{k}}$ is the planet’s Keplerian orbital frequency. The functions F, G, and K control the torque saturation due to viscous and thermal effects. They depend on

[TABLE]

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

[TABLE]

and

[TABLE]

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

[TABLE]

In Eq. 33, $\rho$ is the gas volume density, $\kappa$ is the gas disk opacity, and $\sigma$ is the Stefan-Boltzmann constant. To set the protoplanetary disk opacity we follow Bell & Lin (1994).

Different diffusion timescales (e.g., viscosities) influence the corotation torque as well, where lower levels of viscosity will not allow outward migration even when the radial gradients in entropy are steep enough to promote outward migration (e.g., Baruteau & Masset, 2008; Bitsch et al., 2013a). Different viscosities can therefore result in different migration speeds, which might influence the structure of planetary systems in the inner disk (e.g., Bitsch et al., 2013a).

$\gamma_{\text{eff}}$ depends on

[TABLE]

as

[TABLE]

F, G, and K take the form

[TABLE]

[TABLE]

and

[TABLE]

Here, p is either $p_{\nu}$ or $p_{\chi}$ as defined above.

The co-orbital torque reduction is

[TABLE]

where ${e_{\text{f}}}$ is given by Fendyke & Nelson (2014) as

[TABLE]

We nte that for sufficiently large eccentricities, $\Delta_{\text{L}}$ may assume negative values favoring unrealistic reversal of migration. In extreme situations, where $\Delta_{\text{L}}<0,$ we impose $\Delta_{\text{L}}$ = $\Delta_{\text{C}}$. This approach is slightly different from that in Izidoro et al. (2017), but both approaches produce qualitatively equivalent results.

We follow Papaloizou & Larwood (2000) and Cresswell & Nelson (2008) to write the migration timescale used in our N-body code:

[TABLE]

where L and ${\rm\Gamma_{tot}}$ are the planet’s orbital angular momentum and the total type-I torque, respectively.

The damping of orbital eccentricity and inclination due to gas tidal interaction follow the prescriptions of Papaloizou & Larwood (2000) and Tanaka & Ward (2004) reformulated by Cresswell & Nelson (2006, 2008). ${\rm t_{e}}$ and ${\rm t_{i}}$ are the orbital eccentricity and inclination damping timescales, respectively:

[TABLE]

[TABLE]

In Eq. 42 and 43,

[TABLE]

with ${\rm M_{\odot}}$, ${\rm a_{p}}$, ${\rm m_{p}}$, ${\rm i}$, and ${\rm\leavevmode\nobreak\ e}$ being the solar mass, the planet’s semimajor axis, its mass, its orbital inclination, and its eccentricity, respectively.

Finally, in order to incorporate all these effects into our simulations, we follow Papaloizou & Larwood (2000) and Cresswell & Nelson (2006, 2008) and implement the following artificial accelerations into FLINTSTONE accounting for migration ($\mathbf{a}_{m}$), orbital eccentricity ($\mathbf{a}_{e}$), and inclination damping ($\mathbf{a}_{i}$):

[TABLE]

[TABLE]

and

[TABLE]

where ${\mathbf{r}}$, ${\mathbf{v}}$, and ${\mathbf{k}}$ are the planet’s heliocentric distance, velocity, and the unit vector in the z-direction. The formula of $\mathbf{a}_{i}$ includes a factor of two correcting for a typo in (Cresswell & Nelson, 2006, 2008).