Growth after the streaming instability: from planetesimal accretion to pebble accretion

TL;DR

This study investigates how planetesimals grow into super-Earths via streaming instability, focusing on the effects of initial size distributions, pebble flux, disk turbulence, and pebble properties.

Contribution

It demonstrates that diverse initial planetesimal size distributions enable super-Earth formation within disk lifetimes, highlighting the importance of dynamical friction and pebble flux.

Findings

Mono-dispersed 400 km planetesimals cannot form Earth-mass planets.

Poly-dispersed and bimodal distributions facilitate super-Earth formation.

Higher pebble flux, lower turbulence, and larger Stokes number favor super-Earth growth.

Abstract

Streaming instability is a key mechanism in planet formation, clustering pebbles into planetesimals. It is triggered at a particular disk location where the local volume density of solids exceeds that of the gas. After their formation, planetesimals can grow by feeding from other planetesimals in the birth ring as well as by accreting inwardly drifting pebbles from the outer disk. To investigate the growth of planetesimals at a single location by the streaming instability, we test the conditions under which super-Earths are able to form within the lifetime of the gaseous disk. We modify the \texttt{Mercury} N-body code to trace the growth and dynamical evolution of a swarm of planetesimals at the ice line for a solar-mass star. Three distributions of planetesimal sizes are investigated: (i) a mono-dispersed population of 400 km radius planetesimals, (ii) a poly-dispersed populations of planetesimals from 200 km up to 1000 km, (iii) a bimodal distribution with a single runaway body and a swarm of smaller, 100 km size planetesimals. The mono-disperse population of 400 km size planetesimals cannot form $\gtrsim$ Earth mass protoplanets. Their velocity dispersions are quickly excited, which suppresses both planetesimal and pebble accretion. Planets can form from the poly-dispersed and bimodal distributions. In these circumstances, the two-component nature damps the random velocity of the large embryo by small planetesimals' dynamical friction, allowing the embryo to accrete pebbles efficiently when it approaches $10^{-2}$ Earth mass. We find that super-Earth planets are preferred to form when the pebble mass flux is higher, the disk turbulence is lower, or the Stokes number of the pebbles is higher.

Tables

Table 1. Table 1: Simulations set-up in Sect.  3.1 , Sect.  3.2 and Sect.  3.3 . The mass and semi-major axis of the most massive planet are given at t = 1 𝑡 1 t=1 Myr.

Name	initial planetesimal size	inclination	$m_p\left(M_{\oplus }\right)$	$a_p$ (AU)
run_md_test	mono-size	No	$4.0$	$0.1$
run_md_1	mono-size	Yes	$0.004$	$2.7$
run_md_2	mono-size	Yes	$0.009$	$2.7$
run_md_3	mono-size	Yes	$0.02$	$2.6$
run_pd_1	poly-size	Yes	$1.3$	$1.7$
run_pd_2	poly-size	Yes	$1.6$	$1.5$
run_pd_3	poly-size	Yes	$1.4$	$1.8$
run_sp_1	one embryo + small planetesimals	Yes	$1.3$	$1.8$
run_sp_2	one embryo + small planetesimals	Yes	$1.4$	$1.6$
run_sp_3	one embryo + small planetesimals	Yes	$1.3$	$1.7$

Table 2. Table 2: Simulations for the parameter study in Sect.  4 . The last two columns represent the final mass and the time when the planet reaches the inner 0.1 0.1 0.1 AU.

Name (description)	${\dot{M}}_{\mathrm{peb}}$ ($M_{\oplus }{\mathrm{Myr}}^{-1}$)	${\alpha }_t$	${\tau }_s$	$\dot{a}/a$	$m_{p,f}\left(M_{\oplus }\right)$	$t_f\left(\mathrm{Myr}\right)$
run_fid (fiducial)	$100$	${10}^{-3}$	$0.1$	yes	$4.0$	$1.2$
run_nmig (no migration)	$100$	${10}^{-3}$	$0.1$	no	$3.8$	$1.2$
run_alpha (low turbulence)	$100$	${10}^{-4}$	$0.1$	yes	$4.1$	$0.9$
run_tau (low Stokes number)	$100$	${10}^{-3}$	$0.03$	yes	$4.8$	$1.5$
run_hpeb (high pebble flux)	$200$	${10}^{-3}$	$0.1$	yes	$6.8$	$0.7$
run_lpeb (low pebble flux)	$50$	${10}^{-3}$	$0.1$	yes	$2.3$	$2.1$

Equations

$${\Sigma }_{\mathrm{gas}}={\Sigma }_{gas0}{\left(\frac{r}{1\mathrm{AU}}\right)}^{-1},$$

$$h_{\mathrm{gas}}=h_{gas0}{\left(\frac{r}{1\mathrm{AU}}\right)}^{1/4},$$

$$T=\frac{\mu h_{\mathrm{gas}}^{2}G{M}_{\star }}{r{R}_{\mathrm{g}}}=T_0{\left(\frac{r}{1\mathrm{AU}}\right)}^{-1/2},$$

$$\eta =-\frac{h_{\mathrm{gas}}^2}{2}\frac{\partial \log P_{\mathrm{gas}}}{\partial \log r}={\eta }_0{\left(\frac{r}{1\mathrm{AU}}\right)}^{1/2},$$

$${\mathbf{a}}_{\mathrm{drag}}=-\left(\frac{3C_{\mathrm{D}}{\rho }_{\mathrm{gas}}}{8R_{\mathrm{p}}{\rho }_{\bullet }}\right)v_{\mathrm{rel}}{\mathbf{v}}_{\mathrm{rel}},$$

$$\begin{array}{rl}{t}_{\mathrm{m}}\simeq & 0.5\left(\frac{M_{\star }}{m_{\mathrm{p}}}\right)\left(\frac{M_{\star }}{{\Sigma }_{\mathrm{gas}}{a}_{\mathrm{p}}^2}\right)\left(\frac{h_{\mathrm{gas}}^2}{{\Omega }_{\mathrm{K}}}\right)\simeq 5\times 10^5{\left(\frac{m_{\mathrm{p}}}{1{\mathrm{M}}_{\oplus }}\right)}^{-1}\\ & {\left(\frac{{\Sigma }_{gas0}}{500\mathrm{gc}{\mathrm{m}}^{-1}}\right)}^{-1}{\left(\frac{h_{gas0}}{3.3\times 10^{-2}}\right)}^2\left(\frac{a_{\mathrm{p}}}{2.7\mathrm{AU}}\right)\mathrm{yr},\end{array}$$

$${\mathbf{a}}_{\mathrm{m}}=-\frac{\mathbf{v}}{t_{\mathrm{m}}},{\mathbf{a}}_{\mathrm{e}}=-2\frac{(\mathbf{v}\cdot \mathbf{r})\mathbf{r}}{r^2{t}_{\mathrm{e}}},{\mathbf{a}}_{\mathrm{i}}=-\frac{{\mathbf{v}}_{\mathrm{z}}}{t_{\mathrm{i}}},$$

$$H_{\mathrm{peb}}=\sqrt{\frac{{\alpha }_{\mathrm{t}}}{{\alpha }_{\mathrm{t}}+{\tau }_{\mathrm{s}}}}{H}_{\mathrm{gas}}.$$

$${\epsilon }_{2D}=\frac{0.32}{\eta }\sqrt{\frac{m_{\mathrm{p}}\Delta v}{M_{\star }{v}_{\mathrm{K}}{\tau }_s}}{f}_{\mathrm{set}}$$

$${\epsilon }_{3D}=0.39\frac{m_{\mathrm{p}}}{\eta h_{\mathrm{peb}}{M}_{\star }}{f}_{\mathrm{set}}^2$$

$$f_{\mathrm{set}}=\exp \left[-0.5(\Delta v/v_{\ast }{)}^2\right]$$

