Growth after the streaming instability: from planetesimal accretion to pebble accretion
Beibei Liu Chris W. Ormel Anders Johansen



From planetesimal accretion to pebble accretion



Lund Observatory, Department of Astronomy and Theoretical Physics, Lund University, Box 43, 22100 Lund, Sweden.  
Anton Pannekoek Institute (API), University of Amsterdam, Science Park 904,1090GE Amsterdam, The Netherlands

bbliu@astro.lu.se,c.w.ormel@uva.nl,anders@astro.lu.se
 





Streaming instability is a key mechanism in planet formation, clustering pebbles into planetesimals with the help of self-gravity. 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 into protoplanets 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 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.



We modify the  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.


   
 The mono-disperse population of 400 km size planetesimals cannot form ≳ 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. 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. 

 

    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.
  





Growth after the streaming instability
     Beibei Liu 1,2, Chris W. Ormel2, Anders Johansen 1
    December 30, 2023
=======================================================







§ INTRODUCTION
 




In protoplanetary disks,  micron-sized dust grains  coagulate into pebbles of mm-cm sizes  <cit.>. But  further growth is suppressed by bouncing or fragmentation due to the increasing compactification in collisions <cit.>. In addition, these pebbles also drift too fast compared to their growth <cit.> so that the particles cannot cross the meter size barrier even if they would stick perfectly <cit.>, unless the pebbles can remain fluffy  during the growth <cit.>. The subsequent growth of these pebbles to planetesimals is still not well understood in planet formation theory (see  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 <cit.>.
 The characteristic  size of these planetesimals  approximate a few hundred kilometres  <cit.>.

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 <cit.> 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 (dm_ p/ dt ∝ m_ p^γ with γ>1), which is termed `runaway growth' <cit.>. 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 dm_ p/ dt ∝ m_ p^γ with γ<1 <cit.>. This phase is characterized by a decreasing mass ratio of  two adjacent massive runaway bodies <cit.>. 

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 <cit.>. This is known as pebble accretion  <cit.>. 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 <cit.>; and second,  a large flux of pebbles grow and drift inward from the outer regions of disks <cit.>. Pebble accretion can be classified into 2D/3D regimes <cit.>. When the pebble accretion radius is larger than the pebble scale height, the accretion is in the 2D  regime, where dm_ p/ dt ∝ m_ p^2/3 <cit.>. 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 3D regime is reduced compared to 2D, and d m_ p/ dt ∝ m_ p <cit.>.  

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 ε_ PA <cit.>, defined as the number of pebbles that are accreted divided by the total number of pebbles that the disk supplies. 
Recently we have computed  ε_ PA  under general circumstances <cit.>.
For instance, when the eccentricity and inclination of the planet become high, ε_ 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 (, the H_2O ice line).  <cit.> 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 method). Furthermore, <cit.> 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, <cit.> concluded that their born size should be large (≳ 100 km), while <cit.> argued that planetesimals starting from sub-km-sized still cannot be ruled out. 
<cit.> 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 ≳ 100km. In the context of combined planetesimal and pebble accretion,  <cit.> 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  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 method, we outline our model and the implementation of the N-body code.  
In threesizes three initial size distributions of planetesimals are investigated, including a mono-dispersed population in mono-dispersed, a poly-dispersed population in poly-dispersed, and a single  runaway body plus a swarm of small planetesimals in two-component. In parameter we  investigate the influence of different  parameters in the pebble accretion dominated growth phase (phase B). The  key results are summarized  in summary.




§ 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 <cit.>.  For instance,  <cit.> 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_ ice =2.7    AU based on the disk model in disk) in this paper. But the following results and applications  can be scaled to other locations where the streaming instability condition is realized, , the inner edge of the  zone <cit.> or a distant location due to the FUV photoevaporation <cit.>. 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 carton.
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 (Δ r = η r, see carton and further discussion in planetesimal initial condition). 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). 



 §.§ Disk model


The surface density and the aspect ratio of the gas disk are assumed to be   
 
    Σ_ gas =  Σ_ gas0(  r/1 )^-1,
   
 
    h_ gas =  h_ gas0(  r/1 )^1/4,

where Σ_ gas0 and h_ 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 <cit.>.
Here Σ_ gas0 = 500  g cm^-2 and h_ gas0 = 0.033  <cit.> are adopted as the default values in this paper. 
Therefore, the disk temperature T and  η are 
 
    T =  μ h_ gas^2  G M_⋆/r R_ g  = T_0 (  r/1 )^-1/2,

and 
 
    η =  -h_ gas^2/2∂ log P_ gas/∂ log r = η_0 (  r/1 )^1/2,

where G is the gravitational constant, R_ g = 8.31 ×10^7 erg   mol^-1  K^-1 is the gas constant, M_⋆ is the stellar mass, μ=2.34 is the mean molecular weight and P_ gas is the gas pressure in the disk.  In the above equations  T_0 = 280  K and  η_0=1.5 × 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 (Σ_ g∝ r^-1) is more consistent with disk observations <cit.>. The fiducial value of Σ_ gas0 is chosen to be the same as <cit.>, and the pebble flux is calibrated accordingly in pebble. 
Based on the above disk model, the ice line (T∼ 170  K) is located at  r_ ice = 2.7 AU.



 §.§ 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 <cit.> 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 <cit.> 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 gas.  In pebble, we implement the pebble accretion prescriptions by <cit.> and <cit.>, 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.



  §.§.§ Gas drag and type I migration torque


Small embryos and planetesimals experience  the aerodynamic gas drag <cit.>, 

    a_ drag =  - (  3 C_ Dρ_ gas/8  R_ pρ_∙) v_ relv_ rel,

where the drag coefficient C_ D =0.5, v_ rel  is the relative velocity between the planetesimal and the gas, v_rel = v-v_gas, where v_gas equals v_ K(1-η) in the azimuthal direction,  is the Keplerian velocity, ρ_ gas is the local gas density, ρ_∙ and   R_ p are the internal density (assumed to be 1.5  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, <cit.>). The characteristic migration timescale  for a planet on a circular orbit is <cit.>
 
    t_ m≃     0.5 ( M_⋆/ m_ p)  ( M_⋆/Σ_ gas a_ p^2  )  (  h_ gas^2 /Ω_ K)
      ≃ 5×10^5( m_ p/ 1 )^-1
       ( Σ_ gas0/500  g cm^-1)^-1( h_ gas0/3.3 × 10^-2)^2( a_ p/2.7  )   yr,

where m_ p and a_ p are the mass and the semimajor axis of the planet, respectively, and Ω_ K is the Keplerian angular velocity at the planet location.



The accelerations acting on planets due to type I migration torque are 
   
    a_ m = -v/t_ m,  a_ e = -2 (v·r) r/ r^2 t_ e,  a_ i = - v_ z/ t_ i,

 where r = (x,y,z), v = (v_ x,v_ y,v_ z)  are the position and the velocity vectors of the planet. In the above expression, t_ m, t_ e and t_ i are the characteristic type I migration, eccentricity and inclination damping timescales from Eqs. (13), (11) and (12) of <cit.>. 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 <cit.>, the dynamical torque <cit.>, and the heating torque from planet gas accretion <cit.> 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 parameter (see Table 2). In that case we neglect the semimajor axis damping  (set a_ m =0) but still consider the eccentricity and inclination damping in amig 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.
 
 


  §.§.§ Pebble accretion
 

The pebble-sized particles drift inwards across the protoplanetary disk. The radial drift velocity is v_ r = -2 η v_ Kτ_ s/(1 + τ_ s^2) <cit.> where  is the dimensionless stopping time (Stokes number). The drift speed is determined by the aerodynamical size of the pebble  and  η.  Based on <cit.> and <cit.>, 
the pebble mass flux (Ṁ_ peb = 2π r  v_ rΣ_ 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 Ṁ_ peb= 100   M_⊕/Myr (consistent  with Eq. (14) of <cit.>) and neglect the time dependence for simplicity. A lower Ṁ_ 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 (ρ_ peb/ρ_ gas >1) and the streaming instability will be triggered over the entirely disk. The above 
inequality can be written as   Σ_ peb/H_ peb > Σ_ gas/H_ gas where the scale height of the pebble is given by <cit.>
 
    H_ peb = √(α_ t/α_ t + ) H_ gas.
  
The pebble flux therefore needs to be smaller than 2π r  v_ rΣ_ gas H_ peb/H_ gas. Our default value is well below this limit. 
We also explore two additional pebble fluxes    Ṁ_ peb = 200  Myr^-1 and 50  Myr^-1  in parameter.


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 <cit.>, which is correlated with the gas disk accretion rates <cit.>.  For a typical T Taurs star with a disk accretion rate of 10^-8  M_⊙  yr^-1, the dust mass is ∼ 100 (Fig.1 of <cit.>), consistent with our adopted fiducial value.


 In Hpeb α_ t is the coefficient of turbulent gas diffusivity, which is different from the concept of turbulent viscosity α_ν.  We nevertheless note that the above two quantities are approximated the same when the turbulence is driven by magnetorotational instability <cit.>.  The value of α_ t can be constrained from the molecule line broadening measurements <cit.> and the level of dust settling <cit.>.
We adopt a  fiducial value of α_ t  =10^-3 and test the influence of a lower turbulent disk  (α_ t  =10^-4) in parameter.

For pebbels we adopt a fiducial aerodynamical size of =0.1. This chosen value is consistent with the sophisticated dust coagulation model <cit.> and disk observations <cit.>.  A lower =0.03 is also explored in parameter.


 

A fraction of pebbles will be accreted onto a planetesimal when pebbles drift through its orbit. The pebble accretion efficiency (ε_ PA = Ṁ_ PA/ Ṁ_ peb) is taken from  <cit.> and <cit.>. 
This quantity depends on the disk properties (, η, α_ t) and the planet properties (m_ p, a, e, i).  In the limit of 2D and 3D pebble accretion these are given by

    ε_2D = 0.32/η√(m_ pΔ v/M_⋆ v_ Kτ_s) f_set

and

    ε_3D = 0.39 m_ p/η h_peb M_⋆ f_set^2

respectively, where Δ v is the relative velocity between the planet and the pebble, h_peb= H_ peb/r is the aspect ratio of the pebble disk.  In eps-3D we have already assume that i< h_ peb. The above expressions include a modulation factor,  

    f_set = exp[ -0.5 (Δ v/v_∗)^2 ]

where v_∗=(M_p/τ_s M_⋆)^1/3 v_ K [A more general form of f_ set can be found in Eq. (35) of <cit.>, 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_ set≪1 and pebble accretion fails <cit.>.

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 

    Ṁ_ i,peb = {   Ṁ_ peb     i = 1 
       Ṁ_ peb∏_k=1^i-1(1-ε_ k,PA ),         i≥ 2 
    .
 
where i is ordered for bodies from the furthest to closest in terms of semimajor axis and ε_i,PA is the pebble accretion efficiency of the i^th body.
  
We  focus on the formation of protoplanets with m_ p ≳ 1 (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 <cit.>. Therefore, we do not implement  the termination of pebble accretion.





   


 §.§ Planetesimal initial condition



In the streaming instability mechanism, pebbles accumulate into dense filaments. The typical width of the filament  is Δ r ≃η r_ ice <cit.>.  The threshold condition to trigger gravitational collapse requires ρ_ peb≃ρ_ gas.  The total solid mass available to build planetesimals is therefore 2π r_ iceΔ r  Σ_ peb(r_ ice)  =2π r_ iceΔ r  Σ_ gas(r_ ice) H_ peb/H_ gas.  <cit.> 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_⊕. 






§ 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 <cit.>. 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 of pebbles ( 100 of Ṁ_ peb  over 1 Myr), what is the mass growth of planetesimals for the above three different initial size distributions.    





 §.§ 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 (, planet inclination) shaping the planetesimal + pebble accretion  process.

    






  §.§.§ Effect of inclinations on  planetesimals growth


As mentioned in planetesimal initial condition,  the streaming instability  produces in total 0.039  M_⊕ planetesimals at r_ ice. In the mono-dispersed initial conditions,  we assume that all  planetesimals are 400  km in radius. (6.7 × 10^-5 in mass). Therefore, N=580  planetesimals are generated in this compact ring belt (Δ r = η 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_ ice-Δ r /2 to r_ ice + Δ 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_ 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 <cit.>. The planetesimals would be lifted out of the midplane and acquire an inclination distribution.
 

In hydrodynamic simulations <cit.> 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 ≪η).  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).   




2d3d 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 2d3da the mass growth is dominated by planetesimal-planetesmial collision at the beginning (the sudden jump in mass of grey lines in 2d3d). Since the pebble accretion rate is an increasing function of the planet mass,  when the mass approaches  10^-2 M_⊕, the growth is driven by pebble accretion (smooth growth in mass  in 2d3d). Only a few bodies survive after 1 Myr, and the mass of the dominant body is  4  M_⊕. 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 ∼10^-3 to 10^-2 emerges. 
For the realistic case when planetesimals are on inclined orbits  (2d3db),  however, much fewer mergers occur and the mass growth remains modest. The mass of the largest body remains < 10^-2, and there is no efficient pebble accretion at the end of the simulation. 

 



  §.§.§ Excitation of eccentricities and inclinations


eccinc 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 (δ v) of planetesimals is smaller than their Hill velocity (v_ H =R_ H) <cit.>, where the mutual Hill radius  R_ H = (2m_ p/3M_⋆)^1/3 a.  After a few hundred years their eccentricities get excited to values larger than the Hill velocity (δ v > v_ H, dashed line in eccinc), 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_ rms≃ e_ rms/2 <cit.>.


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 <cit.>,
 
    τ_  vs≃   3 × 10^5( i/0.005) [  (e/0.01)^2 + ( i/0.005)^2  ]^3/2(  R_ p/400   km)^-3
       (  Σ_ p/10  g cm^-2)^-1(  ρ_∙/1.5  g cm^-3)^-1(  a/2.7 )^-1/2  yr.

The damping timescale of eccentricity from the gas drag is <cit.>,
 
    τ_  drag≃      10^7(  R_ p/400   km) (  Σ_ gas0/500  g cm^-2)^-1(  h_ gas0/0.033)  
       [ ( e/10^-2)^2 +( i/5 × 10^-3)^2 +  ( η/2.5 × 10^-3)^2 ]^-1/2
       (  ρ_∙/1.5  g cm^-3) (  a/2.7 )^11/4  yr.

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 <cit.> 
 
    τ_ tidal≃( M_⋆/ m_ p)  ( M_⋆/Σ_ gas a_ p^2  )    (  h_ gas^4 /Ω_ K)
       ≃ 7 × 10^7(  R_ p/400   km)^-3
       (  Σ_ gas0/500  g cm^-2)^-1(  h_ gas0/0.033)^4    (  ρ_∙/1.5  g cm^-3)^-1(  a/2.7 )^3/2  yr.

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). 
From the above equations we obtain τ_ vs≪τ_ drag  <  τ_ 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 2d3db? 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 (eps-2Deps-3D), less massive bodies accrete pebbles more slowly.  Second, pebble accretion is quenched due to a high eccentricity and inclination. Specifically, from fset0 it follows that pebble accretion requires encounters to have a sufficiently low relative velocity

    Δ v ≲ v_∗ = ( M_p/M_⋆/τ_ s)^1/3 v_ K.

Inserting Δ v ≃ ev_ K and e ∼ 10^-2 we find that pebble accretion becomes only significant when the mass of the planet approaches e^3 τ_ s M_⋆≈ 10^-2 M_⊕.  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_ rms> h_ peb) <cit.>. In eccinc we find that this happens at t≃ 5 × 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”.




 
   
  
 


  §.§.§ Spreading of the planetesimal belt


In expansion we plot the time evolution of the planetesimals' semimajor axes  in  run_md_1. We define the full width of the planetesimal belt as Δ W = 2√(∑_i (a_i -r_ ice)^2/N). The initial width is η  r_ 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 expansion 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. <cit.> obtain an analytical expression of the viscosity due to the mutual gravitational scattering based on equal size planetesimals. Substituting their viscosity (ν in their Eq. (18)) into Δ W = √(ν t),  we have   
 
    Δ W =  C_ w√(( N   Ω_ K t /e^2  ) ) R_ H^2/a,
 
where  the prefactor C_ w^3 = C_ fit  [24  I(β)/π^2ln(Λ ^2 + 1)], I(β)=0.2,  t is the time, Λ=  i̅_  rms( e̅_ rms^2 +  i̅_ rms^2)/3, and  i̅_ rms = a i_ rms/R_ H and e̅_ rms = a e_ rms/R_ H. The numerical factor of C_ fit = 30 is calibrated with our numerical simulations and is ≃ 3 times larger than <cit.>'s result. This formula is obtained under the assumption that planetesimals are in the dispersion-dominated accretion regime. 
We show in expansion that our analytical expression (width) agrees well with the simulation (black line).
For fixed total mass Nm_ p and a, Δ W increases with time (∝ t^1/3), the initial size of the planetesimal (∝ 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 <cit.>.
 
  
 It is clearly seen from 2d3db that the growth of planetesimals is slower than 2d3da 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 of pebbles drift through. 



 §.§ Poly-dispersed initial conditions with large planetesimals









In the previous section we found that when the sizes of planetesimals are all 400   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. <cit.> find that the initial size can be described by a single power law.  Recently <cit.> 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 <cit.>).
From  <cit.> simulations, the number fraction of planetesimals with m_ p>m is given by their Eq.(19),
 
    N_>(m)/N_ tot = ( m/ m_ min)^-pexp[  ( m_ min/m_ e)^β  - ( m/ m_ e)^β].
 
Following this distribution, we adopt their  parameters: p = 0.6, β = 0.35.  The minimum mass and  exponential characteristic mass are fitting parameters, which are given as  m_ min =10^-6 (100 km), m_ e=1.6×10^-5 (250 km), respectively.  SImass 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 planetesimals. Here we simulate the most massive 470 bodies from this population, corresponding to the two-thirds of the total mass (blue brach in SImass). 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 poly-disperseda while the semi-major axis of the largest body is shown in  poly-dispersedb.  In poly-dispersedc 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 (SImass).  Energy equipartition leads to a configuration that the massive planetesimals have low random velocities while the low-mass planetesimals have high random velocities  <cit.>.  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 poly-disperseda) can accrete more neighbouring planetesimals. Since its eccentricity and inclination remains low,  once the mass  is beyond  a  few 10^-3, 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.  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   ∼ 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_ p and low Δ v) facilitate the transition to the pebble accretion dominated phase (phase B)  and promote the formation of a  super-Earth planet after ≃1 Myr.

   



 §.§ 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 <cit.>. For the single-site planetesimal formation scenario, the pebble density is locally enhanced at a particular disk location (, the ice line). But how strong the enhancement  is determined by complicated  physical processes such as water vapor diffusion and condensation <cit.>. Taken these factors and uncertainties into account, we also consider the case of an initially small planetesimal size.  Nevertheless, we argue in runaway that in this circumstance runaway growth of planetesimals proceeds fast, producing the desired two-component mass distribution that promotes the formation of massive planets. 



  §.§.§ Runaway growth in a narrow planetesimal belt



The characteristic runaway planetesimal accretion timescale is <cit.>

    τ_ rg    =  C_ rgρ_∙ R_0/Σ_ p 
      = 10^5(  C_ rg/0.1)  (  R_0/100  km) 
       (   Σ_ p/10  g cm^-2)^-1( a/ 2.7 )^3/2   yr.
 
where C_ rg≃0.1 is a numerical-corrected prefactor from <cit.>, and ρ_∙ =1.5   g  cm^-3 is the internal density of the planetesimal. 
Small planetesimals have three advantages in terms of  mass growth.
First, based on runaway 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 (tau_vs), which also boosts the planetesimal and pebble accretion. 
Third, as shown by in spread, the spreading of their semimajor axes  increases with the size of the planetesimal (R_0 ∝ R_ H in width). 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    km  in size (≃ 10^-6 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 (d N/dm_ p∝ m_ p^-2.5  <cit.>, most of the mass remains in small planetesimals. As discussed in poly-dispersed, 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 ∼ 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 poly-disperseda).
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 τ_ vs  (tau_vs) equals  the runaway timescale τ_ rg (runaway) [This criterion is different from <cit.>. See discussion in <cit.>. ]. When τ_ vs > τ_ rg, the stirring of the random velocities is smaller compared to the accretion, which is the key feature of runaway growth. When τ_ vs < τ_ 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 <cit.>,
 
    R_ rg/oli≃     850 (C_ rg/ 0.1)^3/7(R_0/ 100   km)^3/7(a/ 2.7  )^5/7(Σ_ p/  1   g  cm^-2)^2/7 km.
 
Substituting runaway into width, we obtain that when the massive body grows into a transition radius, the planetesimal belt width is   

    Δ W  =  √(C_ rg C_ w^3 R_ H^3/2R_0).
 
 We note that R_ H is the Hill radius of the planetesimal with size of R_0. 
From  Σ_ p =  M_ tot,pl/2π  a Δ W  and  and  width2,  we calculate that the belt width extends to Δ W =0.029 AU and the planetesimal surface density reduces to Σ_ p = 2  g cm^-2 after the runaway growth (∼ 10^5 yr). 
Substituting the above valuses into Rtr, we obtain, 

    R_ rg/oli≃   1000 (C_ rg/ 0.1)^3/7(R_0/ 100   km)^1/7(M_ tot,pl/ 0.039 )^2/7  km.
 
We therefore conclude that the feeding zone of this embryo (≃ 10  R_ H, <cit.>) is similar to the width of the planetesimal belt, supporting one single embryo assumption.








  §.§.§ The superparticle approach


Since the total planetesimal mass is 0.039  and small planetesimals dominate the mass, there are  in total  35000 100-km-sized planetesimals (≃10^-6) 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 <cit.> and references therein), in which N_ 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_ sp times higher (m_ sp= N_ 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 ≪  m_ sp≪ m_ p, where m_ 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_ sp,  while  the dynamical evolution  may become artificially stochastic for a too large N_ 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 (δ v ≃ v_ esc = √(2G m_0/R_0)),  , e ≃ e_ esc= v_ esc/v_ K. The eccentricities and inclinations are adopted to follow the Rayleigh distributions, where  e_0 =2 i_0= e_ esc. Their semi-major axes are randomly initialized  from r_ ice(1-η/2) to r_ ice(1+ η/2).       In addition, since pebble accretion is extremely inefficient for 100 km size planetesimals (the inequality (<ref>) is far from satisfied), it is  reasonable to neglect pebble accretion of the superparticles.  The embryo is initially placed at r_ ice.  



  §.§.§ Results



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

 

The disk and pebble parameters in this (run_sp) are taken to be identical to  run_pd in poly-dispersed. The mass of the embryo is chosen to be the transition value  (m_ p = 10^-3 when R_ p = R_ rg/oli  in  Rtr2).
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. 


runaway 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 (, m̅ (t) =∑_i=1^3 m_ 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, the type I migration becomes important.  The planet would take  ∼ 0.5 Myr to migrate into the inner region of the disk. 
Eventually, a 1.3  planet forms within 1 Myr. 




  §.§.§ Timescale analysis


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

    τ_ PA,2D =    m_ p/ε_ PA,2DṀ_ peb≃ 9 × 10^4 f_set^-1(  m_ p/0.05 )^1/3
       ( /0.1)^1/3( η/2.5 × 10^-3)  ( Ṁ_ peb/  100  M_⊕/Myr)^-1 yr.

where we use eps-2D and Δ v is assumed to be dominated by the Keplerian shear since the eccentricity of the embryo is  insignificant due to dynamical friction. From fset0 the above f_set evaluates  as

    f_set   = exp[-0.07   (η/2.5× 10^-3)^2   (m_ p/0.01 )^-2/3( /0.1) ^2/3].

When the planet is ≳ 10^-2, f_ set≃ 1.  
From fset, we clearly see that  f_ set becomes smaller when a planet is less massive  or their pebble-planet relative velocity is higher. 

Similarly, in the 3D limit (eps-3D), the pebble accretion timescale becomes
 
    τ_ PA,3D =    m_ p/ε_ PA,3DṀ_ peb≃ 9 × 10^4 f_set^-2(  h_ peb/ 4.2 × 10^-3)  
       ( η/2.5 × 10^-3)  ( Ṁ_ peb/  100  M_⊕/Myr)^-1 yr.

From the planet mass dependence on  tau_PA2dtau_PA3d, we know that the pebble accretion is in 3D  when the planet  mass is low, and it transitions to 2D accretion when the planet becomes more massive. The transition mass between 2D and 3D is ≃ 0.05 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_ p =3× 10^-4) 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 analysis, where the grey,  blue,  light blue and  magenta  lines represent growth timescales obtained from the simulation,  the pebble accretion in 3D (with and without f_ set) and in 2D, respectively. We can identify the following different growth stages.



 



When the mass of the embryo is lower than ∼ 10^-3, 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 analysis.  The pebble accretion is inefficient  in this phase because of the low embryo's mass (small m_ p in tau_PA3d) and hence the weak setting accretion (small f_ set^2  in tau_PA3d). When the mass of the embryo becomes larger, 3D 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. With increasing mass f_ set increases. When f_set approaches unity, the growth timescale becomes independent of the planet mass (see tau_PA3d and the light blue in analysis). Finally,  growth enters  the 2D  pebble accretion regime when the mass approaches 0.1. In that case growth timescale starts to increase,  τ_ PA,2D∝ m_ p^1/3 (magenta line in analysis), since ε_ PA,2D no longer scales linearly with the planet mass.






    



 §.§ 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_ rg/oli∝ R_0^1/7, Rtr2), 
and 2) since the orbital expansion is less significant for smaller planetesimals (Δ W ∝ R_0 from width2) and Σ_ p∝Δ W^-1,
 the runaway accretion timescale is strongly correlated with the initial size ( t_ rg∝ R_0^2 from runaway). 
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 runaway 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  <cit.>. This is because f_ set becomes smaller when  η increases with r.  Therefore, both two effects further limit the growth of the protoplanets from a large size planetesimal (, 400 km) when the planetesimal formation occurs at a large distance.


 
Our findings are similar to <cit.>. 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 (, <cit.> and <cit.>).  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. 







§ 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.  As illustrated in runaway,  the pebble accretion dominates the growth and the planetesimals' contribution can be neglected after m_ p≳ 10^-3.  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 ȧ/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 parameter  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 and  6.8, 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  in run_fid while  in run_nmig the planet always stays at the ice line and attains 3.8. 
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_ p≳ 1 does not form within 1 Myr from these two configurations.  In  run_tau the planet attains 4.8 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 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  and migrates into the inner disk cavity at t = 0.7 Myr in run_hpeb, the planet mass is only 0.02  in run_lpeb.  We find  a super linear correlation between the planet mass (m_ p)  and  the integrated pebble flux (Ṁ_peb t).
 A factor of 4 change in  Ṁ_ pebt results in more than two order of magnitude growth in M_ p. 
 This is due to the fact that pebble accretion efficiency (ε_ PA) increases with m_ p   in both 2D and 3D 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 (Hpeb).  The transition mass  from 3D  to 2D pebble accretion  will also become lower. Therefore, in run_alpha the embryo starts efficient 2D 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 α_ 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 tau_PA2dtau_PA3d transition mass  from 3D to 2D in run_tau is 0.15,  three times higher than the fiducial run. In this case starting with 10^-3, the planet grows a significant fraction of its mass in the slow, 3D pebble accretion regime. On the other hand, when the planet enters 2D pebble accretion, the accretion is faster when the Stokes number is lower (tau_PA2d). 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.


 


    


§ 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  (, 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, , typically following the surface density distribution of the gas.  In simulations their planetesimal surface density remains constant by suppling material from the neighbouring region <cit.>. 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  ≃ 10^-3. 

We have modified the  N-body code to perform  simulations of the mass growth and orbital evolution of these planetesimals (method).  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 <cit.> and <cit.> that accounts for the disk  parameters (, η, α and Ṁ_ peb) and planet properties (m_ 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 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 (mono-dispersed).   
    
  * 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, 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 (poly-dispersed).   
    
  * Alternatively, protoplanets also form out of their birth ring when the initial size of the planetesimals is small (, 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 (two-component) 
    
  * 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 (parameter).       


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

To address these issues, a global disk model is needed. Using a novel Lagrangian approach, such a model has just been developed <cit.>. 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 <cit.>.

 



§ CONVERGENCE TEST FOR SUPERPARTICLE SIMULATIONS



We show the simulations with different mass of the superparticle (m_ sp = N_ sp m_ 0) for run_sp starting with m_ p = 3 × 10^-4 in sp_timescale. Three simulations with randomlized initial conditions are performed for each set of N_ sp.  The results are shown in compare where red, green represent  N_ sp = 25 and 50, respectively. The thick line represents the mean value averaged from three individual simulations, for instance, m̅ (t) =∑_i=1^3 m_ i(t)/3 whereas the light region marks the range between the minimum and maximum values from the three simulations.
 
compare 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_ sp.   In comparea the difference between the mean mass  for the above N_ sps 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 comparec. 
  
 
 To summarize, results from the above three tested N_ sp are in general agreement with each other. 

 








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.







aa



