© ESO 2021

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 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, 2015; Boley & Ford 2013; Boley et al. 2014; 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, 2018; Liu et al. 2015, 2016; Izidoro et al. 2017; Ormel et al. 2017; Unterborn 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); Goldreich & Schlichting (2014); Deck & Batygin (2015); Xu et al. (2018); Pichierri et al. (2018); 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. 2011b; 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 ≲10% of stable and ≳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 thelegitimacy 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, 2017; Simon et al. 2016). 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. 2015a,b; Bitsch et al. 2015b; 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 ofterrestrial 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 M_⊕ – lead to the slow growth of protoplanetary embryos that do not migrate substantially during the gas disk lifetime. These embryos aretypically 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− 480 M_⊕) 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. 2015a; 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 Sect. 4 we discuss our results in light of those of Paper I. In Sect. 5 we present the results of our simulations dedicated to modeling the growth and long-term dynamical evolution of hot super-Earth systems. In Sect. 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 Sect. 7 we place the Solar System in the context of our model. In Sect. 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 usethe 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. 2015a)¹. 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 $${\dot{M}}_{\text{gas}}=3\pi \alpha h^2{r}^2{\Omega }_{\text{k}}{\Sigma }_{\text{gas}},$$(1)

where h = H∕r is the disk aspect ratio with H being the disk scale height, α the dimensionless viscosity parameter (Shakura & Sunyaev 1973), r the distance to the star, ${\Omega }_{\text{k}}=\sqrt{G{M}_{\odot }/r^3}$ the orbital the Keplerian frequency, and Σ_gas the gas surface density.

Following Hartmann et al. (1998) and Bitsch et al. (2015a), the relationship between the disk and star age and the gas accretion rate onto the star is given by $$\text{log}\left(\frac{{\dot{M}}_{\text{gas}}}{M_{\odot }{\text{yr}}^{-1}}\right) =-8-1.4\log \left(\frac{t_{\text{disk}}+{10}^5 \text{yr}}{{10}^6 \text{yr}}\right).$$(2)

Finally, the hydrostatic equilibrium yields $$T=h^2\frac{G{M}_{\odot }}{r}\frac{\text{upmu}}{\mathcal{R}},$$(3)

where t_isk is the disk age, T is the temperature at the gas disk midplane, μ is the gas mean molecular weight set to 2.3 g mol⁻¹, G is the gravitational constant, M_⊙ is the mass of the star set equal to 1 solar mass, and $\mathcal{R}$ is the ideal gasconstant. We use t_start to represent the disk age at the starting time of our simulations.

From Eqs. (1)–(3) we obtain the gas surface density (Σ_gas) of the disk by using the fits of the disk temperature at the midplane provided in Bitsch et al. (2015a). The disk metallicity and α-viscosity parameter are set to 1% and to α = 5.4 × 10⁻³, 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 yr 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 (≫500 yr).

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 ~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 ~ 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 $$\Re =\tan \text{h}\left(\frac{r-r_{\text{in}}}{0.05r_{\text{in}}}\right),$$(4)

where r is the heliocentric distance and r_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⁻⁴. 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, 2018; Armitage et al. 2016; Dra̧żkowska & Alibert 2017; Carrera et al. 2017; Schoonenberg & Ormel 2017). Other models suggest planetesimal formation takes place just inside the snow line (Ida & Guillot 2016) or even around 1 AU (Drąż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 1. In simulations of Models I and II, planetary seeds are initially distributed past 0.7 AU. Because at t_disk = 3 Myr the snow line is at 0.7 AU, simulations with t_start = 3 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-memberscenarios. 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 t_start = 0.5 Myr. In the second scenario, planetesimals are assumed to form late and our simulations start with t_start = 3.0 Myr (see also Paper III). We note that different t_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. 2015b).

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 and 2 g cm⁻³, 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. 2015a). 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 $${\dot{M}}_{\text{core}}=\pi R_{\text{acc}}^2{\rho }_{\text{p},\text{mid}}\overline{S}\delta v,$$(5)

where ρ_p,mid is the midplane density of pebbles, related to the pebble surface density layer Σ_peb via $${\rho }_{\text{p},\text{mid}}={\Sigma }_{\text{peb}}/(\sqrt{2\pi }{H}_{\text{peb}}),$$(6)

where $H_{\text{peb}}=H_{\text{gas}}\sqrt{{\alpha }_{\text{set}}/{\tau }_{\text{f}}}$ is the pebble scale height (Youdin & Lithwick 2007), τ_f is the Stokesnumber, and α_set is the dimensionless viscosity parameter representing disk midplane turbulence which determines the vertical settling of pebbles. In our nominal simulations, α_set = α, although some studies argue for α_set ≪ α (Zhu & Stone 2014). R_acc denotes the accretion radius of theplanet, δv = v_rel + ΩR_acc with v_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 $\overline{S}$ is defined as the mean pebble density normalized by the pebble density in the midplane, namely: $$\overline{S}=\frac{1}{\pi R_{\text{acc}}^2}{\displaystyle {\int }_{z_0-R_{\text{acc}}}^{z_0+R_{\text{acc}}}2}\exp [-z^2/(2H_{\text{peb}}^2)]\sqrt{R_{\text{acc}}^2-{(z-z_0)}^2}\text{dz},$$(7)

where z₀ 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 $${\dot{M}}_{\text{peb}}=2\pi r_{\text{g}}\frac{\text{d}{r}_{\text{g}}}{\text{d}t}Z{\Sigma }_{\text{g}}(r_{\text{g}}).$$(8)

Here, r_g represents the heliocentric distance at which dust particles have grown to pebble sizes and start drifting inwards by gas drag, and Σ_g (r_g) is the gas surface density at the pebble production line location r_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_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 $$r_{\text{g}}={\left(\frac{3}{16}\right)}^{1/3}{(G{M}_{\star })}^{1/3}{({ϵ}_{\text{D}}Z)}^{2/3}{t}^{2/3},$$(9)

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

where M_⋆ is the stellar mass, which we set to 1M_⊙, G is the gravitational constant, and ϵ_D = 0.05 is associatedwith 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 ~77 AU, ~ 250 AU, and ~ 350 AU, respectively². 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 ~ 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_peb = 1, the total mass in pebbles between 50 and 90 AU is only about 10 M_⊕ at 1 Myr. Observations suggest that the total mass in dust and pebbles in the outer parts (50–90 au) of the HL Tau disk, which is about 1 Myr old, is between 80 and 280  M_⊕ (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_peb in order to increase the pebble surface density (or decrease in some extreme cases; see Eq. (13)). We note that S_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 Sect. 2.4 in Lambrechts & Johansen 2014) the stokes number may be written as $${\tau }_{\text{f}}=\frac{\sqrt{3}}{8}\frac{{ϵ}_{\text{P}}}{\eta }\frac{{\Sigma }_{\text{peb}}(r_{\text{P}})}{{\Sigma }_{\text{g}}(r_{\text{P}})},$$(11)

where Σ_g(r_P) and Σ_peb (r_P) are the gas and pebble surface densities at the location of the planets r_P. We represent the radial pressure support of the disk through the dimensionless parameter $$\eta =-\frac{1}{2}{\left(\frac{H}{r}\right)}^2\frac{\partial \ln P}{\partial \ln r}.$$(12)

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 Σ_peb can be calculated using the pebble flux and Stokes number (Eqs. (8) and (11)) as $${\Sigma }_{\text{peb}}=\sqrt{\frac{2S_{\text{peb}}{\dot{M}}_{\text{peb}}{\Sigma }_{\text{g}}}{\sqrt{3}\pi {ϵ}_{\text{P}}{r}_{\text{P}}{v}_{\text{K}}}},$$(13)