$${\dot{M}}_{\mathrm{i},\mathrm{peb}}=\{\begin{array}{rlr}& {\dot{M}}_{\mathrm{peb}}& i=1\\ & {\dot{M}}_{\mathrm{peb}}\prod _{k=1}^{i-1}(1-{\epsilon }_{\mathrm{k},\mathrm{PA}}),& i\ge 2\end{array}$$

$$\begin{array}{rl}& {\tau }_{\mathrm{vs}}\simeq 3\times 10^5\left(\frac{i}{0.005}\right){\left[{\left(\frac{e}{0.01}\right)}^2+{\left(\frac{i}{0.005}\right)}^2\right]}^{\frac{3}{2}}{\left(\frac{R_{\mathrm{p}}}{400\mathrm{km}}\right)}^{-3}\\ & {\left(\frac{{\Sigma }_{\mathrm{p}}}{10\mathrm{g}\mathrm{c}{\mathrm{m}}^{-2}}\right)}^{-1}{\left(\frac{{\rho }_{\bullet }}{1.5\mathrm{g}\mathrm{c}{\mathrm{m}}^{-3}}\right)}^{-1}{\left(\frac{a}{2.7\mathrm{AU}}\right)}^{-1/2}yr.\end{array}$$

$$\begin{array}{rl}{\tau }_{\mathrm{drag}}\simeq & 10^7\left(\frac{R_{\mathrm{p}}}{400\mathrm{km}}\right){\left(\frac{{\Sigma }_{gas0}}{500\mathrm{g}\mathrm{c}{\mathrm{m}}^{-2}}\right)}^{-1}\left(\frac{h_{gas0}}{0.033}\right)\\ & {\left[{\left(\frac{e}{10^{-2}}\right)}^2+{\left(\frac{i}{5\times 10^{-3}}\right)}^2+{\left(\frac{\eta }{2.5\times 10^{-3}}\right)}^2\right]}^{-1/2}\\ & \left(\frac{{\rho }_{\bullet }}{1.5\mathrm{g}\mathrm{c}{\mathrm{m}}^{-3}}\right){\left(\frac{a}{2.7\mathrm{AU}}\right)}^{11/4}yr.\end{array}$$

$$\begin{array}{rl}& {\tau }_{\mathrm{tidal}}\simeq \left(\frac{M_{\star }}{m_{\mathrm{p}}}\right)\left(\frac{M_{\star }}{{\Sigma }_{\mathrm{gas}}{a}_{\mathrm{p}}^2}\right)\left(\frac{h_{\mathrm{gas}}^4}{{\Omega }_{\mathrm{K}}}\right)\simeq 7\times 10^7{\left(\frac{R_{\mathrm{p}}}{400\mathrm{km}}\right)}^{-3}\\ & {\left(\frac{{\Sigma }_{gas0}}{500\mathrm{g}\mathrm{c}{\mathrm{m}}^{-2}}\right)}^{-1}{\left(\frac{h_{gas0}}{0.033}\right)}^4{\left(\frac{{\rho }_{\bullet }}{1.5\mathrm{g}\mathrm{c}{\mathrm{m}}^{-3}}\right)}^{-1}{\left(\frac{a}{2.7\mathrm{AU}}\right)}^{3/2}yr.\end{array}$$

$$\Delta v\lesssim v_{\ast }={\left(\frac{M_p/M_{\star }}{{\tau }_{\mathrm{s}}}\right)}^{1/3}{v}_{\mathrm{K}}.$$

$$\Delta W=C_{\mathrm{w}}\sqrt[3]{\left(N{\Omega }_{\mathrm{K}}t/e^2\right)}{R}_{\mathrm{H}}^2/a,$$

$$\frac{N_{>}(m)}{N_{\mathrm{tot}}}={\left(\frac{m}{m_{\min }}\right)}^{-p}\exp \left[{\left(\frac{m_{\min }}{m_{\mathrm{e}}}\right)}^{\beta }-{\left(\frac{m}{m_{\mathrm{e}}}\right)}^{\beta }\right].$$

$$\begin{array}{rl}{\tau }_{\mathrm{rg}}& =C_{\mathrm{rg}}\frac{{\rho }_{\bullet }{R}_0}{{\Sigma }_{\mathrm{p}}{\Omega }_{\mathrm{K}}}=10^5\left(\frac{C_{\mathrm{rg}}}{0.1}\right)\left(\frac{R_0}{100\mathrm{km}}\right)\\ & {\left(\frac{{\Sigma }_{\mathrm{p}}}{10\mathrm{g}\mathrm{c}{\mathrm{m}}^{-2}}\right)}^{-1}{\left(\frac{a}{2.7\mathrm{AU}}\right)}^{3/2}yr.\end{array}$$

$$\begin{array}{rl}{R}_{rg/oli}\simeq & 850{\left(\frac{C_{\mathrm{rg}}}{0.1}\right)}^{3/7}{\left(\frac{R_0}{100\mathrm{km}}\right)}^{3/7}{\left(\frac{a}{2.7\mathrm{AU}}\right)}^{5/7}{\left(\frac{{\Sigma }_{\mathrm{p}}}{1\mathrm{g}\mathrm{c}{\mathrm{m}}^{-2}}\right)}^{2/7}km.\end{array}$$

$$\Delta W=\sqrt{\frac{C_{\mathrm{rg}}{C}_{\mathrm{w}}^3{R}_{\mathrm{H}}^3}{2R_0}}.$$

$$\begin{array}{r}{R}_{rg/oli}\simeq 1000{\left(\frac{C_{\mathrm{rg}}}{0.1}\right)}^{3/7}{\left(\frac{R_0}{100\mathrm{km}}\right)}^{1/7}{\left(\frac{M_{\mathrm{tot},\mathrm{pl}}}{0.039{\mathrm{M}}_{\oplus }}\right)}^{2/7}km.\end{array}$$

$$\begin{array}{rl}{\tau }_{\mathrm{PA},2D}=& \frac{m_{\mathrm{p}}}{{\epsilon }_{\mathrm{PA},2D}{\dot{M}}_{\mathrm{peb}}}\simeq 9\times 10^4{f}_{\mathrm{set}}^{-1}{\left(\frac{m_{\mathrm{p}}}{0.05{\mathrm{M}}_{\oplus }}\right)}^{1/3}\\ & {\left(\frac{{\tau }_{\mathrm{s}}}{0.1}\right)}^{1/3}\left(\frac{\eta }{2.5\times 10^{-3}}\right){\left(\frac{{\dot{M}}_{\mathrm{peb}}}{100{\mathrm{M}}_{\oplus }/Myr}\right)}^{-1}yr.\end{array}$$

$$\begin{array}{rl}{f}_{\mathrm{set}}& =\exp \left[-0.07{\left(\frac{\eta }{2.5\times 10^{-3}}\right)}^2{\left(\frac{m_{\mathrm{p}}}{0.01{\mathrm{M}}_{\oplus }}\right)}^{-2/3}{\left(\frac{{\tau }_{\mathrm{s}}}{0.1}\right)}^{2/3}\right].\end{array}$$

$$\begin{array}{rl}{\tau }_{\mathrm{PA},3D}=& \frac{m_{\mathrm{p}}}{{\epsilon }_{\mathrm{PA},3D}{\dot{M}}_{\mathrm{peb}}}\simeq 9\times 10^4{f}_{\mathrm{set}}^{-2}\left(\frac{h_{\mathrm{peb}}}{4.2\times 10^{-3}}\right)\\ & \left(\frac{\eta }{2.5\times 10^{-3}}\right){\left(\frac{{\dot{M}}_{\mathrm{peb}}}{100{\mathrm{M}}_{\oplus }/Myr}\right)}^{-1}yr.\end{array}$$

Full text

*Aims. * To investigate the growth of planetesimals into protoplanets at a single location by the streaming instability. For a solar-mass star, we test the conditions under which super-Earths are able to form within the lifetime of the gaseous disk.

*Methods. * We modify the Mercury N-body code to trace the growth and dynamical evolution of a swarm of planetesimals at a distance of $2.7$ AU from the star. The code simulates gravitational interactions and collisions among planetesimals, gas drag, type I torque, and pebble accretion. Three distributions of planetesimal sizes are investigated: (i) a mono-dispersed population of $400$ km radius planetesimals, (ii) a poly-dispersed populations of planetesimals from $200$ km up to $1000$ km, (iii) a bimodal distribution with a single runaway body and a swarm of smaller, $100$ km size planetesimals.

*Results. * The mono-disperse population of $400$ km size planetesimals cannot form $\gtrsim$ Earth mass protoplanets. Their eccentricities and inclinations are quickly excited, which suppresses both planetesimal accretion and pebble accretion. Planets can form from the poly-dispersed and bimodal distributions. In these circumstances, it is the two-component nature that damps the random velocity of the large embryo by small planetesimals’ dynamical friction, allowing the embryo to accrete pebbles efficiently when it approaches $10^{-2}\ \rm M_{\oplus}$. Accounting for migration, close-in super-Earth planets form. Super-Earth planets are preferred to form when the pebble mass flux is higher, the disk turbulence is lower, or the Stokes number of the pebbles are higher.

*Conclusions. * For the single site planetesimal formation scenario, a two-component mass distribution with a large embryo and small planetesimals promotes planet growth, first by planetesimal accretion and then by pebble accretion of the most massive protoplanet. Planetesimal formation at single locations such as ice lines naturally leads to super-Earth planets by the combined mechanisms of planetesimal accretion and pebble accretion.

1 Introduction

In protoplanetary disks, micron-sized dust grains coagulate into pebbles of mm-cm sizes (Dominik & Tielens, 1997; Birnstiel et al., 2012; Krijt et al., 2016; Pérez et al., 2015; Tazzari et al., 2016). But further growth is suppressed by bouncing or fragmentation due to the increasing compactification in collisions (Güttler et al., 2010; Zsom et al., 2010). In addition, these pebbles also drift too fast compared to their growth (Weidenschilling, 1977) so that the particles cannot cross the meter size barrier even if they would stick perfectly (Birnstiel et al., 2012; Lambrechts & Johansen, 2014), unless the pebbles can remain fluffy during the growth (Okuzumi et al., 2012; Kataoka et al., 2013). The subsequent growth of these pebbles to planetesimals is still not well understood in planet formation theory (see Johansen et al. 2014 for a review).

The streaming instability mechanism provides a promising solution by concentrating drifting pebbles due to a locally enhanced solid-to-gas ratio. Once the threshold of solid-to-gas ratio is satisfied, the pebble clumps can directly collapse into planetesimals (Youdin & Goodman, 2005; Johansen et al., 2007, 2009; Bai & Stone, 2010). The characteristic size of these planetesimals approximate a few hundred kilometres (Johansen et al., 2012, 2015; Simon et al., 2016; Schäfer et al., 2017; Simon et al., 2017; Abod et al., 2018).

The subsequent growth after planetesimal formation by streaming instability has not been well studied. These newly born planetesimals would interact with each other. In the classical planetesimal accretion scenario (see Raymond et al. (2014); Izidoro & Raymond (2018) for reviews), the gravitational interactions among these planetesimals lead to orbital crossings, scatterings and collisions. Initially, the velocity dispersions of planetesimals are not strongly excited and remain modest. In this stage the accretion is super linear (${\rm d}m_{\rm p}/{\rm d}t\propto m_{\rm p}^{\gamma}$ with $\gamma>1$), which is termed ‘runaway growth’ (Greenberg et al., 1978; Wetherill & Stewart, 1989; Ida & Makino, 1993; Kokubo & Ida, 1996; Rafikov, 2004; Ida & Lin, 2004). It means that the massive body has a faster accretion rate and therefore will get more massive quickly. However, this stage cannot last forever. Since the growing massive bodies would stir the random velocities of neighbouring small planetesimals, the accretion rates of the massive bodies slow down and turn into a self-limiting mode. This phase is called ‘oligarchic growth’ where ${\rm d}m_{\rm p}/{\rm d}t\propto m_{\rm p}^{\gamma}$ with $\gamma<1$ (Kokubo & Ida, 1998). This phase is characterized by a decreasing mass ratio of two adjacent massive runaway bodies (Lissauer, 1987; Kokubo & Ida, 2000; Thommes et al., 2003; Ormel et al., 2010).

In addition to planetesimal accretion, planetesimals formed by streaming instability can accrete inward drifting pebbles. A planetesimal may capture a fraction of the pebbles which cross its orbit (Ormel & Klahr, 2010; Lambrechts & Johansen, 2012). This is known as pebble accretion (see recent reviews by Johansen & Lambrechts, 2017; Ormel, 2017). Even if only a fraction of pebbles are able to be accreted by planets, the pebble accretion rate can still be high for two reasons. First, the accretion cross section is significantly enhanced by gas drag (Ormel & Klahr, 2010); and second, a large flux of pebbles grow and drift inward from the outer regions of disks (Birnstiel et al., 2012; Lambrechts & Johansen, 2014). Pebble accretion can be classified into 2D/3D regimes (Ormel & Klahr, 2010; Morbidelli et al., 2015). When the pebble accretion radius is larger than the pebble scale height, the accretion is in the $2$D regime, where ${\rm d}m_{\rm p}/{\rm d}t\propto m_{\rm p}^{2/3}$ (Lambrechts & Johansen, 2014, Hill regime). On the other hand, when the pebble scale height exceeds the pebble accretion radius, only pebbles with heights smaller than the accretion radius can be accreted. Therefore, the accretion rate in this $3$D regime is reduced compared to $2$D, and ${\rm d}m_{\rm p}/{\rm d}t\propto m_{\rm p}$ (Ida et al., 2016).

In general, the efficacy of the pebble accretion mechanism to grow planet(esimals) depends on many variables related to the properties of the disk, pebble, and planet (the eccentricity, the inclination, and the mass of the planet, the pebble size, the disk turbulence, etc.) A key quantity is the pebble accretion efficiency $\varepsilon_{\rm PA}$ (Guillot et al., 2014; Lambrechts & Johansen, 2014), defined as the number of pebbles that are accreted divided by the total number of pebbles that the disk supplies. Recently we have computed $\varepsilon_{\rm PA}$ under general circumstances (Liu & Ormel, 2018; Ormel & Liu, 2018). For instance, when the eccentricity and inclination of the planet become high, $\varepsilon_{\rm PA}$ drops significantly compared to planets on coplanar and circular orbits, because pebbles are approaching at too high velocity to the planet.

Progress in planet formation requires an improved understanding under which conditions the mass growth is dominated by accreting pebbles or planetesimals. In this work, our goal is to investigate the growth of planetesimals after their formation by streaming instability at a single disk location (e.g., the H2O ice line). Hansen (2009) already proposed that the architecture of the Solar system’s terrestrial planets can be explained when planetesimals grow in a narrow annulus. In his model the width of the annulus is $0.3$ AU, much wider than our planetesimal forming zone (see Sect. 2). Furthermore, Hansen (2009) focused on the planetesimal accretion in a gas-free environment. Our work instead considers the growth just after the streaming instability in gas-rich disk phase.

Constrained by the size distribution of planetesimals in the asteroid belt, Morbidelli et al. (2009) concluded that their born size should be large ($\gtrsim 100$ km), while Weidenschilling (2011) argued that planetesimals starting from sub-km-sized still cannot be ruled out. Kenyon & Bromley (2010) studied the formation of ice planets beyond $30$ AU and found that protoplanets grow more efficiently with smaller planetesimal sizes. Motivated by streaming instability simulations, our adopted initial sizes are typical $\gtrsim 100$km. In the context of combined planetesimal and pebble accretion, Johansen et al. (2015) studied the growth of asteroids using a statistical approach, and concluded that massive protoplanets or even super-Earths can form by a combination of pebble accretion, planetesimal accretion and giant impacts.

In order to study planet formation from a narrow ring of planetesimals, we employ direct N-body techniques. Growth can be classified into two phases: (A) the planetesimal accretion dominated phase and (B) the pebble accretion dominated phase. The N-body approach is necessary to treat phase A and the transition from phase A to phase B. The Mercury N-body code has been modified to include gas drag, type I torque and pebble accretion. Three different types of initial planetesimal size distributions are investigated. We find that a two-component mass distribution (large embryo + small planetesimals) is needed to grow a massive planet. This condition could arise either from the high mass tail distribution of planetesimals formed by streaming instability or be the result of runaway growth of a population of small planetesimals.

The paper is structured as follows. In Sect. 2, we outline our model and the implementation of the N-body code. In Sect. 3 three initial size distributions of planetesimals are investigated, including a mono-dispersed population in Sect. 3.1, a poly-dispersed population in Sect. 3.2, and a single runaway body plus a swarm of small planetesimals in Sect. 3.3. In Sect. 4 we investigate the influence of different parameters in the pebble accretion dominated growth phase (phase B). The key results are summarized in Sect. 5.

2 Method

The hypothesis of this paper is that planetesimals only form at a specific location by streaming instability, which requires a locally enhanced pebble density (Carrera et al., 2015; Yang & Johansen, 2014; Yang et al., 2017, 2018). For instance, Ros & Johansen (2013); Schoonenberg & Ormel (2017); Dra̧żkowska & Alibert (2017) have proposed that the ice line could be such a place since the water vapor inside the ice line will diffuse back and re-condense onto ice pebbles, enriching the solid to gas density ratio. We focus on planetesimals that form at the ice line ($r_{\rm ice}=2.7\ \rm AU$ based on the disk model in Sect. 2.1) in this paper. But the following results and applications can be scaled to other locations where the streaming instability condition is realized, e.g., the inner edge of the zone (Chatterjee & Tan, 2014; Hu et al., 2018) or a distant location due to the FUV photoevaporation (Carrera et al., 2017). In order to generalize our results, we did not include the specific ice line effects such as a reduced pebble size and pebble flux inside of the ice line due to sublimation.

This scenario is illustrated in Figure 1. We consider the situation where streaming instability operates to quickly spawn planetesimals at the initial time of our simulations. A population of planetesimals has formed in a narrow annulus of relative width equal to the normalized pressure gradient ($\Delta r=\eta r$, see Fig. 1 and further discussion in Sect. 2.3). These planetesimals can grow further in two ways: by coagulation among themselves (planetesimal accretion) or by sweeping up pebbles that drift in from the outer disk (pebble accretion).

2.1 Disk model

The surface density and the aspect ratio of the gas disk are assumed to be

[TABLE]

[TABLE]

where $\Sigma_{\rm gas0}$ and $h_{\rm gas0}$ are the gas surface density and the aspect ratio at $1$ AU and $r$ is the distance to the central star. The aspect ratio index of 1/4 assumes an optically thin stellar irradiated disk. For simplicity, we neglect that in the early stages disks might be hotter due to viscous accretion (Garaud & Lin, 2007; Bitsch et al., 2015). Here $\Sigma_{\rm gas0}=500\rm\ g\,cm^{-2}$ and $h_{\rm gas0}=0.033$ (Hayashi, 1981) are adopted as the default values in this paper. Therefore, the disk temperature $T$ and $\eta$ are

[TABLE]

and

[TABLE]

where $G$ is the gravitational constant, $R_{\rm g}$ $=$ $8.31\times 10^{7}\rm erg\,mol^{-1}\,K^{-1}$ is the gas constant, $M_{\star}$ is the stellar mass, $\mu=2.34$ is the mean molecular weight and $P_{\rm gas}$ is the gas pressure in the disk. In the above equations $T_{0}=280\rm\ K$ and $\eta_{0}=1.5\times 10^{-3}$ are the values at $1$ AU.

Although the adopted surface density power-law index is shallower than minimum mass solar nebula model, this profile ($\Sigma_{\rm g}\propto r^{-1}$) is more consistent with disk observations (Andrews et al., 2009). The fiducial value of $\Sigma_{\rm gas0}$ is chosen to be the same as Lambrechts & Johansen (2014), and the pebble flux is calibrated accordingly in Sect. 2.2.2. Based on the above disk model, the ice line ($T\sim 170\rm K$) is located at $r_{\rm ice}=2.7$ AU.

2.2 Simulation setup

We use numerical N-body simulations to study the planetesimal growth after the streaming instability during the gas-rich disk phase. We have adopted the Mercury code (Chambers, 1999) and used the Bulirsch-Stoer integrator. An initial timestep of $3$ days and the integration accuracy parameter of $10^{-12}$ are chosen.

For the N-body part, collisions between bodies are treated as inelastic mergers that conserve the linear momentum. The fragmentation/restitution (Leinhardt & Stewart, 2012; Mustill et al., 2018) is not taken into account in this work. Since eccentricities are initially low, the perfect merger assumption is appropriate. For instance, the impact velocity of $400$-km-sized planetesimals is lower than the escape velocity when their eccentricities are lower than a few times $10^{-2}$. Fragmentation among planetesimals will become important after embryos form and stir the planetesimals more vigorously. However, by then pebble accretion is already expected to commence, rendering planetesimal fragmentation irrelevant.

The effects of the disk gas on the planetesimals/planets, such as gas drag, type I torque, eccentricity and inclination damping are taken into account by applying effective forces in Sect. 2.2.1. In Sect. 2.2.2, we implement the pebble accretion prescriptions by Liu & Ormel (2018) and Ormel & Liu (2018), which accounts for the effects of the planet’s eccentricity, inclination and the disk turbulence. The modified code uniquely simulates planet-planet, planet-disk and planet-pebble interactions.

2.2.1 Gas drag and type I migration torque

Small embryos and planetesimals experience the aerodynamic gas drag (Adachi et al., 1976),

[TABLE]

where the drag coefficient $C_{\rm D}=0.5$, $v_{\rm rel}$ is the relative velocity between the planetesimal and the gas, $\bm{v}_{\mathrm{rel}}=\bm{v}-\bm{v}_{\mathrm{gas}}$, where $\bm{v}_{\mathrm{gas}}$ equals $v_{\rm K}(1-\eta)$ in the azimuthal direction, $v_{\rm K}$ is the Keplerian velocity, $\rho_{\rm gas}$ is the local gas density, $\rho_{\bullet}$ and $R_{\rm p}$ are the internal density (assumed to be $1.5\rm\ g\,cm^{-3}$, with half water and half silicate rock) and the physical radius of the planetesimal.

Large embryos and low-mass planets feel the gravitational torques from the disk gas (called type I migration, Goldreich & Tremaine (1979); Kley & Nelson (2012); Baruteau et al. (2014)). The characteristic migration timescale for a planet on a circular orbit is (Cresswell & Nelson, 2008)

[TABLE]

where $m_{\rm p}$ and $a_{\rm p}$ are the mass and the semimajor axis of the planet, respectively, and $\Omega_{\rm K}$ is the Keplerian angular velocity at the planet location.

The accelerations acting on planets due to type I migration torque are

[TABLE]

where $\bm{r}=(x,y,z)$, $\bm{v}=(v_{\rm x},v_{\rm y},v_{\rm z})$ are the position and the velocity vectors of the planet. In the above expression, $t_{\rm m}$, $t_{\rm e}$ and $t_{\rm i}$ are the characteristic type I migration, eccentricity and inclination damping timescales from Eqs. (13), (11) and (12) of Cresswell & Nelson (2008). We note that the torque prescription is based on an isothermal disk and planets always migrate inward. The effect of the unsaturated corotation torque due to the thermal diffusion in a radiative disk (Paardekooper et al., 2010, 2011; Bitsch et al., 2015; Brasser et al., 2017), the dynamical torque (Paardekooper, 2014; McNally et al., 2018), and the heating torque from planet gas accretion (Benítez-Llambay et al., 2015; Masset, 2017; Chrenko et al., 2017) are not taken into account in this work. The mentioned migration prescriptions are implemented in all our model runs except for run_sp_nmig in Sect. 4 (see Table 2). In that case we neglect the semimajor axis damping (set $\bm{a_{\rm m}=0}$) but still consider the eccentricity and inclination damping in Eq. (7) due to the type I torque. This setup would mimic the situation when the planet is situated at a net zero-torque location in the disk where the negative Lindblad torque and positive corotation torque are cancelled out.

2.2.2 Pebble accretion

The pebble-sized particles drift inwards across the protoplanetary disk. The radial drift velocity is $v_{\rm r}=-2\eta v_{\rm K}\tau_{\rm s}/(1+\tau_{\rm s}^{2})$ (Weidenschilling, 1977) where $\tau_{\rm s}$ is the dimensionless stopping time (Stokes number). The drift speed is determined by the aerodynamical size of the pebble $\tau_{\rm s}$ and $\eta$. Based on Lambrechts & Johansen (2014) and Schoonenberg et al. (2018), the pebble mass flux ($\dot{M}_{\rm peb}=2\pi rv_{\rm r}\Sigma_{\rm peb}$) is proportional to the disk pebble surface density, but weakly dependent on time. In the default model we adopt a constant pebble flux of $\dot{M}_{\rm peb}=100\ \rm M_{\oplus}/Myr$ (consistent with Eq. (14) of Lambrechts & Johansen (2014)) and neglect the time dependence for simplicity. A lower $\dot{M}_{\rm peb}$ means an intrinsic metal-deficient/less massive disk and the planetesimals/planets take longer time to grow by pebble accretion. On the other hand, the pebble flux cannot be too high as the density ratio between pebbles and gas would otherwise exceed unity ($\rho_{\rm peb}/\rho_{\rm gas}>1$) and the streaming instability will be triggered over the entirely disk. The above inequality can be written as $\Sigma_{\rm peb}/H_{\rm peb}>\Sigma_{\rm gas}/H_{\rm gas}$ where the scale height of the pebble is given by (Youdin & Lithwick, 2007)

[TABLE]

The pebble flux therefore needs to be smaller than $2\pi rv_{\rm r}\Sigma_{\rm gas}H_{\rm peb}/H_{\rm gas}$. Our default value is well below this limit. We also explore two additional pebble fluxes $\dot{M}_{\rm peb}=200\ \rm M_{\oplus}Myr^{-1}$ and $50\ \rm M_{\oplus}Myr^{-1}$ in Sect. 4.

The pebble mass can be measured from (sub)millimeter dust continuum emission from the young protoplanetary disks. The inferred values are from a few Earth mass to a few hundreds of Earth mass (Ricci et al., 2010; Andrews et al., 2013; Ansdell et al., 2017), which is correlated with the gas disk accretion rates (Manara et al., 2016). For a typical T Taurs star with a disk accretion rate of $10^{-8}\ \rm M_{\odot}\,yr^{-1}$, the dust mass is $\sim 100\ \rm M_{\oplus}$ (Fig.1 of Manara et al. (2016)), consistent with our adopted fiducial value.

In Eq. (8) $\alpha_{\rm t}$ is the coefficient of turbulent gas diffusivity, which is different from the concept of turbulent viscosity $\alpha_{\nu}$. We nevertheless note that the above two quantities are approximated the same when the turbulence is driven by magnetorotational instability (Johansen & Klahr, 2005; Zhu et al., 2015; Yang et al., 2018). The value of $\alpha_{\rm t}$ can be constrained from the molecule line broadening measurements (Flaherty et al., 2015, 2017) and the level of dust settling (Pinte et al., 2016). We adopt a fiducial value of $\alpha_{\rm t}=10^{-3}$ and test the influence of a lower turbulent disk ($\alpha_{\rm t}=10^{-4}$) in Sect. 4.

For pebbels we adopt a fiducial aerodynamical size of $\tau_{\rm s}=0.1$. This chosen value is consistent with the sophisticated dust coagulation model (Birnstiel et al., 2012) and disk observations (Tazzari et al., 2016). A lower $\tau_{\rm s}=0.03$ is also explored in Sect. 4.

A fraction of pebbles will be accreted onto a planetesimal when pebbles drift through its orbit. The pebble accretion efficiency ($\varepsilon_{\rm PA}=\dot{M}_{\rm PA}/\dot{M}_{\rm peb}$) is taken from Liu & Ormel (2018) and Ormel & Liu (2018). This quantity depends on the disk properties ($\tau_{\rm s}$, $\eta$, $\alpha_{\rm t}$) and the planet properties ($m_{\rm p},a,e,i$). In the limit of 2D and 3D pebble accretion these are given by

[TABLE]

and

[TABLE]

respectively, where $\Delta v$ is the relative velocity between the planet and the pebble, $h_{\mathrm{peb}}=H_{\rm peb}/r$ is the aspect ratio of the pebble disk. In Eq. (10) we have already assume that $i<h_{\rm peb}$. The above expressions include a modulation factor,

[TABLE]

where $v_{\ast}=(M_{p}/\tau_{s}M_{\star})^{1/3}v_{\rm K}$ 111A more general form of $f_{\rm set}$ can be found in Eq. (35) of Ormel & Liu (2018), considering the pebble’s velocity distribution in a three-dimensional turbulent disk.. Physically, when the pebble-planet encounter is too fast compared to the coupling time between the gas and the pebble, gas drag is no longer effective to aid the planet to capture pebbles. Therefore, $f_{\rm set}\ll 1$ and pebble accretion fails (Visser & Ormel, 2016).

In the multi-planetesimal system we consider the filtering of flux when pebbles drift through these planetesimals. Therefore, the pebble flux of the body $i$ is given by

[TABLE]

where $i$ is ordered for bodies from the furthest to closest in terms of semimajor axis and $\varepsilon_{i,\mathrm{PA}}$ is the pebble accretion efficiency of the $i^{\mathrm{th}}$ body.

We focus on the formation of protoplanets with $m_{\rm p}\gtrsim 1\ \rm M_{\oplus}$ (progenitors of super-Earths). In our simulations the forming planets are still not massive enough to inverse the local gas pressure gradient and truncate the inward drift of pebbles (Lambrechts et al., 2014; Bitsch et al., 2018; Ataiee et al., 2018). Therefore, we do not implement the termination of pebble accretion.

2.3 Planetesimal initial condition

In the streaming instability mechanism, pebbles accumulate into dense filaments. The typical width of the filament is $\Delta r\simeq\eta r_{\rm ice}$ (Yang & Johansen, 2016; Li et al., 2018). The threshold condition to trigger gravitational collapse requires $\rho_{\rm peb}\simeq\rho_{\rm gas}$. The total solid mass available to build planetesimals is therefore $2\pi r_{\rm ice}\Delta r\Sigma_{\rm peb}(r_{\rm ice})=2\pi r_{\rm ice}\Delta r\Sigma_{\rm gas}(r_{\rm ice})H_{\rm peb}/H_{\rm gas}$. Simon et al. (2016) find that the planetesimal generation efficiency approximates $50\%$ (half of the pebbles convert into planetesimals, see their Fig. 10a). We note that this value might be also dependent on the local metallicity, Stokes number of the pebbles, and the disk turbulence. However, due to a lack of detailed streaming instability simulations exploration, we still use this number to approximate the total mass in planetesimals. With our adopted disk parameters, this amounts to $0.039\ M_{\oplus}$.

3 Scenarios with different initial sizes

Many works have simulated streaming instability numerically, finding that it spawns planetesimals of typical size of several hundreds kilometers, albeit with a considerable uncertainty regarding the precise shape of the size distribution (Johansen et al., 2015; Simon et al., 2016; Schäfer et al., 2017; Abod et al., 2018). However, this initial distribution affects the planet growth, as it will determine the duration of the planetesimal-dominated growth phase (phase A). Here we consider three scenarios for the initial planetesimal distribution in the following subsections: (i) a mono-dispersed population of big planetesimals; (ii) a population of planetesimals of various sizes; and (iii) a two-component population of one large body among many small planetesimals. The first and third scenario reflects standard practice where planetesimals typically start from a fixed size, whereas the second more closely follows the streaming instability results. While the first and second case can be integrated directly, the third scenario employs a superparticle approach to handle the large number of small bodies with N-body techniques.

In this section we investigate the question that provided in total $100\ \rm M_{\oplus}$ of pebbles ( $100\ \rm M_{\odot}yr^{-1}$ of $\dot{M}_{\rm peb}$ over $1$ Myr), what is the mass growth of planetesimals for the above three different initial size distributions.

3.1 Mono-dispersed initial conditions

In this scenario we consider that all planetesimals start with an equal size ($400$ km in radius). The mass growth and dynamical evolution are simulated after the formation at a single site location. We conduct both simulations when the planetesimals orbits are coplanar and inclined. In the former case planetesimals are all initiated at $z=0$ and, by construction, remain in the disk midplane. This (as we will see) is not realistic, but comparing the results between these two gives us important insights on key factors (e.g., planet inclination) shaping the planetesimal + pebble accretion process.

3.1.1 Effect of inclinations on planetesimals growth

As mentioned in Sect. 2.3, the streaming instability produces in total $0.039\ M_{\oplus}$ planetesimals at $r_{\rm ice}$. In the mono-dispersed initial conditions, we assume that all planetesimals are $400\rm\ km$ in radius. ($6.7\times 10^{-5}\ \rm M_{\oplus}$ in mass). Therefore, $N=580$ planetesimals are generated in this compact ring belt ($\Delta r=\eta r$ in width). These bodies are injected into the simulation one by one after each $0.1$ orbital period. Their initial semimajor axes are uniformly distributed from $r_{\rm ice}-\Delta r/2$ to $r_{\rm ice}+\Delta r/2$. The initial eccentricities are assumed to follow a Rayleigh distribution, $p(e)=e/e_{0}\exp[-e^{2}/2e_{0}^{2}]$.

We conduct two sets of simulations. The first idealised set considers that planetesimals are in coplanar orbits ($i_{\rm p}=0$, $run\_md\_test$ in Table 1). However, this configuration ($i_{0}=0$) is an unphysical case. Realistically, although the initial inclinations are tiny, they are not zero due to the stochastic fluctuation driven by the disk turbulence (Ida & Lin, 2008; Gressel et al., 2011; Yang et al., 2012; Okuzumi & Ormel, 2013). The planetesimals would be lifted out of the midplane and acquire an inclination distribution.

In hydrodynamic simulations (Simon et al., 2016; Schäfer et al., 2017) the planetesimals are generated from pebble filaments, their random fluctuation is at least smaller compared to the size of the filament ($e_{0}$, $i_{0}\ll\eta$). Therefore, for the second set we assume that their inclinations also follow a Rayleigh distribution with $i_{0}=e_{0}/2$. Two different initial values for $e_{0}$ are numerically tested, $10^{-5}$ and $10^{-6}$. We find that since the orbit of planetesimals are readily excited, the results are insensitive to the initial values. In this non-coplanar configuration, three different simulations are performed with the randomlized initial semi-major axes and orbital phase angles (run_md_1 to run_md_3 in Table 1).

Figure 2 shows the mass growth of planetesimals when they are in co-planar orbits (panel a) and inclined orbits (panel b). We clearly find that mergers (red dots) occur more frequently when the orbits of the planetesimals are coplanar. In Fig. 2a the mass growth is dominated by planetesimal-planetesmial collision at the beginning (the sudden jump in mass of grey lines in Fig. 2). Since the pebble accretion rate is an increasing function of the planet mass, when the mass approaches $10^{-2}\ M_{\oplus}$, the growth is driven by pebble accretion (smooth growth in mass in Fig. 2). Only a few bodies survive after $1$ Myr, and the mass of the dominant body is $4\ M_{\oplus}$. Even though the accretion timescale in this configuration is artificially short, it clearly illustrates that the growth initially proceeds in a planetesimal accretion-dominated phase (phase A) and transitions to a rapid pebble accretion-dominated phase (phase B) when a massive body of \sim$$10^{-3}\ \rm M_{\oplus} to $10^{-2}\ \rm M_{\oplus}$ emerges. For the realistic case when planetesimals are on inclined orbits (Fig. 2b), however, much fewer mergers occur and the mass growth remains modest. The mass of the largest body remains $<10^{-2}\ \rm M_{\oplus}$, and there is no efficient pebble accretion at the end of the simulation.

3.1.2 Excitation of eccentricities and inclinations

Fig. 3 shows the time evolution of the root-mean-square (rms) of the planetesimals’ eccentricities (orange) and inclinations (cyan) in run_md_1. After a few orbits the eccentricities are quickly excited while the inclinations still remain low. This behavior occurs in the shear-dominated regime, when the velocity dispersion ($\delta v$) of planetesimals is smaller than their Hill velocity ($v_{\rm H}=R_{\rm H}\Omega_{\rm K}$) (Ida, 1990), where the mutual Hill radius $R_{\rm H}=(2m_{\rm p}/3M_{\star})^{1/3}a$. After a few hundred years their eccentricities get excited to values larger than the Hill velocity ($\delta v>v_{\rm H}$, dashed line in Fig. 3), the scattering transitions to the isotropic, dispersion-dominated regime. In that case the rms inclination excites to a value equals to half of the eccentricity, $i_{\rm rms}\simeq e_{\rm rms}/2$ (Ida et al., 1993).

Based on the disk model, we calculate the excitation and damping timescale of the velocity dispersion of the planetesimals. The timescale of the viscous stirring (eccentricity excitation) is given by (Ida, 1990; Kokubo & Ida, 2000),

[TABLE]

The damping timescale of eccentricity from the gas drag is (Adachi et al., 1976),

[TABLE]

This damping is inefficient for large planetesimals, and is a strongly increasing function of the semi-major axis of the planet $a$. Tidal interactions (type I torque) between an embryo and the gas also damps its eccentricity and inclination. The damping terms is expressed as (Artymowicz, 1993)

[TABLE]

The type I damping is inversely proportional to the mass of the embryo, and is negligible when the size of the body is smaller than $1000$ km ($10^{-3}\ \rm M_{\oplus}$). From the above equations we obtain $\tau_{\rm vs}\ll\tau_{\rm drag}<\tau_{\rm tidal}$. We conclude that both gas drag and type I damping are ineffective to circularize the orbit of planetesimals. Their eccentricities and inclinations keep increasing with time.

The next question is why pebble accretion is inefficient in the simulation in Fig. 2b? There are two reasons. First, fewer planetesimal collisions result in less massive planetesimals. Since the pebble accretion rate correlates with the mass of the body (Eqs. (9) and (10)), less massive bodies accrete pebbles more slowly. Second, pebble accretion is quenched due to a high eccentricity and inclination. Specifically, from Eq. (11) it follows that pebble accretion requires encounters to have a sufficiently low relative velocity

[TABLE]

Inserting $\Delta v\simeq ev_{\rm K}$ and $e\sim 10^{-2}$ we find that pebble accretion becomes only significant when the mass of the planet approaches $e^{3}\tau_{\rm s}M_{\star}\approx 10^{-2}M_{\oplus}$. But from the above discussion we know that gas drag and type I damping are ineffective to reduce the random motions of $400$-km-sized planetesimals. In addition, a high planetesimal inclination further reduces pebble accretion when it becomes larger than the aspect ratio of the pebble disk ($i_{\rm rms}>h_{\rm peb}$) (Levison et al., 2015). In Fig. 3 we find that this happens at $t\simeq 5\times 10^{5}$ yr. This means after that the planetesimals exceed the pebble layer during part of its orbits, reducing the amount of pebbles they can “eat”.

3.1.3 Spreading of the planetesimal belt

In Fig. 4 we plot the time evolution of the planetesimals’ semimajor axes in $\it run\_md\_1$. We define the full width of the planetesimal belt as $\Delta W=2\sqrt{\sum_{i}(a_{i}-r_{\rm ice})^{2}/N}$. The initial width is $\eta r_{\rm ice}$ (blue area). After their injection, planetesimals in this compact zone experience mutual gravitational interactions; scatterings excite their eccentricities and, in order to conserve the angular momentum, the range of the semi-major axes also expands. In Fig. 4 we find that the width increases steadily with time and becomes $0.05$ AU after $10^{5}$ yr.

This orbital expansion can be described by a diffusion process. Ohtsuki & Tanaka (2003) obtain an analytical expression of the viscosity due to the mutual gravitational scattering based on equal size planetesimals. Substituting their viscosity ($\nu$ in their Eq. (18)) into $\Delta W=\sqrt{\nu t}$, we have

[TABLE]

where the prefactor $C_{\rm w}^{3}=C_{\rm fit}[24{I(\beta)/\pi^{2}\ln(\Lambda^{2}+1)]}$, $I(\beta)=0.2$, $t$ is the time, $\Lambda=\bar{i}_{\rm rms}(\bar{e}_{\rm rms}^{2}+\bar{i}_{\rm rms}^{2})/3$, and $\bar{i}_{\rm rms}=ai_{\rm rms}/R_{\rm H}$ and $\bar{e}_{\rm rms}=ae_{\rm rms}/R_{\rm H}$. The numerical factor of $C_{\rm fit}=30$ is calibrated with our numerical simulations and is $\simeq 3$ times larger than Ohtsuki & Tanaka (2003)’s result. This formula is obtained under the assumption that planetesimals are in the dispersion-dominated accretion regime. We show in Fig. 4 that our analytical expression (Eq. (17)) agrees well with the simulation (black line). For fixed total mass $Nm_{\rm p}$ and $a$, $\Delta W$ increases with time ($\propto t^{1/3}$), the initial size of the planetesimal ($\propto R_{0}$), and decreases with the eccentricity ($e^{-2/3}$).

Since the ring belt expands over time, the surface density of the planetesimal also decreases. Therefore, the accretion of the planetesimals formed from a narrow ring will be longer than the classical /oligarchic accretion that assumes an ‘infinite’ width of the planetesimal disk (Kokubo & Ida, 2000).

It is clearly seen from Fig. 2b that the growth of planetesimals is slower than Fig. 2a when considering the non-zero inclinations. We find that the runs in Table 1 with different initial randomness have a similar growth trend. The masses of the most massive planetesimals from run_md_1 to run_md_3 are all far below an Earth mass.

In conclusion, when starting from a narrow ring belt of $400$-km-sized planetesimals, the growth of planetesimals is suppressed by the self-excitation of the planetesimals. Planetesimal accretion is slow and pebble accretion is inefficient. Formation of planeterary embryos, let alone super-Earth planets, is impossible in this scenario, even though $100\ \rm M_{\oplus}$ of pebbles drift through.

3.2 Poly-dispersed initial conditions with large planetesimals

In the previous section we found that when the sizes of planetesimals are all $400\ \rm km$, their growth is inefficient at the ice line. Then the question is: what is the realistic sizes of the planetesimals? In this section we consider a poly-dispersed population of planetesimals from the streaming instability simulations. In this case, planetesimals evolve into two-component mass distribution (one large embryo + small planetesimals). The large body is able to start rapid pebble accretion to eventually grow into a massive planet.

The initial mass function of planetesimals by the streaming instability has been investigated by numerous authors in detail. Simon et al. (2016, 2017) find that the initial size can be described by a single power law. Recently Schäfer et al. (2017) suggests that this initial mass distribution is better fitted by a power law plus an shallow exponential decay. The power law is for low and intermediate size planetesimals, while the exponential decay tail represents the high mass cutoff of planetesimals (Fig. 4 in Schäfer et al. (2017)). From Schäfer et al. (2017) simulations, the number fraction of planetesimals with $m_{\rm p}>m$ is given by their Eq.(19),

[TABLE]

Following this distribution, we adopt their parameters: $p=0.6$, $\beta=0.35$. The minimum mass and exponential characteristic mass are fitting parameters, which are given as $m_{\rm min}=10^{-6}\ \rm M_{\oplus}$ ($100$ km), m_{\rm e}$$=1.6\times 10^{-5}\ \rm M_{\oplus} ($250$ km), respectively. Figure 5 illustrates the mass distribution of planetesimals for all mass range (red + blue) and for the simulated mass branch (blue). There are in total $0.039\ \rm M_{\oplus}$ planetesimals. Here we simulate the most massive $470$ bodies from this population, corresponding to the two-thirds of the total mass (blue brach in Fig. 5). The simulated bodies span two order of magnitude in mass and a factor of $5$ in size, from $200$ km to $1000$ km. For simplicity, we neglect the lower one-thirds of small planetesimals. These bodies, if presented, would only modestly contribute to the planetesimal accretion and dynamical friction.

Other initial conditions are chosen to be the same as the simulation of equal $400$-km-sized planetesimals (Table 1). The mass growth of all planetesimals is shown in Fig. 6a while the semi-major axis of the largest body is shown in Fig. 6b. In Fig. 6c the eccentricity (solid) and inclination (dashed) evolution of the largest body is represented by the black line and the rms values of the planetesimals by the red line.

For the poly-dispersed population of planetesimals simulated here, most mass is in medium-sized planetesimals of $300$-$400$ km in size (Fig. 5). Energy equipartition leads to a configuration that the massive planetesimals have low random velocities while the low-mass planetesimals have high random velocities (Ida, 1990; Goldreich et al., 2004). Therefore, the random velocity of the massive planetesimal is lower than the rms velocity dispersion due to dynamical friction. At the beginning, the largest planetesimal (black line in Fig. 6a) can accrete more neighbouring planetesimals. Since its eccentricity and inclination remains low, once the mass is beyond a few $10^{-3}\ \rm M_{\oplus}$, the largest planetesimals will begin to accrete pebbles very efficiently. The mass of the largest planetesimal then increases further, rendering pebble accretion more efficient compared to lower-mass planetesimals. Therefore, this massive embryo dominates the subsequent growth and dynamical evolution of the system.

Rather than stop the simulation at $1$ Myr, we continue it until the most massive planet migrated into the inner $0.1$ AU ( the assumed location of the inner disk edge) in this illustrated run (run_pd_1). We find that the planet reaches the inner disk edge at $1.2$ Myr when the mass is $3.9\ \rm M_{\oplus}$. The second largest body only grow its mass rapidly when the largest one migrated out of the ring. In the end of the simulation the second largest body is still one order of magnitude less massive than the largest one.

To conclude, the presence of $\sim 10^{3}$ km planetesimals in a sea of predominantly a few hundreds km bodies will significantly speed up planet formation. This poly-dispersed population results in a low velocity dispersion of the massive body. Both effects (high $m_{\rm p}$ and low $\Delta v$) facilitate the transition to the pebble accretion dominated phase (phase B) and promote the formation of a super-Earth planet after $\simeq 1$ Myr.

3.3 Two-component initial conditions as a result of runaway growth of small planetesimals

As demonstrated above, a poly-dispersed population of planetesimals containing initially a few very large bodies can form planets rapidly. Such an initial planetesimal mass distribution would be a consequence of forming planetesimals directly from the streaming instability. Generally, a lower pebble surface density yields smaller planetesimals in streaming instability simulations (Johansen et al., 2015; Simon et al., 2016). For the single-site planetesimal formation scenario, the pebble density is locally enhanced at a particular disk location (e.g., the ice line). But how strong the enhancement is determined by complicated physical processes such as water vapor diffusion and condensation (Schoonenberg & Ormel, 2017). Taken these factors and uncertainties into account, we also consider the case of an initially small planetesimal size. Nevertheless, we argue in Sect. 3.3.1 that in this circumstance runaway growth of planetesimals proceeds fast, producing the desired two-component mass distribution that promotes the formation of massive planets.

3.3.1 Runaway growth in a narrow planetesimal belt

The characteristic runaway planetesimal accretion timescale is (Ormel et al., 2010)

[TABLE]

where C_{\rm rg}$$\simeq$$0.1 is a numerical-corrected prefactor from Ormel et al. (2010), and $\rho_{\bullet}=1.5\ \rm g\,cm^{-3}$ is the internal density of the planetesimal. Small planetesimals have three advantages in terms of mass growth. First, based on Eq. (19) the runaway growth timescale is proportional to the initial size of the planetesimals ($R_{0}$). It takes less time to form a runaway body (before the onset of the oligarchic growth) when starting with smaller planetesimals. Second, the eccentricity and inclination excitation by self-stirring is less severe for smaller planetesimals (Equation (13)), which also boosts the planetesimal and pebble accretion. Third, as shown by in Sect. 3.1.3, the spreading of their semimajor axes increases with the size of the planetesimal ($R_{0}$ $\propto$ $R_{\rm H}$ in Eq. (17)). For the same total planetesimal mass, smaller planetesimals expands their orbits less, resulting in a higher surface density. Therefore, the subsequent accretion also becomes faster.

Accounting for these effects, we assume that the planetesimals are all $100\ \rm km$ in size ($\simeq 10^{-6}\ \rm M_{\oplus}$ in mass). These planetesimals would undergo a faster runaway growth and less orbital expansion compared to $400$-km- sized bodies. Since the runaway growth will end up with a steep mass distribution ($dN/dm_{\rm p}\propto m_{\rm p}^{-2.5}$ (Wetherill & Stewart, 1993; Kokubo & Ida, 1996, 2000; Ormel et al., 2010; Morishima et al., 2008), most of the mass remains in small planetesimals. As discussed in Sect. 3.2, the inclination and eccentricity of the large body could therefore be damped by small planetesimals through dynamical friction. The planet can accrete more pebbles when the orbit remains nearly coplanar and circular.

The key difference between the classical runaway and oligarchic growth scenario and our scenario is that we here consider planetesimal formation at a single location. In this study we propose a clean and simplistic physical picture. The hypothesis is that after the planetesimal runaway growth phase, the system can be well-described by a two components: a single big embryo of radius $\sim 10^{3}$ km and a swarm of $100$ km size planetesimals, which dominate the total mass. We reasonably neglect bodies in between these two are dynamically insignificant and cannot compete with the runaway body in terms of growth (also shown in Fig. 6a). Altogether, the situation resembles the classical runaway/oligarchy transition, except that, in our case, there is only one single “oligarch”.

Why only one big embryo instead of multiple oligarchs? The transition from the runaway to the oligarchic growth takes place when the viscous stirring timescale $\tau_{\rm vs}$ (Eq. (13)) equals the runaway timescale $\tau_{\rm rg}$ (Eq. (19)) 222This criterion is different from Kokubo & Ida (1998). See discussion in Ormel et al. (2010). . When $\tau_{\rm vs}>\tau_{\rm rg}$, the stirring of the random velocities is smaller compared to the accretion, which is the key feature of runaway growth. When $\tau_{\rm vs}<\tau_{\rm rg}$, the eccentricity growth is faster than the accretion, and the growth transitions to the oligarchic regime. The transition size of the planet between two regimes is given by Eq. (11) of Ormel et al. (2010),

[TABLE]

Substituting Eq. (19) into Eq. (17), we obtain that when the massive body grows into a transition radius, the planetesimal belt width is

[TABLE]

We note that $R_{\rm H}$ is the Hill radius of the planetesimal with size of $R_{0}$. From $\Sigma_{\rm p}=M_{\rm tot,pl}/2\pi a\Delta W$ and and Eq. (21), we calculate that the belt width extends to $\Delta W=0.029$ AU and the planetesimal surface density reduces to $\Sigma_{\rm p}=2\rm\ g\,cm^{-2}$ after the runaway growth ($\sim 10^{5}$ yr). Substituting the above valuses into Eq. (20), we obtain,

[TABLE]

We therefore conclude that the feeding zone of this embryo ($\simeq 10\ R_{\rm H}$, Kokubo & Ida (1998)) is similar to the width of the planetesimal belt, supporting one single embryo assumption.

3.3.2 The superparticle approach

Since the total planetesimal mass is $0.039\ \rm M_{\oplus}$ and small planetesimals dominate the mass, there are in total $35000$ 100-km-sized planetesimals (\simeq$$10^{-6}\ \rm M_{\oplus}) in this case. It is prohibitively computationally to simulate the interactions among all these bodies with an N-body scheme. Therefore, we adopt the superparticle approach to mimic the dynamics of small bodies (see Levison et al. (2015); Raymond et al. (2016) and references therein), in which $N_{\rm sp}$ small planetesimals are clustered as one superparticle. The superparticle still feels the aerodynamical gas drag as if it was a single $100$-km-sized planetesimal, but the mass of the particle is $N_{\rm sp}$ times higher ($m_{\rm sp}=N_{\rm sp}m_{0}$). These superparticles gravitationally interact with the embryo but not with each other. This is based on the fact that the collision timescale among small planetesimals is much longer than the embryo-planetesimal collision. In principle, the approach requires $m_{0}\ll m_{\rm sp}\ll m_{\rm p}$, where $m_{\rm p}$ is the embryo’s mass. The latter inequality is needed to ensure that dynamical friction operates correctly. The simulation would be very time-consuming for a too small $N_{\rm sp}$, while the dynamical evolution may become artificially stochastic for a too large $N_{\rm sp}$ (low number of superparticles).

We here assume that the velocity dispersion of small planetesimals is excited by the planetesimals themselves, and that it is close to their escape velocity ($\delta v\simeq v_{\rm esc}=\sqrt{2Gm_{0}/R_{0}}$), i.e., $e\simeq e_{\rm esc}=v_{\rm esc}/v_{\rm K}$. The eccentricities and inclinations are adopted to follow the Rayleigh distributions, where $e_{0}=2i_{0}=e_{\rm esc}$. Their semi-major axes are randomly initialized from $r_{\rm ice}(1-\eta/2)$ to $r_{\rm ice}(1+\eta/2)$. In addition, since pebble accretion is extremely inefficient for $100$ km size planetesimals (the inequality (16) is far from satisfied), it is reasonable to neglect pebble accretion of the superparticles. The embryo is initially placed at $r_{\rm ice}$.

3.3.3 Results

We have conducted simulations with different masses of the superparticle, $N_{\rm sp}=25$ and $50$, respectively. The results are in good agreement with each other (see Appendix for details). The results presented in Sects. 3.3.3 and 3.3.4 are based on the $N_{\rm sp}=50$ run.

The disk and pebble parameters in this (run_sp) are taken to be identical to * run_pd* in Sect. 3.2. The mass of the embryo is chosen to be the transition value ($m_{\rm p}=10^{-3}\ \rm M_{\oplus}$ when $R_{\rm p}=R_{\rm rg/oli}$ in Eq. (22)). Three individual simulations (run_sp_1 to run_sp_3) are conducted, varying only the initial positions and orbital phase angles of the planetesimals. In addition, we also conduct a separate simulation for the growth of a single embryo without small planetesimals. This comparison enables us to isolate the effect of N-body dynamics (planetesimal accretion). Simulations are terminated when $t=1$ Myr.

Figure 7 shows the mass growth and orbital evolution of the embryos in panel (a) and (b), respectively. For the superparticle approach, the thick red lines represent the mean values averaging three runs (e.g., ${\bar{m}}(t)=\sum_{i=1}^{3}m_{\rm i}(t)/3$), whereas the light red areas represent the spreading among these runs. It can be seen that the difference in mass and semi-major axis among the three simulations is small. The black dashed line represents the growth of a single embryo (assuming it is on a circular and coplanar orbit) purely by pebble accretion. We find that when starting with the transition size embryo, the planetesimal accretion are already modest and the mass growth is mainly driven by pebble accretion.

When the mass of the embryo is beyond $0.1\ \rm M_{\oplus}$, the type I migration becomes important. The planet would take $\sim 0.5$ Myr to migrate into the inner region of the disk. Eventually, a $1.3\ \rm M_{\oplus}$ planet forms within $1$ Myr.

3.3.4 Timescale analysis

In this section we analyse the growth timescale of the embryo. The pebble accretion timescale in the 2D regime is given by

[TABLE]

where we use Eq. (9) and $\Delta v$ is assumed to be dominated by the Keplerian shear since the eccentricity of the embryo is insignificant due to dynamical friction. From Eq. (11) the above $f_{\mathrm{set}}$ evaluates as

[TABLE]

When the planet is $\gtrsim 10^{-2}\ \rm M_{\oplus}$, $f_{\rm set}\simeq 1$. From Eq. (24), we clearly see that $f_{\rm set}$ becomes smaller when a planet is less massive or their pebble-planet relative velocity is higher.

Similarly, in the 3D limit (Eq. (10)), the pebble accretion timescale becomes

[TABLE]

From the planet mass dependence on Eqs. (23) and (25), we know that the pebble accretion is in $3$D when the planet mass is low, and it transitions to $2$D accretion when the planet becomes more massive. The transition mass between $2$D and $3$D is $\simeq 0.05\ \rm M_{\oplus}$ for the adopted parameters.

In order to numerically capture the transition from the planetesimal dominated accretion to the pebble dominate accretion (phase A to phase B), we conduct the superparticle approach simulations starting with a less massive embryo ($m_{\rm p}=3\times 10^{-4}\ \rm M_{\oplus}$) here. The other conditions are the same as $run\_{sp}$ in Table 1. Again, three individual superparticle approach simulations are performed with randomly varied initial positions and orbital phase angles of the planetesimals.

These timescale are shown in Fig. 8, where the grey, blue, light blue and magenta lines represent growth timescales obtained from the simulation, the pebble accretion in $3$D (with and without $f_{\rm set}$) and in $2$D, respectively. We can identify the following different growth stages.

When the mass of the embryo is lower than $\sim 10^{-3}\ \rm M_{\oplus}$, the mass growth is dominated by planetesimal accretion (phase A). In this phase the growth time is relatively long, and decreases with the growing mass of the big body. Since the timescale calculation is affected by the stochasticity of the collisions, the average value among three different runs (grey line) is given in Fig. 8. The pebble accretion is inefficient in this phase because of the low embryo’s mass (small $m_{\rm p}$ in Eq. (25)) and hence the weak setting accretion (small $f_{\rm set}^{2}$ in Eq. (25)). When the mass of the embryo becomes larger, $3$D pebble accretion is more efficient than planetesimal accretion (phase B). The transition mass from the planetesimal accretion dominated regime to the pebble accretion dominated regime is close to $10^{-3}\ \rm M_{\oplus}$. With increasing mass $f_{\rm set}$ increases. When $f_{\mathrm{set}}$ approaches unity, the growth timescale becomes independent of the planet mass (see Eq. (25) and the light blue in Fig. 8). Finally, growth enters the $2$D pebble accretion regime when the mass approaches $0.1\ \rm M_{\oplus}$. In that case growth timescale starts to increase, $\tau_{\rm PA,2D}\propto m_{\rm p}^{1/3}$ (magenta line in Fig. 8), since $\varepsilon_{\rm PA,2D}$ no longer scales linearly with the planet mass.

3.4 Discussion

Comparing the cases of initial $100$-km-sized and $400$-km-sized planetesimals, our result suggests that starting with smaller planetesimals is more optimal for forming planets. The reason is as follows: 1) the runaway body’s size is insensitive to the initial planetesimal size ($R_{\rm rg/oli}\propto R_{0}^{1/7}$, Eq. (22)), and 2) since the orbital expansion is less significant for smaller planetesimals ($\Delta W\propto R_{0}$ from Equation (21)) and $\Sigma_{\rm p}\propto\Delta W^{-1}$, the runaway accretion timescale is strongly correlated with the initial size ( $t_{\rm rg}\propto R_{0}^{2}$ from Eq. (19)). For the mono-dispersed population of $400$-km-sized planetesimals, planet growth fails mainly because it takes too long to evolve into a two-component mass distribution. However, in the case of $100$-km-sized planetesimals, this configuration can be realized within a fraction of the disk lifetime. In that case, the eccentricity and inclination of the large embryo remains low due to the dynamical friction of small planetesimals. It is therefore possible for the largest embryo to accrete pebbles efficiently and grow into an Earth-mass planet. As a result, massive planets could form from planetesimals with a smaller initial size.

The above discussed preferable size of the planetesimal ($R_{0}$) to grow protoplanets is dependent on the formation location of planetesimals. When the planetesimals form further out, from Eq. (19) the runaway planetesimal growth timescale increases. On the other hand, the starting size of the embryo for the efficient pebble accretion also increases with the distance (Visser & Ormel, 2016). This is because $f_{\rm set}$ becomes smaller when $\eta$ increases with $r$. Therefore, both two effects further limit the growth of the protoplanets from a large size planetesimal (e.g., $400$ km) when the planetesimal formation occurs at a large distance.

Our findings are similar to Levison et al. (2015). They report that large core formation by pebble accretion is only feasible after a period of classical planetesimal runaway growth. The system then ends up with a few large embryos and a swarm of small planetesimals. The embryos stir the random velocities of all small planetesimals, preventing any further growth by pebble accretion. In their case, planetesimals were born over the entire disk region and assumed to follow the same power-law as the disk gas, which is based on the classical planet formation scenario (e.g., Ida & Lin (2004) and Mordasini et al. (2009)). The key difference with this work is that we consider the planetesimals to form at a single site disk location with a narrow radial width. Nevertheless, the ‘common ground’ between both works is that a two-component mass distribution for planetesimals is essential to promote further growth from planetesimals into protoplanets.