where r_P denotes the semi-major axis of the planet, and v_K = Ω_Kr_P. S_peb is a nondimensional linear scaling factor of the pebble flux (see Eq. (8)). In Eq. (13), ϵ_P = 0.5 represents the pebble sticking efficiency under the assumption of near-perfect sticking (Lambrechts & Johansen 2014). S_peb = 1 corresponds to an integrated pebble flux I_peb ≈ 194 M_⊕ beyond the snow line (see Fig. 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 $${\tau }_{\text{f},\text{rock}}=\frac{{\rho }_{\text{rock}}{\text{Rocky}}_{\text{peb}}}{{\Sigma }_{\text{g}}},$$(14)

(e.g., Lambrechts & Johansen 2012) where Rocky_peb and ρ_rock define the size and bulk density of silicate pebbles, respectively. The bulk density of silicate pebbles is assumed to be 5.5 g cm⁻³ 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 (Ṁ_peb,rock = 0.5Ṁ_peb). We calculate the pebble surface density inside the snow line using the following equation. $${\Sigma }_{\text{peb},\text{rock}}=\frac{{\dot{M}}_{\text{peb},\text{rock}}}{2\pi r{v}_{\text{peb},\text{rock}}},$$(15)

where the radial velocity of rocky pebbles is given by (e.g., Ormel & Klahr 2010; Bitsch et al. 2018a) $$v_{r_{\text{peb},\text{rock}}}=-\frac{2.0\eta v_{\text{k}}}{{\tau }_{\text{f},\text{rock}}+{\tau }_{\text{f},\text{rock}}^{-1}}-\frac{\nu }{r},$$(16)

where ν = αH²Ω_k is the gas kinematic viscosity.

In our nominal simulations, we set Rocky_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 pebblesinside the snow line. We have verified that in our disk model, at the very early stages of the disk (t_disk ≈ 0.0 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 t_start = t_disk = 0.5 Myr or 3 Myr, both 1mm and 1cm silicate pebbles beyond ~0.2 AU (the starting location of our innermost seeds) are in the Epstein regime of gas-particle coupling. The top panel of Fig. 2 shows the gas and pebble surface densities as a function of orbital distance as the disk evolves. The bottom panel of Fig. 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 massfluxes 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 planetonly accretes a fraction f_acc of the flux of pebbles Ṁ_peb drifting pastits orbit: $$f_{\text{acc}}=\frac{{\dot{M}}_{\text{core}}}{{\dot{M}}_{\text{peb}}}.$$(17)

The pebble flux arriving at interior planets is thus reduced by this fraction f_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 $$M_{\text{iso}}=25f_{\text{fit}}{M}_{\text{E}}+\frac{{\Pi }_{\text{crit}}}{\lambda }{M}_{\text{E}},$$(18)

with λ ≈ 0.00476∕f_fit, ${\Pi }_{\text{crit}}=\frac{\alpha }{2{\tau }_{\text{f}}}$, and $$f_{\text{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]\text{},$$(19)

where α₃ = 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 slow³ 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 M_⊕ and in disks with opacities lower than ~ 0.01 cm² 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 M_⊕ during the gas disk phase.

Fig. 1 Pebble flux in our nominal disk model. The integrated mass in drifted pebbles (Ipeb) is shown on the left-hand vertical axis as a function of disk age (tdisk). Fluxes are calculated at 5 AU. The right-hand vertical axis shows the instantaneous pebble flux decreasing as the disk dissipates. Simulations of Paper I and III also feature the decay of the pebble flux. All our simulations start with a tdisk of at least 0.5 Myr, thus the pebble production line is already past the initial positions of our outermost seeds (~ 60 AU; see Table 1). This plot is generated considering Speb = 1 in Eq. (13) but a simple rescaling of these curves accounts for other considered pebble fluxes. Both curves correspond to the pebble flux beyond the snow line. The pebble flux inside the snow line is reduced by a factor of two because of the mass sublimation of the ice component of the pebbles when they cross the snow line. The gas flux and the integrated gas flux in our gas disk model are shown by the green and blue lines, respectively. We note that in the plot we rescaled the gas flux by a factor of 1/30, for clear comparison with the pebble flux.

Fig. 2 Top: gas and pebble surface densities as a function of orbital distance at different disk ages. Bottom: pebble sizes and Stokes number as a function of orbital distance at different disk ages. In both panels, Rockypeb = 1 mm and Speb = 1, which correspond to our nominal values. The sharp drop in the Stokes number in the inner regions of the disk is due to the snow line transition where we set the size of silicate pebbles to 1 mm.

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_peb = 0.2, 0.4, 1, 2.5, 5, and 10; see Table 1). 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 ~ 15 M_⊕ – during the gas disk phase. The second phase of long-term integrations beyond the gas disk phase is discussed in Sect. 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 114 M_⊕ leads to the formation of classical terrestrial planets while simulations with integrated pebble fluxes of 190 M_⊕ (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_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 r_in ≃ 0.3 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 r_in ≃ 0.1 AU in our simulations of Sect. 5. However, here we set r_in ≃ 0.3 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.

Fig. 3 Evolution of an example simulation from model I. The panels represent snapshots of the growth of protoplanetary embryos in a simulation with tstart = 0.5 Myr, Speb = 1, and a size of the rocky pebbles equal to Rocky_peb = 0.1 cm. The vertical and horizontal axes represent mass and semi-major axis, respectively. The location of the disk snow line (blue vertical dashed line) and the pebble isolation mass (black dashed line) are also shown for reference as the disk evolves. Planetary embryos growing inside the snow line accrete silicate pebbles while those outside the snow line accrete icy pebbles. The gray solid lines delimit the region of outward migration. The disk inner edge is shown at ~ 0.3 AU, where planets are trapped as well. The color of each dot gives the ice mass fraction. The size of each dot scales as m1∕3, where m is the planetary mass. The exact time evolution is shown in Fig. 4.

3.1 Model I

Figure 3 shows the growth and dynamical evolution of protoplanetary embryos in one simulation of Model I (see Table 1 for the definition of our models) with t_start = 0.5 Myr and S_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 typicallylarger in the cold regions of the disk in our model. At 1 Myr, the largest planetary embryo outside the snow line is about 1 M_⊕. At 2 Myr, the mass of the largest planetary embryo is about ~7 M_⊕. 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. 2015a). 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 Fig. 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 Fig. 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 ~10 M_⊕.

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 t = 0.5 Myr in Fig. 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 0.7 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 Fig. 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 1 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. 2015b). At the end of the gas disk phase, this simulation produced six planets with masses larger than 2 M_⊕ 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., ~ 1−1.5 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 Fig. 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 Fig. 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 ~ 1.5 Myr protoplanetary embryos have already reached masses of greater than ~6 M_⊕ 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 Fig. 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 ~ 50 M_⊕. Thus, the final masses of planetary objects in the simulations of Figs. 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 Fig. 5 – with masses larger than ~ 20 M_⊕ – 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 thisissue more carefully. In this work, whenever a planetary embryo becomes more massive than ~ 15 M_⊕ 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 ≤15 M_⊕ during the gas disk phase, in Sect. 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 ~ 15 M_⊕ 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_peb = 0.2 to 2.5. Each panel of Fig. 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_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_peb = 0.2 and S_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 thelong-term dynamical evolution of a similar population of rocky embryos). Earth-mass or more-massive planets at 1–2 AU produced for S_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_peb = 0.2 and t_start = 0.5 Myr, for example, features an integrated pebble flux of ~0.2 × 194 M_⊕ = 39 M_⊕. A simulation with S_peb = 1 and t_start = 0.5 Myr results in an integrated pebble flux of ~1 × 194 M_⊕ (see Fig. 1). In both cases, only a fraction of the integrated pebble flux is used to build planets. In Fig. 6 these numbers are 17, 25, 32, and 30% for S_peb = 0.2, 0.4, 1, and 2.5, respectively.