4 Growth and migration in the pebble accretion dominated regime

In this section, we particularly focus on growth in the pebble accretion dominated regime (phase B). We start the mass of the embryo from $10^{-3}\ \rm M_{\oplus}$. As illustrated in Fig. 7, the pebble accretion dominates the growth and the planetesimals’ contribution can be neglected after $m_{\rm p}\gtrsim 10^{-3}\ \rm M_{\oplus}$. The relative velocity of the massive embryo is very low due to the dynamical friction by small planetesimals and later on type I tidal damping when it grows larger. It is therefore justified to assume the embryo is on a circular and coplanar orbit. We conduct simulations for a single embryo and neglect all planetesimals, investigating the role of various disk and pebble parameters as shown in Table 2. Simulations are terminated when the planet has migrated inside of 0.1 AU.

In Table 2, the fiducial run (run_fid) is adopted to be the same disk and pebble values as in the previous sections. From run_nmig to run_lpeb, only one parameter is varied compared to the fiducial run_fid. For instance, in run_nmig we assume the planet is at the zero-torque location, where $\dot{a}/a=0$. In this case, the planet does not undergo type I migration. Therefore, the planet accretes materials in-situ at the ice line. For a comparison between run_nmig and run_fid, we will gain a knowledge of the effect of migration on the growth of the planet. For this purpose, we stop run_nmig at the same time when the planet in run_fid migrates inside of $0.1$ AU. For a comparison between run_fid and the other individual runs (run_alpha, run_tau, run_hpeb, run_lpeb) , we can understand the effect of disk turbulence, pebble size and pebble mass flux on the planet growth.

We find in Fig. 9 that in a low turbulent disk (run_alpha, magenta) and in a disk with high pebble flux (run_hpeb, thick red), the mass growth is faster than the fiducial disk (red). In these two circumstances the planets reach $4.1\ \rm M_{\oplus}$ and $6.8\ \rm M_{\oplus}$, respectively, when they migrate inside of $0.1$ AU within $1$ Myr. We find the mass growth in run_nmig is very similar to the fiducial run but slightly less efficient. At $t=1.2$ Myr, the final planet mass is $4.0\ \rm M_{\oplus}$ in run_fid while in run_nmig the planet always stays at the ice line and attains $3.8\ \rm M_{\oplus}$. It is also clearly seen that when the pebble flux is lower (run_lpeb, light red), or the Stokes number is lower (run_tau, orange), the mass growth is slower than the fiducial run. The protoplanet with $m_{\rm p}\gtrsim 1\ \rm M_{\oplus}$ does not form within $1$ Myr from these two configurations. In run_tau the planet attains $4.8\ \rm M_{\oplus}$ when it arrives at the inner edge of the disk at $t=1.5$ Myr, while in run_lpeb it grows into a $2.3\ \rm M_{\oplus}$ planet at $t=2.1$ Myr.

The planet grows much faster in a high pebble flux disk. The effect of pebble mass flux on planet growth is intuitive since pebble accretion benefits from a high pebble flux (a massive pebble disk). When the planet reaches $6.8\ \rm M_{\oplus}$ and migrates into the inner disk cavity at $t=0.7$ Myr in run_hpeb, the planet mass is only $0.02\ \rm M_{\oplus}$ in run_lpeb. We find a super linear correlation between the planet mass ($m_{\rm p}$) and the integrated pebble flux ($\dot{M}_{\mathrm{peb}}t$). A factor of $4$ change in $\dot{M}_{\rm peb}t$ results in more than two order of magnitude growth in $M_{\rm p}$. This is due to the fact that pebble accretion efficiency ($\varepsilon_{\rm PA}$) increases with $m_{\rm p}$ in both $2$D and $3$D regime. The planet can initially accrete more pebbles in a high pebble flux disk. It becomes more massive due to this faster growth and would further accrete even a higher fraction of pebbles. This positive feedback promotes the rapid growth of the planet in a high pebble flux disk.