Figure 6 also shows that our nominal pebble flux S_peb = 1 results in the formation of planets that do indeed migrate to the inner edge of the disk (assumed at ~ 0.3 AU in these simulations). Planets inside 0.7 AU in simulations with S_peb = 1 have masses lower than ~10−15 M_⊕ (see also Fig. 4). A further factor of 2.5 increase in the pebble flux results in the formation of multiple planets with masses as large as ~40−50 M_⊕ (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_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_peb = 2.5 are more massive than the Solar System ice giants by a factor of a few. Although not shown in Fig. 6, for completeness, we also inspected the results of simulations considering higher scaling pebble fluxes of S_peb = 5 and 10. In these simulations, the final planets reaching the inner edge of the disk are as massive as 100 M_⊕ 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.

Fig. 4 Dynamical evolution of protoplanetary embryos in a simulation of Model I where tstart = 0.5 Myr, Speb = 1 and the size of the rocky pebbles is Rockypeb = 0.1 cm (same simulation as in Fig. 3). Top left-hand panel: evolution of the semi-major axes. Top right-hand panel: evolution of the eccentricities. Bottom left-hand and right-hand panels: mass growth and the orbital inclinations, respectively. The dashed vertical line shows the end of the gas disk dispersal (corresponding to Ṁ ~ 10−9M⊙ yr−1 in Eq. (2)). Colored lines show the final planets inside 0.7 AU. The gray line shows final planets and leftover protoplanetary embryos with orbits outside 0.7 AU. Leftover embryos are those with masses smaller than ~ 0.1 M⊕. Finally, the black line shows collided or ejected objects over the course of the simulation.

Fig. 5 Dynamical evolution of a simulation similar to the one shown in Fig. 4, but using Speb = 2.5. The larger pebble flux leads to faster growth and thus larger planets via pebble accretion. Even before the disk dissipates, planetary embryos collide and form massive bodies.

Fig. 6 Final masses of protoplanetary embryos in simulations of Model I where tstart = 0.5 Myr and the size of the rocky pebbles is Rockypeb = 0.1 cm for different pebble fluxes Speb. Each panel shows the outcome of ten different simulations with slightly different initial conditions. Each final planetary object is represented by a colored dot where the color represents its final ice mass fraction. Planetary objects belonging to the same simulation are connected by lines. The integrated pebble fluxes are ~ 39 M⊕, ~ 78 M⊕, ~ 194 M⊕, and ~ 485 M⊕ for Speb = 0.2,  0.4,  1,  and 2.5, respectively. The efficiency of pebble accretion (fraction of integrated pebble flux used to build planets) is about 17, 25, 32, and 30% for Speb = 0.2, 0.4, 1, and 2.5, respectively.

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 t_start = 3.0 Myr, i.e., towards the end of the lifetime of the disk (5 Myr). For the nominal pebble flux scaling S_peb = 1, the integrated drifted pebble flux for a disk with t_start = 0.5 Myr is about ~194 M_⊕. In a disk with t_start = 3.0 Myr and S_peb = 1 the total integrated pebble flux is about ~30 M_⊕ (see Fig. 1). This is the total reservoir of mass available for planetary growth.

We note that the results from Models I and II shown in Figs. 6 and 7, respectively, are not dramatically different: with S_peb = 0.2 and t_start = 0.5 Myr versus S_peb = 1 and t_start = 3.0 Myr. For example, the largest planets in the top left-hand panel of Fig. 7 are around 1 AU and have masses of between 1 and 2 M_⊕. Planets around 1 AU in the top left-hand panel of Fig. 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 (~ 39 M_⊕ for Model-I with S_peb = 0.2 and ~ 30 M_⊕ for Model II with S_peb = 1.0; note that they have different 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 Fig. 6 (S_peb = 0.4) with the top right panel of Fig. 7 (S_peb = 2.5). We note that these two cases have also comparable integrated pebble fluxes (~ 75 M_⊕). In the top right-hand panel of Fig. 6, the typical final mass of planetary embryos around 10 AU is sub-Mars mass while those in top right-hand panel of Fig. 7 have masses of about ~ 2 M_⊕, at least a factor of 20 larger.

Figure 7 shows that increasing pebble flux from S_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 Fig. 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_peb = 10 (see bottom left-hand panel of Fig. 7) promotes the formation of planetary embryos with masses greater than ~ 20 M_⊕. This reinforces our previous finding of Fig. 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).

Fig. 7 Same as Fig. 6 but using Model II where tstart = 3.0 Myr and the size of the rocky pebbles is Rockypeb = 0.1 cm for different pebble fluxes Speb, as annotatedin each panel. As these simulations start at tstart = 3.0 Myr, the integrated pebble fluxes are ~30 M⊕, ~ 75 M⊕, ~ 150 M⊕, and ~ 300 M⊕ for Speb = 1,  2.5,  5,  and 10, respectively. The efficiency of pebble accretion is about 49, 43, 39, and 33% for Speb = 1,  2.5,  5,  and 10, respectively.

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 1). 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 Fig. 8 the size of the silicate pebble is Rocky_peb = 0.1 cm while in Fig. 8 this is Rocky_peb = 1 cm. In both cases t_start = 0.5 Myr.

In simulations of Fig. 8, small planetary embryos are initially distributed between ~ 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 theinner 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 ~ 2 M_⊕ protoplanetary embryos becoming stranded between 1 and 2 AU (see e.g., top left-hand panel of Fig. 8 and the shape of themigration zone in Fig. 3).