When the disk turbulence is lower, the pebble scale height becomes smaller (Eq. (8)). The transition mass from $3$D to $2$D pebble accretion will also become lower. Therefore, in run_alpha the embryo starts efficient $2$D pebble accretion early than the fiducial run. Eventually, it takes less time to grow into a super-Earth planet in a lower turbulent disk.

The Stokes number has an effect opposite to the turbulent $\alpha_{\rm t}$ in terms of the pebble scale height. Lower Stokes number pebbles mean they are more tightly coupled to the gas, and the pebble scale height becomes larger. Therefore, from Eqs. (23) and (25) transition mass from $3$D to $2$D in run_tau is $0.15\ \rm M_{\oplus}$, three times higher than the fiducial run. In this case starting with $10^{-3}\ \rm M_{\oplus}$, the planet grows a significant fraction of its mass in the slow, $3$D pebble accretion regime. On the other hand, when the planet enters $2$D pebble accretion, the accretion is faster when the Stokes number is lower (Eq. (23)). Balancing these two effects, we find that the planet growth is slightly slower in run_tau compared to the fiducial run.

The effect of pebble mass flux on the planet growth is super linear. A high pebble disk mass benefits the formation of a massive planet. A less turbulent disk and a large Stokes number pebbles also promotes pebble accretion and the formation of a massive planet.