Figure 8 shows that simulations with S_peb = 1 failed to form multiple planets with masses larger than a few M_⊕ of any composition anchored at the inner edge of the disk. Increased pebble fluxes (S_peb = 2.5, 5, and 10) successfully promote the formation and delivery of icy planetary embryos with masses of ~ 5−10 M_⊕ to the inner edge of the disk at 0.3 AU. However, only simulations with S_peb = 10 (botton right-hand panel of Fig. 8) produce final rocky protoplanetary embryos with masses larger than ~ 1 M_⊕, 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 (Rocky_peb = 1 cm) for different pebble fluxes. As expected, comparing the results of Fig. 8 and Fig. 9 it is clear that the growth of rocky planetary embryos is far more efficient when Rocky_peb = 1 cm for all pebble fluxes. Interestingly, only simulations with S_peb ≳ 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_peb = 10 produce at least a few planetary embryos as large as ~20 M_⊕. 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 (Rocky_peb = 1 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 0.01 M_⊕ planetary seed at 1 AU growing in a 0.5 Myr-old disk where S_peb = 5 and Rocky_peb = 1 cm is very similar to that of an identical seed in a simulation with S_peb = 5, Rocky_peb = 0.1 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 α_set ≃ α∕10 (see Sect. 2.3). We discuss this issue again in Sect. 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_peb and t_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 Figs. 8 and 9), seeds starts only inside the snow line where pebbles are typically smaller than beyond the snow line (even when Rocky_peb = 1.0 cm) and, consequently, the efficiency of pebble accretion is reduced. Second, in simulations of Model-III with Rocky_peb = 0.1 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_peb = 5 and 10 result in a large mass in pebbles available for planetary growth (if t_start = 0.5 Myr they imply ~970 M_⊕ and ~1940 M_⊕, 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 reportedpebble 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.

Fig. 8 Same as Fig. 7 except that we now use Model III, where the embryos are initially spread from 0.3 to 2.0 AU. The integrated pebble fluxes are ~ 194 M⊕, ~ 485 M⊕, ~ 970 M⊕, and ~ 1940 M⊕ for Speb = 1,  2.5,  5,  and 10, respectively. The efficiency of pebble accretion is about 4.3, 2.1, 1.6, and 1.4% for Speb = 1,  2.5,  5,  and 10, respectively.Clearly, the most massive bodies in each simulation have a large water fraction.

Fig. 9 Same as Fig. 8, except that we now use rocky pebbles of Rockypeb = 1 cm. This leads to enhanced growth of the inner embryos in contrast to Fig. 8. The efficiency of pebble accretion is about 5.6, 4.8, 4.1, and 3.4% for Speb = 1,  2.5,  5,  and 10, respectively.

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 Sect. 4 below). The rocky planets at ~1 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 Sect. 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 × 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 Sects. 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 Sect. 7. Paper I can therefore be seen as a subcase where some process, for example theformation 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 aimto model the formation of super-Earth systems like those observed. As shown in the previous section, different combinations of parameters (t_start, S_peb, and Rocky_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 ~ 1 M_⊕ and ~ 15 M_⊕. 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 2, which also presents the integrated pebble flux of each setup, which depends on S_peb and t_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 r_in = 0.1 AU rather than at ~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 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 1. After the gas disk dispersal, simulations are continued for an additional ~ 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 t_start = 3.0 Myr, S_peb = 5, and Rocky_peb = 0.1 cm. The dynamical evolution during the gas disk phase of the growing planetary embryos is qualitatively similar to those in Figs. 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 ~ 7 M_⊕. After gas dispersal, the simulation is continued for an additional 45 Myr in a gas-free scenario but the resonant chain remains dynamicallystable. The right-hand panels of Fig. 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_start = 3.0 Myr, S_peb = 1, and Rocky_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 Fig. 10. After the gas disk phase (see vertical grayline 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 Figs. 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 ~ 10 M_⊕. There is no clear radial mass ranking in planetary systems produced in our simulations (in contrast, see Fig. 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 Fig. 12 that unstable systems are relatively more spread and have fewer but larger planets. The results of Fig. 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.

Fig. 10 Growth and long-term dynamical evolution of protoplanetary embryos in a simulation of Model II where tstart = 3.0 Myr, Speb = 5, and the size of the rocky pebbles is Rockypeb = 0.1 cm. The top left panel: evolution of semi-major axes. Top right panel: evolution of the eccentricities. Bottom left and right panels: mass growth and the orbital inclinations, respectively. Colored lines show the final planets orbiting inside 0.7 AU. The gray lines show final planets and leftover protoplanetary embryos with orbits outside 0.7 AU. Leftover embryos are those with masses smaller than ~ 0.1 M⊕. Finally, the black lines show collided or ejected objects over the course of the simulation. The dashed vertical line shows the instant of the gas disk dispersal. After gas disk dispersal the resonant chains of super-Earths anchored at the inner edge of the disk remain dynamically stable during the total integration time of 50 Myr.

Fig. 11 Same as Fig. 10, except that we show here a case where the system undergoes a dynamical instability.

Fig. 12 Planetary systems produced in different simulations in a diagram showing semi-major axis versus mass after 50 Myr of integration. Left panel: stable planetary systems. Middle panel: final architecture of selected dynamically unstable systems. Right panel: systems where all planet masses (or m sin (i) in the case of radial velocity detections) are estimated to be lower than 22 M⊕. Only planets within 0.7 AU are shown even though Earth-mass planets may exist in these systems farther out. In these simulations, Rockypeb = 0.1 cm.

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 2 produce a much higherrate 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 ~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 ~ 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 Sect. 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 ~ 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 Fig. 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 linesin Fig. 13). The purple dashed line of Fig. 14 shows that Model I with t_start = 3.0 Myr and S_peb = 5 produces the most compact planetary systems at the end of the gas disk phase. The middle panel of Fig. 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 t_start = 0.5 Myr, S_peb = 10, and Rocky_peb = 0.1 cm, in which fewer than half of the planetary systems become dynamically unstable after gas dispersal (yellow dashed line in Fig. 13). The period ratio distribution of planet pairs for this set of simulations shows that they are theleast 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 dynamicallyunstable 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.

Fig. 13 Cumulative distributions of the last collision epoch in all our unstable systems also with Rockypeb = 0.1 cm. The distributions only account for collisions happening after the gas disk dispersal.

Fig. 14 Dynamical architecture of stable and unstable planetary systems anchored at the disk inner edge at the end of the gas disk dispersal (before the onset of dynamical instabilities). From left to right, the panels show the period ratio distribution, the distribution of the number of planets, and planet mass distributions. These distributions were calculated considering planets inside 1 AU with masses larger than 0.1M⊕.

6 Matching observations

We now evaluate how well our simulations match observations. In Sect. 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 similarmasses 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 (Sect. 6.2) and then the Kepler dichotomy (Sect. 6.3). Finally, we discuss the rocky versus icy nature of super-Earths (Sect. 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 et al. 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. 2011b; 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 2016). 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 consequenceof 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.

Fig. 15 Dynamical architecture of stable and unstable planetary systems anchored at the disk inner edge at the end of the simulations (50 Myr). From left to right, the top panels: cumulative distributions of period ratio, mass, and orbital eccentricity. From left to right, the bottom panels: distributions of orbital inclination and number of planets. Only planets with semi-major axes smaller than 0.7 AU and masses larger than 1 M⊕ are considered. In all considered models Rockypeb = 0.1 cm. The Kepler planet mass distribution inferred from the mass–radius relationship of Wolfgang et al. (2016) is shown as a thick gray line. The uncertainty in the mass distribution is shown by the thin gray lines (M ± σM).

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 R_⊕. 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 R_⊕ 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 Sect. 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 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 (P ≲ 200 days) and with masses larger than 1 M_⊕ 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 4 R_⊕, 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 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 Fig. 15 shows the period ratio distribution of adjacent Kepler planets (R < 4 R_⊕). 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 Fig. 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 ~ 10 M_⊕ and a median value of ~3−4 M_⊕ for all scenarios of Table 2. However, the mass distributions of unstable systems are dramatically different.

The blue dashed line representing unstable systems of Model II with t_start = 3.0 Myr, S_peb = 5, and Rocky_peb = 0.1 cm produced final planets with a median mass of ~7.5 M_⊕. The red dashed line representing unstable systems of Model I with t_start = 0.5 Myr, S_peb = 1, and Rocky_peb = 0.1 cm shows a median planet mass of about ~15 M_⊕, a factor of two larger. This result becomes easier to understand by revisiting Fig. 14. Although the mass distributions of Fig. 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 Fig. 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 $M=2.6{\left(R/R_{\oplus }\right)}^{1.3}$ with a standard deviation ${\sigma }_{\text{M}}=\sqrt{4.41+1.5\left({\text{R/R}}_{\oplus }-1\right)}$. Because we impose a cutoff of R < 4 R_⊕ in the Kepler sample, the maximum mass of Kepler planets inferred from Wolfgang’s mass–radius relation is ~ 18.9 M_⊕. With this relation, the median planet mass in our Kepler sample is 6.4 M_⊕. The gray lines in the top-middle panel of Fig. 15 show the expected mass (M; thick line) and mass dispersion (M ± σ_M; thin lines) distributions of Kepler planets.

The planet-mass distribution of unstable systems of model I (t_start = 0.5 Myr, S_peb = 1, and Rocky_peb = 0.1 cm; red dashed line in Fig. 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 18.9 M_⊕. Planet masses in unstable systems of model I with t_start = 3.0 Myr and S_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 18.9 M_⊕, in good agreement with the expected masses of Kepler planets. Figure 15 also shows a clear trend. Planet pairs in Model I with t_start = 0.5 Myr typically have larger period ratios than planet pairs composed of lower mass planets (Model I with t_start = 3.0 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 dynamicallyexcited. Statistical analyses have inferred the orbital eccentricity and inclination distributions of Kepler planets (Lissauer et al. 2011b; Kane et al. 2012; Tremaine & Dong 2012; Figueira et al. 2012; Fabrycky et al. 2014; Plavchan et al. 2014; Ballard & Johnson 2016; 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 σ_e = 0.1 (e.g., Moorhead et al. 2011)and σ_i = 1.5° (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 Fig. 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 Fig. 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 Fig. 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 ≿ 10 days (which corresponds to orbital distances greater than ~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 ~ 0.1 AU (P ≈ 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_⋆ ∕a, where R_⋆ 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 detected⁴ we store theorbital details of each detect planet and from the detected planet(s) we create a new observed planetary system.

The top left-hand panel of Fig. 17 shows the period ratio distribution of stable and unstable detected systems of Model I with t_start = 0.5 Myr. In the whole of our analysis, we only consider adjacent planet pairs with 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 Fig. 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. 17) because the orbital separation between adjacent planets and mutual inclinations are significantly smaller than in the case of Model I with t_start = 0.5 Myr.

Fig. 16 Orbital period distribution of stable and unstable planetary systems anchored at the disk inner edge at the end of the simulations (50 Myr). The color and style of each line representing our simulations correspond to those shown in the legend of Fig. 15. The thick gray line shows the period distribution of Kepler planets before correction for transit probabilities and detection biases. The thin gray line shows the debiased period distribution following Fressin et al. (2013). Top x-axis: corresponding orbital distance calculated assuming a solar-mass star.

Fig. 17 Left-hand panels: period ratio distribution of simulations of Model-I with tstart = 0.5 Myr (red), Model-I with tstart = 3.0 Myr (purple), and Model-III with tstart = 0.5 Myr (yellow). In all cases Rockypeb=0.1 cm. The black lines show the period ratio distribution of detected planet pairs in synthetic observations of our real systems. Solid lines show stable systems. Dashed lines show unstable systems. The gray line shows the Kepler sample. Right-hand panels: p-values of Kolmogorov-Smirnof tests between the Kepler sample and simulated detection of samples mixing different fractions of stable and unstable systems.

Matching Kepler observations with a mix of stable and unstable systems

Figure 17 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 Fig. 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 Fig. 17 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 Fig. 17 shows that our simulated observations match the Kepler sample if less than ~ 1−5% of systems remain stable. This also holds true for simulated observations of Model I with t_start = 3.0 Myr and Model III (right-second-from-top and bottom panels of Fig. 17). 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 (≲25%) than the one obtained here (≲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 P_out∕P_in < 3 and with masses lower than 18.9 M_⊕. We have relaxed these cutoffs in our analysis because our simulations produced lower mass planets than those of Izidoro et al. (2017; compare Fig. 15 with Fig. 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 (Fig. 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 t_start = 0.5 Myr and S_peb = 1 and model II with t_start = 3.0 Myr and S_peb = 5 are the least favored models overall. The former model produces planets that are systematically too massive (top-middle panel of Fig. 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 t_start = 3.0 Myr and S_peb = 5 and model III with t_start = 0.5 Myr and S_peb = 10.

Each of our favored models has its own caveats. The masses of planets in model I with t_start = 3.0 Myr and S_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. 2011b; 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 2016; 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 N_det = 1. Izidoro et al. (2017) matched the Kepler dichotomy assuming that Kepler planets comprise about ≲ 10% of stable and ≳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 18 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 t_start = 0.5 Myr have two or more planets inside 0.7 AU (red dashed line in the left-hand panel of Fig. 18). The two stable systems of Model I with t_start = 0.5 Myr have six and eight planets in their chains.

Simulated detections of unstable systems of Model I with t_start = 0.5 Myr provide a very good match to the Kepler multiplicity. Simulated observations of unstable systems rise the fraction of single planets from ~ 10%⁵ to about ~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 ≲ 5% of stable systems with ≳95% of unstable ones. Figure 18 shows that mixing 1% stable and 99% unstable systems provides an almost perfect match to the Kepler dichotomy in Model I with t_start = 0.5 Myr. The middle and right-hand panels of Fig. 18 show that Model I with t_start = 3.0 Myr and Model III with t_start = 0.5 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 Fig. 18 show a deficit at 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 Fig. 18). We also reiterate thefact that Model II fails to match Kepler observations period ratio distribution (Fig. 17).

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 19 shows again (as in Fig. 18) the observed planet multiplicity distribution in synthetic transit observationsof a mix of 1% stable and 99% unstable systems (dashed green histograms). The colors filling each N_det dashed green histogram bin show the fractional contributions as a function of true system multiplicity. In the left panel, where t_start = 0.5 Myr, 80% of all the observed systems show a single transit. However, only ~ 5% of these systems have only one planet inside 0.7 AU, ~39% of systems have two planets, ~28% have three planets, and ~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 t_start = 3.0 Myr, the peak at N_det = 1 constitutes ~ 65% of the systems, where ~ 11% of the systems have two planets, ~14% have three planets, ~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 t_start = 3.0 Myr and S_peb = 5 and Model III with t_start = 0.5 Myr and S_peb = 10 do not perfectly match the dichotomy, there may be several false positives in the 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 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. 2015b, 2017). Therefore, Model I with t_start = 3.0 Myr and Model III with t_start = 0.5 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 R_⊕ 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 2) grew mainly from pebbles originating from beyond the snow line. Thus, their final ice-mass fractions are as high as 50% (e.g., Figs. 6 and 7). Also, as suggested by Fig. 8, only Model III simulations with large refractory pebbles (Rocky_peb = 1 cm) or extremely large pebble fluxes produced multiple inner planets with low water-ice contents.

Figure 20 shows the final planetary systems produced in 50 simulations of Model III with t_start = 0.5 Myr, S_peb = 10, and Rocky_peb = 0.1 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, Fig. 20 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 Fig. 20 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 ~ 9 M_⊕ 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 Fig. 20) the innermost planets in several systems are typically rocky. In some cases, rocky super-Earths end up on orbits that are immediatelyadjacent 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 Rocky_peb = 0.1 cm the largest rocky close-in super-Earths in stable systems have typical masses of 2 M_⊕ or lower. Observations of photoevaporated planets show rocky super-Earths up to at least 5 M_⊕. unstable systems of Model III with Rocky_peb = 0.1 cm do produce some ~5 M_⊕ rocky super-Earths but we can produce them far more easily in Model III with Rocky_peb = 1 cm.

We also note that in our unstable systems of Fig. 20, 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 ~1−2 M_⊕ 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 Fig. 20, is at least qualitatively consistent with the findings of Zeng et al. (2019), which suggest that super-Earths larger than 2 R_⊕ 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 Fig. 1) are all rocky. This is roughly consistent with the results of Fig. 20. Our results in Fig. 20 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.

Fig. 18 Number of planets detected by transit in synthetic observations of planetary systems produced in our N-body simulations of Model I with tstart = 0.5 Myr and Speb = 1 (red), ModelI with tstart = 3.0 Myr and Speb = 5 (purple), Model II with tstart = 3.0 Myr and Speb = 5 (blue), and Model III with tstart = 0.5 Myr and Speb = 10 (yellow). Forall cases, Rockypeb = 0.1 cm. Solid lines show stable systems, and dashed lines show unstable systems. The black lines show the results of our synthetic observations. The densely dashed green lines show multiplicity distribution combining different real fractions of stable and unstable systems. We note that the y-axes combine linear and log scaling.

Fig. 19 Number of planets detected by transit in synthetic observations of planetary systems of two sets of simulations of Model I. The colors filling the dashed green histogram bins show the fractional contributions as a function of true system multiplicity. Left and right panels: simulations with tstart = 0.5 Myr, and tstart = 3.0 Myr, respectively. In both cases we use Rockypeb = 0.1 cm. The dashed green histograms show detection distributions considering different real mixing ratios of stable and unstable systems, as in Fig. 18. The over-plotted red numbers show the size of each colored sub-bar representing fractional contributions. We note that the y-axes combine linear and log scaling.

Fig. 20 Planetary systems produced in Model III with tstart = 0.5 Myr, Speb = 10, and Rockypeb = 0.1 cm from 50 simulations at the end of the simulations (50 Myr). Two panels are shown. Left-hand panel: stable systems (50%). Right-hand panel: unstable ones (50%). Planets are represented by dots. The sizes of the dots scale with mass as m1∕3. The color ofeach dot indicates its ice-mass fraction.

Fig. 21 Final masses of protoplanetary embryos in simulations where the gas disk is not evolving and the snow line is kept fixed. The size of the rocky pebbles is Rockypeb = 0.1 cm for different pebble fluxes Speb. Each panel shows the outcome of five different simulations with slightly different initial conditions. Each final planetary object is represented by a colored dot where the color represents its final ice mass fraction. Planetary objects belonging to the same simulation are connected by lines.

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 to ~ 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 Sect. 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 21 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 (t_start = 0.5 Myr). We keep theinitial disk structure throughout the entire simulation, and therefore the snow line is kept fixed at about ~ 3 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 ~ 1 AU.

Figure 21 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 22 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 Fig. 22.

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., ~ 10 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 Rocky_peb = 0.1 cm and S_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.

Fig. 22 Final masses of protoplanetary embryos in simulations where type-I migration is neglected but the disk evolves as in our nominal simulations. The size of rocky pebbles is Rockypeb = 0.1 cm for different pebble fluxes Speb. Each panel shows the outcome of five different simulations with slightly different initial conditions. Each final planetary object is represented by a colored dot where the color represents its final ice mass fraction. Planetary objects belonging to the same simulation are connected by lines.

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

Our goal forthis 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 Fig. 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 (Rocky_peb = 1 cm) and S_peb = 5. We do not re-use the results of Fig. 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_in = 0.1 AU. We note that we stop the simulations of Fig. 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 Fig. 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 23 shows our results.

The upper-panels of Fig. 23 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 Fig. 23 shows that 88% of these systems became dynamically unstable after gas dispersal. This is different from simulations of Model III with Rocky_peb = 0.1 cm where only ~ 50% of the systems became unstable⁶. This results in a much better though still imperfect match to observations (see bottom-right panel of 23).

The bottom-left panel of Fig. 23 shows the period ratio distribution of planet pairs of stable and unstable systems andalso 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 t_start = 0.5 Myr, Rocky_peb = 0.1 cm, and S_peb = 10 (see bottompanels of Fig. 17 and left-panel of 18). 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 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 Sect. 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 ~ 0.1 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, Fig. 24 shows all unstable systems produced in simulations of Model III with Rocky_peb = 1 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 Fig. 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 layerwith 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.

Fig. 23 Statistics of simulations of Model III with tstart = 0.5 Myr, Speb = 5, and Rockypeb = 1 cm. Solid lines show stable systems, and dashed lines show unstable systems. The black lines show the results of our synthetic observations. The thick gray lines represent the Kepler planets. Upper-left panel: period ratio distribution of simulations and observations. Upper-right panel: KS-tests between the Kepler period ratio sample and our simulated detections of distributions mixing different real fractions of stable and unstable systems. Bottom-left panel: planet multiplicity distribution.The densely dashed green line shows the detected multiplicity distribution for a mix of 1% stable and 99% unstable systems. Bottom-right panel: distribution of the last collision epoch after gas dispersal. About 88% of these systems became unstable after gas dispersal.

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 ~ 8% of the planetary systems have the innermost planet (>1 R_⊕) 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. 2015b,a).

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 Rocky_peb = 0.1 cm and S_peb = 1 (see top-left of Fig. 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.

Fig. 24 Planetary systems produced in Model III with tstart = 0.5 Myr, Speb = 5, and Rockypeb = 1 cm from 46 unstable systems at the end of the simulations (300 Myr). The sizes of the dots scale with mass as m1∕3. The color ofeach dot indicates its ice-mass fraction.

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 ~ 39 M_⊕ – have final typical masses of ~ 2 M_⊕ or lower. Simulations where the total pebble reservoir is of ~194 M_⊕ 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 ~ 194 M_⊕ to about ~ 485 M_⊕ is enough to bifurcate the evolution of our planetary systems from one leading to typical super-Earth-mass planets (≲ 15 M_⊕) 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 isnot 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 ~ 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 ≲2% stable with ≳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 R_⊕ 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 757 448-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, 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, areassociated with shadowed regions, where the disk’s aspect ratio H/r decreases with orbital distance (Bitsch et al. 2013b, 2014; Bitsch et al. 2015a; 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 Γ_tot is the sum of the Lindblad torque Γ_L and the corotation torque Γ_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 $${\Gamma }_{\text{tot}}={\Gamma }_{\text{L}}{\Delta }_{\text{L}}+{\Gamma }_{\text{C}}{\Delta }_{\text{C}},$$(A.1)

where Δ_L and Δ_C are rescaling factors accounting for torque reduction due to the orbital eccentricity and inclination of the planet. The Lindblad torque is $${\Gamma }_{\text{L}}=(-2.5-1.7\beta +0.1x)\frac{{\Gamma }_{\text{0}}}{{\gamma }_{\text{eff}}},$$(A.2)

where β 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): $$x=-\frac{\partial ln {\Sigma }_{\text{gas}}}{\partial ln r}, \beta =-\frac{\partial ln T}{\partial ln r}.$$(A.3)

The Lindblad torque reduction Δ_L (Cresswell & Nelson 2008) is $$\begin{array}{lll}{\Delta }_{\text{L}}\hfill & =\hfill & [P_{\text{e}}+\frac{P_{\text{e}}}{\left|P_{\text{e}}\right|}\times \{0.07\left(\frac{i}{h}\right)+0.085{\left(\frac{i}{h}\right)}^4\hfill \\ \hfill & \hfill & {-0.08\left(\frac{e}{h}\right){\left(\frac{i}{h}\right)}^2\}]}^{-1},\hfill \end{array}$$(A.4)

where $$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},$$(A.5)

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 $$\begin{array}{lll}{\Gamma }_{\text{C}}\hfill & =\hfill & {\Gamma }_{\text{c,hs,baro}}F(p_{\nu })G(p_{\nu })+(1-K(p_{\nu })){\Gamma }_{\text{c,lin,baro}}\hfill \\ \hfill & \hfill & +{\Gamma }_{\text{c,hs,ent}}F(p_{\nu })F(p_{\chi })\sqrt{G(p_{\nu })G(p_{\chi })}\hfill \\ \hfill & \hfill & +\sqrt{(1-K(p_{\nu }))(1-K(p_{\chi })}{\Gamma }_{\text{c,lin,ent}} ,\hfill \end{array}$$(A.6)

where $${\Gamma }_{\text{c,hs,baro}}=1.1\left(\frac{3}{2}-x\right)\frac{{\Gamma }_{\text{0}}}{{\gamma }_{\text{eff}}},$$(A.7) $${\Gamma }_{\text{c,lin,baro}}=0.7\left(\frac{3}{2}-x\right)\frac{{\Gamma }_{\text{0}}}{{\gamma }_{\text{eff}}},$$(A.8) $${\Gamma }_{\text{c,hs,ent}}=7.9\xi \frac{{\Gamma }_{\text{0}}}{{\gamma }_{\text{eff}}^2},$$(A.9)

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

In the previous equations, ξ = β − (γ − 1)x is the negative of the entropy gradient and γ = 1.4 is the adiabatic index. The scaling torque ${\Gamma }_0={(q/h)}^2{\Sigma }_{\text{gas}}{r}^4{\Omega }_{\text{k}}^2$ is calculated at the planet’s location. Here, q is the planet–star mass ratio, and Ω_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 $$p_{\nu }=\frac{2}{3}\sqrt{\frac{r^2{\Omega }_{\text{k}}}{2\pi \nu }{x}_{\text{s}}^3},$$(A.11)

where x_s is the nondimensional half width of the horseshoe region, $$x_{\text{s}}=\frac{1.1}{{\gamma }_{\text{eff}}{}^{1/4}}\sqrt{\frac{q}{h}},$$(A.12)

and $$p_{\chi }=\frac{2}{3}\sqrt{\frac{r^2{\Omega }_{\text{k}}}{2\pi \chi }{x}_{\text{s}}^3},$$(A.13)

where χ is the thermal diffusion coefficient which reads as $$\chi =\frac{16\gamma (\gamma -1)\sigma T^4}{3\kappa {\rho }^2{(hr)}^2{\Omega }_{\text{k}}^2}.$$(A.14)

In Eq. (A.14), ρ is the gas volume density, κ is the gas disk opacity, and σ 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).

γ_eff depends on $$Q=\frac{2\chi }{3h^3{r}^2{\Omega }_{\text{k}}},$$(A.15)

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

F, G, and K take the form $$F(p)=\frac{1}{1+{\left(\frac{p}{1.3}\right)}^2,}$$(A.17) $$G(p)=\{\begin{array}{ll}\frac{16}{25}{\left(\frac{45\pi }{8}\right)}^{\frac{3}{4}}{p}^{\frac{3}{2}},\hfill & \text{if}\$\text{p}<\sqrt{\frac{\text{8}}{\text{45}\pi }}\$\hfill \\ 1-\frac{9}{25}{\left(\frac{8}{45\pi }\right)}^{\frac{4}{3}}{p}^{-\frac{8}{3}},\hfill & \text{otherwise}.\hfill \end{array},$$(A.18)

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

Here, p is either p_ν or p_χ as defined above.

The co-orbital torque reduction is $${\Delta }_{\text{C}}=\exp \left(\frac{e}{e_{\text{f}}}\right)\left\{1-\tan \text{h}\left(\frac{i}{h}\right)\right\},$$(A.20)

where e_f is given by Fendyke & Nelson (2014) as $$e_{\text{f}}=0.5h+\mathrm{0.01.}$$(A.21)

We note that for sufficiently large eccentricities, Δ_L may assumenegative values favoring unrealistic reversal of migration. In extreme situations, where Δ_L < 0, we impose Δ_L = Δ_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: $$t_{\text{m}}=-\frac{L}{{\Gamma }_{\text{tot}}},$$(A.22)

where L and Γ_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). t_e and t_i are the orbital eccentricity and inclination damping timescales, respectively: $$\begin{array}{lll}{t}_e\hfill & =\hfill & \frac{t_{\text{ave}}}{0.780}(1-0.14{\left(\frac{e}{h}\right)}^2+0.06{\left(\frac{e}{h}\right)}^3\hfill \\ \hfill & \hfill & +0.18\left(\frac{e}{h}\right){\left(\frac{i}{h}\right)}^2),\hfill \end{array}$$(A.23) $$\begin{array}{lll}{t}_i\hfill & =\hfill & \frac{t_{\text{wave}}}{0.544}(1-0.3{\left(\frac{i}{h}\right)}^2+0.24{\left(\frac{i}{h}\right)}^3\hfill \\ \hfill & \hfill & +0.14{\left(\frac{e}{h}\right)}^2\left(\frac{i}{h}\right)).\hfill \end{array}$$(A.24)

In Eqs. (A.23) and (A.24), $$t_{\text{wave}}=\left(\frac{M_{\odot }}{m_{\text{p}}}\right)\left(\frac{M_{\odot }}{{\Sigma }_{\text{gas}}{a}^2}\right)h^4{\Omega }_{\text{k}}^{-1},$$(A.25)

with M_⊙, a_p, m_p, i, and 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 (a_m), orbital eccentricity (a_e), and inclination damping (a_i): $$a_m=-\frac{v}{t_m},$$(A.26) $$a_e=-2\frac{(v.r)r}{r^2{t}_e},$$(A.27)

and $$a_i=-2\frac{v_z}{t_i}k,$$(A.28)

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

References

All Tables

All Figures

Fig. 1 Pebble flux in our nominal disk model. The integrated mass in drifted pebbles (Ipeb) is shown on the left-hand vertical axis as a function of disk age (tdisk). Fluxes are calculated at 5 AU. The right-hand vertical axis shows the instantaneous pebble flux decreasing as the disk dissipates. Simulations of Paper I and III also feature the decay of the pebble flux. All our simulations start with a tdisk of at least 0.5 Myr, thus the pebble production line is already past the initial positions of our outermost seeds (~ 60 AU; see Table 1). This plot is generated considering Speb = 1 in Eq. (13) but a simple rescaling of these curves accounts for other considered pebble fluxes. Both curves correspond to the pebble flux beyond the snow line. The pebble flux inside the snow line is reduced by a factor of two because of the mass sublimation of the ice component of the pebbles when they cross the snow line. The gas flux and the integrated gas flux in our gas disk model are shown by the green and blue lines, respectively. We note that in the plot we rescaled the gas flux by a factor of 1/30, for clear comparison with the pebble flux.

In the text

Fig. 2 Top: gas and pebble surface densities as a function of orbital distance at different disk ages. Bottom: pebble sizes and Stokes number as a function of orbital distance at different disk ages. In both panels, Rockypeb = 1 mm and Speb = 1, which correspond to our nominal values. The sharp drop in the Stokes number in the inner regions of the disk is due to the snow line transition where we set the size of silicate pebbles to 1 mm.

In the text

Fig. 3 Evolution of an example simulation from model I. The panels represent snapshots of the growth of protoplanetary embryos in a simulation with tstart = 0.5 Myr, Speb = 1, and a size of the rocky pebbles equal to Rocky_peb = 0.1 cm. The vertical and horizontal axes represent mass and semi-major axis, respectively. The location of the disk snow line (blue vertical dashed line) and the pebble isolation mass (black dashed line) are also shown for reference as the disk evolves. Planetary embryos growing inside the snow line accrete silicate pebbles while those outside the snow line accrete icy pebbles. The gray solid lines delimit the region of outward migration. The disk inner edge is shown at ~ 0.3 AU, where planets are trapped as well. The color of each dot gives the ice mass fraction. The size of each dot scales as m1∕3, where m is the planetary mass. The exact time evolution is shown in Fig. 4.

In the text

Fig. 4 Dynamical evolution of protoplanetary embryos in a simulation of Model I where tstart = 0.5 Myr, Speb = 1 and the size of the rocky pebbles is Rockypeb = 0.1 cm (same simulation as in Fig. 3). Top left-hand panel: evolution of the semi-major axes. Top right-hand panel: evolution of the eccentricities. Bottom left-hand and right-hand panels: mass growth and the orbital inclinations, respectively. The dashed vertical line shows the end of the gas disk dispersal (corresponding to Ṁ ~ 10−9M⊙ yr−1 in Eq. (2)). Colored lines show the final planets inside 0.7 AU. The gray line shows final planets and leftover protoplanetary embryos with orbits outside 0.7 AU. Leftover embryos are those with masses smaller than ~ 0.1 M⊕. Finally, the black line shows collided or ejected objects over the course of the simulation.

In the text

	Fig. 5 Dynamical evolution of a simulation similar to the one shown in Fig. 4, but using Speb = 2.5. The larger pebble flux leads to faster growth and thus larger planets via pebble accretion. Even before the disk dissipates, planetary embryos collide and form massive bodies.
In the text

Fig. 6 Final masses of protoplanetary embryos in simulations of Model I where tstart = 0.5 Myr and the size of the rocky pebbles is Rockypeb = 0.1 cm for different pebble fluxes Speb. Each panel shows the outcome of ten different simulations with slightly different initial conditions. Each final planetary object is represented by a colored dot where the color represents its final ice mass fraction. Planetary objects belonging to the same simulation are connected by lines. The integrated pebble fluxes are ~ 39 M⊕, ~ 78 M⊕, ~ 194 M⊕, and ~ 485 M⊕ for Speb = 0.2,  0.4,  1,  and 2.5, respectively. The efficiency of pebble accretion (fraction of integrated pebble flux used to build planets) is about 17, 25, 32, and 30% for Speb = 0.2, 0.4, 1, and 2.5, respectively.

In the text

Fig. 7 Same as Fig. 6 but using Model II where tstart = 3.0 Myr and the size of the rocky pebbles is Rockypeb = 0.1 cm for different pebble fluxes Speb, as annotatedin each panel. As these simulations start at tstart = 3.0 Myr, the integrated pebble fluxes are ~30 M⊕, ~ 75 M⊕, ~ 150 M⊕, and ~ 300 M⊕ for Speb = 1,  2.5,  5,  and 10, respectively. The efficiency of pebble accretion is about 49, 43, 39, and 33% for Speb = 1,  2.5,  5,  and 10, respectively.

In the text

Fig. 8 Same as Fig. 7 except that we now use Model III, where the embryos are initially spread from 0.3 to 2.0 AU. The integrated pebble fluxes are ~ 194 M⊕, ~ 485 M⊕, ~ 970 M⊕, and ~ 1940 M⊕ for Speb = 1,  2.5,  5,  and 10, respectively. The efficiency of pebble accretion is about 4.3, 2.1, 1.6, and 1.4% for Speb = 1,  2.5,  5,  and 10, respectively.Clearly, the most massive bodies in each simulation have a large water fraction.

In the text

	Fig. 9 Same as Fig. 8, except that we now use rocky pebbles of Rockypeb = 1 cm. This leads to enhanced growth of the inner embryos in contrast to Fig. 8. The efficiency of pebble accretion is about 5.6, 4.8, 4.1, and 3.4% for Speb = 1,  2.5,  5,  and 10, respectively.
In the text

Fig. 10 Growth and long-term dynamical evolution of protoplanetary embryos in a simulation of Model II where tstart = 3.0 Myr, Speb = 5, and the size of the rocky pebbles is Rockypeb = 0.1 cm. The top left panel: evolution of semi-major axes. Top right panel: evolution of the eccentricities. Bottom left and right panels: mass growth and the orbital inclinations, respectively. Colored lines show the final planets orbiting inside 0.7 AU. The gray lines show final planets and leftover protoplanetary embryos with orbits outside 0.7 AU. Leftover embryos are those with masses smaller than ~ 0.1 M⊕. Finally, the black lines show collided or ejected objects over the course of the simulation. The dashed vertical line shows the instant of the gas disk dispersal. After gas disk dispersal the resonant chains of super-Earths anchored at the inner edge of the disk remain dynamically stable during the total integration time of 50 Myr.

In the text

	Fig. 11 Same as Fig. 10, except that we show here a case where the system undergoes a dynamical instability.
In the text

Fig. 12 Planetary systems produced in different simulations in a diagram showing semi-major axis versus mass after 50 Myr of integration. Left panel: stable planetary systems. Middle panel: final architecture of selected dynamically unstable systems. Right panel: systems where all planet masses (or m sin (i) in the case of radial velocity detections) are estimated to be lower than 22 M⊕. Only planets within 0.7 AU are shown even though Earth-mass planets may exist in these systems farther out. In these simulations, Rockypeb = 0.1 cm.

In the text

	Fig. 13 Cumulative distributions of the last collision epoch in all our unstable systems also with Rockypeb = 0.1 cm. The distributions only account for collisions happening after the gas disk dispersal.
In the text

Fig. 14 Dynamical architecture of stable and unstable planetary systems anchored at the disk inner edge at the end of the gas disk dispersal (before the onset of dynamical instabilities). From left to right, the panels show the period ratio distribution, the distribution of the number of planets, and planet mass distributions. These distributions were calculated considering planets inside 1 AU with masses larger than 0.1M⊕.

In the text

Fig. 15 Dynamical architecture of stable and unstable planetary systems anchored at the disk inner edge at the end of the simulations (50 Myr). From left to right, the top panels: cumulative distributions of period ratio, mass, and orbital eccentricity. From left to right, the bottom panels: distributions of orbital inclination and number of planets. Only planets with semi-major axes smaller than 0.7 AU and masses larger than 1 M⊕ are considered. In all considered models Rockypeb = 0.1 cm. The Kepler planet mass distribution inferred from the mass–radius relationship of Wolfgang et al. (2016) is shown as a thick gray line. The uncertainty in the mass distribution is shown by the thin gray lines (M ± σM).

In the text

Fig. 16 Orbital period distribution of stable and unstable planetary systems anchored at the disk inner edge at the end of the simulations (50 Myr). The color and style of each line representing our simulations correspond to those shown in the legend of Fig. 15. The thick gray line shows the period distribution of Kepler planets before correction for transit probabilities and detection biases. The thin gray line shows the debiased period distribution following Fressin et al. (2013). Top x-axis: corresponding orbital distance calculated assuming a solar-mass star.

In the text

Fig. 17 Left-hand panels: period ratio distribution of simulations of Model-I with tstart = 0.5 Myr (red), Model-I with tstart = 3.0 Myr (purple), and Model-III with tstart = 0.5 Myr (yellow). In all cases Rockypeb=0.1 cm. The black lines show the period ratio distribution of detected planet pairs in synthetic observations of our real systems. Solid lines show stable systems. Dashed lines show unstable systems. The gray line shows the Kepler sample. Right-hand panels: p-values of Kolmogorov-Smirnof tests between the Kepler sample and simulated detection of samples mixing different fractions of stable and unstable systems.

In the text

Fig. 18 Number of planets detected by transit in synthetic observations of planetary systems produced in our N-body simulations of Model I with tstart = 0.5 Myr and Speb = 1 (red), ModelI with tstart = 3.0 Myr and Speb = 5 (purple), Model II with tstart = 3.0 Myr and Speb = 5 (blue), and Model III with tstart = 0.5 Myr and Speb = 10 (yellow). Forall cases, Rockypeb = 0.1 cm. Solid lines show stable systems, and dashed lines show unstable systems. The black lines show the results of our synthetic observations. The densely dashed green lines show multiplicity distribution combining different real fractions of stable and unstable systems. We note that the y-axes combine linear and log scaling.

In the text

Fig. 19 Number of planets detected by transit in synthetic observations of planetary systems of two sets of simulations of Model I. The colors filling the dashed green histogram bins show the fractional contributions as a function of true system multiplicity. Left and right panels: simulations with tstart = 0.5 Myr, and tstart = 3.0 Myr, respectively. In both cases we use Rockypeb = 0.1 cm. The dashed green histograms show detection distributions considering different real mixing ratios of stable and unstable systems, as in Fig. 18. The over-plotted red numbers show the size of each colored sub-bar representing fractional contributions. We note that the y-axes combine linear and log scaling.

In the text

	Fig. 20 Planetary systems produced in Model III with tstart = 0.5 Myr, Speb = 10, and Rockypeb = 0.1 cm from 50 simulations at the end of the simulations (50 Myr). Two panels are shown. Left-hand panel: stable systems (50%). Right-hand panel: unstable ones (50%). Planets are represented by dots. The sizes of the dots scale with mass as m1∕3. The color ofeach dot indicates its ice-mass fraction.
In the text

Fig. 21 Final masses of protoplanetary embryos in simulations where the gas disk is not evolving and the snow line is kept fixed. The size of the rocky pebbles is Rockypeb = 0.1 cm for different pebble fluxes Speb. Each panel shows the outcome of five different simulations with slightly different initial conditions. Each final planetary object is represented by a colored dot where the color represents its final ice mass fraction. Planetary objects belonging to the same simulation are connected by lines.

In the text

Fig. 22 Final masses of protoplanetary embryos in simulations where type-I migration is neglected but the disk evolves as in our nominal simulations. The size of rocky pebbles is Rockypeb = 0.1 cm for different pebble fluxes Speb. Each panel shows the outcome of five different simulations with slightly different initial conditions. Each final planetary object is represented by a colored dot where the color represents its final ice mass fraction. Planetary objects belonging to the same simulation are connected by lines.

In the text

Fig. 23 Statistics of simulations of Model III with tstart = 0.5 Myr, Speb = 5, and Rockypeb = 1 cm. Solid lines show stable systems, and dashed lines show unstable systems. The black lines show the results of our synthetic observations. The thick gray lines represent the Kepler planets. Upper-left panel: period ratio distribution of simulations and observations. Upper-right panel: KS-tests between the Kepler period ratio sample and our simulated detections of distributions mixing different real fractions of stable and unstable systems. Bottom-left panel: planet multiplicity distribution.The densely dashed green line shows the detected multiplicity distribution for a mix of 1% stable and 99% unstable systems. Bottom-right panel: distribution of the last collision epoch after gas dispersal. About 88% of these systems became unstable after gas dispersal.

In the text

	Fig. 24 Planetary systems produced in Model III with tstart = 0.5 Myr, Speb = 5, and Rockypeb = 1 cm from 46 unstable systems at the end of the simulations (300 Myr). The sizes of the dots scale with mass as m1∕3. The color ofeach dot indicates its ice-mass fraction.
In the text