5 Summary

Streaming instability is an important mechanism to convert pebbles into planetesimals. It occurs at the location in the disk where the pebble density is locally enhanced (e.g., the ice line). In this work, we have focused on the growth of the planetesimals after they have formed by this mechanism at a single site disk location.

This annular (ring) planetesimal formation scenario differs from the classical planetesimal accretion scenario from two aspects. First, streaming instability generates planetesimals in a narrow ring at a specific location. Most of the mass is in the largest planetesimals of a few hundred km in size. In contrast, in the classical scenario planetesimals form everywhere in the disk, e.g., typically following the surface density distribution of the gas. In simulations their planetesimal surface density remains constant by suppling material from the neighbouring region (Kokubo & Ida, 2000). This means that the orbital spreading is unimportant, in contrast to the ring formation scenario. Second, the planetesimals in our scenario grow their mass by accreting surrounding planetesimals and inwardly drifting pebbles. The mass at which the planet transitions from accreting predominantly planetesimals to pebbles occurs at $\simeq 10^{-3}\ \rm M_{\oplus}$.

We have modified the Mercury N-body code to perform simulations of the mass growth and orbital evolution of these planetesimals (Sect. 2). The code includes the effects of gravitational interactions and collisions among planetsplanetesimals, planet-disk interactions (gas drag and type I torque), pebble accretion based on the calculation of Liu & Ormel (2018) and Ormel & Liu (2018) that accounts for the disk parameters ($\tau_{\rm s}$, $\eta$, $\alpha$ and $\dot{M}_{\rm peb}$) and planet properties ($m_{\rm p},a,e,i$). Simulations with different initial planetesimal sizes and disk parameters are investigated.

The key findings of this study are the following:

Protoplanets cannot emerge from a mono-dispersed population of $400$ km size planetesimals, fuelled by a $100\ \rm M_{\oplus}$ reservoir of pebbles in the outer disk. Although the initial eccentricities and inclinations are very tiny, they soon get excited through gravitational scatterings. Mechanisms such as gas drag, type I damping are not sufficient enough to to damp their random velocities. Both planetesimal and pebble accretion are strongly suppressed when inclinations and eccentricities of planetesimals become moderate. In this circumstance, the growth of the planetesimals is mainly in the slow planetesimal accretion dominated phase (Sect. 3.1). 2. 2.

Protoplanets can form when streaming instability has in addition spawned a population of larger planetesimals. The largest body grows by planetesimal accretion. Soon after it approaches $10^{-2}\ \rm M_{\oplus}$, the growth enters the rapid pebble accretion dominated regime. During this time the random velocity of the largest body remains low through the dynamical friction of small planetesimals. Finally a super-Earth planet can form within $1$ Myr (Sect. 3.2). 3. 3.

Alternatively, protoplanets also form out of their birth ring when the initial size of the planetesimals is small (e.g., $100$ km). These small planetesimals are expected to undergo a runaway planetesimal accretion to form a massive embryo rapidly. In this way, the two-component mass distribution is also achieved. We find that an Earth mass planet can form as well (Sect. 3.3) 4. 4.

Planets grow larger when the pebble mass flux is higher, the disk is less turbulent, the Stokes number of pebbles are larger. In particular, the growth of the planet mass increases super linearly with the disk pebble mass flux (Sect. 4).

Planetesimal accretion and pebble accretion are not two isolated processes in planet formation. From the streaming instability point of view, pebbles are converted into planetesimals, whereafter these planetesimals accrete nearby planetesimals and pebbles at the same time. The total amount of solids in disks is either in small pebbles, or in large planets/planetesimals. Therefore, a certain level of competition exist between planetesimal formation and planet growth by pebble accretion.

In addition, this work only considered the case of a single burst of planetesimal formation by streaming instability. In reality, even at the ice line the streaming instability may be triggered multiple times (episodic bursts) during the gas disk lifetime, because the planets migrate out of their birth ring. Protoplanets then emerge sequentially and a chain of multiple planets forms, as envisioned by Ormel et al. (2017).

To address these issues, a global disk model is needed. Using a novel Lagrangian approach, such a model has just been developed (Schoonenberg et al., 2018). Because of the flexibility of the Lagrangian (particle-oriented) model, it is straightforward to couple it to the N-body model presented in this paper. Such a model can then, potentially, simulate planet formation in its entirety – starting from dust coagulation and ending with a planetary architecture. It can be applied to model “complete” (as far as one can tell) planetary systems, such as those discovered around TRAPPIST-1 (Gillon et al., 2017).

Appendix A Convergence test for superparticle simulations

We show the simulations with different mass of the superparticle ($m_{\rm sp}=N_{\rm sp}m_{\rm 0}$) for $run\_sp$ starting with $m_{\rm p}=3\times 10^{-4}\ \rm M_{\oplus}$ in Sect. 3.3.4. Three simulations with randomlized initial conditions are performed for each set of $N_{\rm sp}$. The results are shown in Fig. 10 where red, green represent $N_{\rm sp}=25$ and $50$, respectively. The thick line represents the mean value averaged from three individual simulations, for instance, ${\bar{m}}(t)=\sum_{i=1}^{3}m_{\rm i}(t)/3$ whereas the light region marks the range between the minimum and maximum values from the three simulations.

Fig. 10 shows the mass, semi-major axis evolution of the embryo and eccentricities of both planetesimals and the embryo. We find that the mass and semi-major axis converge quite well for the above two $N_{\rm sp}$. In Fig. 10a the difference between the mean mass for the above $N_{\rm sp}$s is smaller compared to the spreading among individual runs (stochastic N-body effects). The excitation of planetesimals by the presence of the embryo is also agreed with each other in in Fig. 10c.

To summarize, results from the above three tested $N_{\rm sp}$ are in general agreement with each other.

Acknowledgements.

We thank Chao-Chin Yang, Lixin Li, Carsten Dominik, Melvyn Davies for useful discussions, and Bertram Bitsch, Djoeke Schoonenberg for proofreading the manuscript. We also thank the anonymous referee for their insightful suggestions and comments. Beibei Liu thanks the support of the Netherlands Organization for Scientific Research (NWO; VIDI project 639.042.422), the European Research Council (ERC Consolidator Grant 724687-PLANETESYS) and the Swedish Walter Gyllenberg Foundation. Chris Ormel is supported by the Netherlands Organization for Scientific Research (NWO; VIDI project 639.042.422). Anders Johansen is funded by the Swedish Research Council (grant 2014-5775), the Knut and Alice Wallenberg Foundation (grants 2012.0150, 2014.0017) and the European Research Council (ERC Consolidator Grant 724687-PLANETESYS). The computations are performed on resources provided by the Swedish Infrastructure for Computing (SNIC) at the LUNARC-Centre in Lund.