Context. Streaming instability is a key mechanism in planet formation clustering pebbles into planetesimals with the help of selfgravity. 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. BR  Aims We aim to investigate the growth of planetesimals into protoplanets at a single location through streaming instability. For a solarmass star we test the conditions under which superEarths are able to form within the lifetime of the gaseous disk. BR  Methods We modified the Mercury Nbody 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 were investigated i a monodispersed population of 400 km radius planetesimals ii a polydispersed population 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. BR  Results The monodispersed population of 400 km size planetesimals cannot form protoplanets of a mass greater than that of the Earth. Their eccentricities and inclinations are quickly excited which suppresses both planetesimal accretion and pebble accretion. Planets can form from the polydispersed and bimodal distributions. In these circumstances it is the twocomponent nature that damps the random velocity of the large embryo through the dynamical friction of small planetesimals allowing the embryo to accrete pebbles efficiently when it approaches 10SUP2SUP MSUBSUB. Accounting for migration closein superEarth planets form. SuperEarth planets are likely to form when the pebble mass flux is higher the disk turbulence is lower or the Stokes number of the pebbles is higher. BR  Conclusions For the single site planetesimal formation scenario a twocomponent 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 superEarth planets by the combined mechanisms of planetesimal accretion and pebble accretion.
['methods: numerical', 'planets and satellites: formation', 'Astrophysics - Earth and Planetary Astrophysics']
Growth after the streaming instability. From planetesimal accretion to pebble accretion
{'From planetesimal accretion to pebble accretion': 'Beibei Liu 1 , 2 , Chris W. Ormel 2 , Anders Johansen 1  \n- 1 Lund Observatory, Department of Astronomy and Theoretical Physics, Lund University, Box 43, 22100 Lund, Sweden.\n- 2 Anton Pannekoek Institute (API), University of Amsterdam, Science Park 904,1090GE Amsterdam, The Netherlands\n- e-mail: bbliu@astro.lu.se,c.w.ormel@uva.nl,anders@astro.lu.se  \nFebruary 27, 2019', 'ABSTRACT': "Context. Streaming instability is a key mechanism in planet formation, clustering pebbles into planetesimals with the help of selfgravity. 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.  \nAims. 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.  \nMethods. 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.  \nResults. The mono-disperse population of 400 km size planetesimals cannot form glyph[greaterorsimilar] 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 M ⊕ . 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.  \nConclusions. 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.  \nKey words. methods: numerical planets and satellites: formation", '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'erez et al. 2015; Tazzari et al. 2016). But further growth is suppressed by bouncing or fragmentation due to the increasing compactification in collisions (Guttler 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).  \nThe 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  \napproximate a few hundred kilometres (Johansen et al. 2012, 2015; Simon et al. 2016; Schafer et al. 2017; Simon et al. 2017; Abod et al. 2018).  \nThe 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 ( d m p / d t ∝ m γ p with γ > 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 d m p / d t ∝ m γ p with γ < 1 (Kokubo & Ida 1998). This phase is characterized  \nby a decreasing mass ratio of two adjacent massive runaway bodies (Lissauer 1987; Kokubo & Ida 2000; Thommes et al. 2003; Ormel et al. 2010).  \nIn 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 d m p / d t ∝ m 2 / 3 p (Lambrechts & Johansen 2014, Hill regime). Onthe 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 d m p / d t ∝ m p (Ida et al. 2016).  \nIn 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 (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 ε PA under general circumstances (Liu & Ormel 2018; Ormel & Liu 2018). 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.  \nProgress 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 H 2 O 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.  \nConstrained by the size distribution of planetesimals in the asteroid belt, Morbidelli et al. (2009) concluded that their born size should be large ( glyph[greaterorsimilar] 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 glyph[greaterorsimilar] 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.  \nIn 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.  \nThe 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 polydispersed 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¸˙zkowska &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 ice = 2 . 7 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.  \nThis 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 ( ∆ r = η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  \nΣ gas = Σ gas0 ( r 1 AU ) -1 , (1)  \nh gas = h gas0 ( r 1 AU ) 1 / 4 , (2)  \nFig. 1. Sketch illustrating the initial setup. Planetesimals (purple circles) have formed at a single location by streaming instability. They accrete both among themselves as well as from the inward drifting pebbles (blue). The pebble disk and gas disk are marked in light blue and light green, where the gas and pebble scale heights are H gas , H peb , respectively. The horizontal black arrow represents the radial width ( ηr ) of the forming planetesimal ring. Due to dynamical friction, large planetesimals have low eccentricities and inclinations whereas small planetesimals have high velocity dispersions.  \n<!-- image -->  \nwhere Σ 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 (Garaud & Lin 2007; Bitsch et al. 2015). Here Σ gas0 = 500 g cm -2 and h gas0 = 0 . 033 (Hayashi 1981) are adopted as the default values in this paper. Therefore, the disk temperature T and η are  \nT = µh 2 gas GM glyph[star] rR g = T 0 ( r 1 AU ) -1 / 2 , (3)  \nand  \nη = -h 2 gas 2 ∂ log P gas ∂ log r = η 0 ( r 1 AU ) 1 / 2 , (4)  \nwhere G is the gravitational constant, R g = 8 . 31 × 10 7 erg mol -1 K -1 is the gas constant, M glyph[star] 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.  \nAlthough 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 (Andrews et al. 2009). The fiducial value of Σ 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 ∼ 170K ) is located at r 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.  \nFor 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 appropri-  \nate. 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.  \nThe 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),  \na drag = -( 3 C D ρ gas 8 R p ρ · ) v rel v rel , (5)  \nwhere 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, v K 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.  \nLarge 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)  \nt m glyph[similarequal] 0 . 5 ( M glyph[star] m p )( M glyph[star] Σ gas a 2 p ) ( h 2 gas Ω K ) glyph[similarequal] 5 × 10 5 ( m p 1 M ⊕ ) -1 ( Σ gas0 500 gcm -1 ) -1 ( h gas0 3 . 3 × 10 -2 ) 2 ( a p 2 . 7 AU ) yr , (6)  \nwhere 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.  \nThe accelerations acting on planets due to type I migration torque are  \na m = -v t m , a e = -2 ( v · r ) r r 2 t e , a i = -v z t i , (7)  \nwhere 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 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.  \n2015; 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 a 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 r = -2 ηv K τ s / (1 + τ 2 s ) (Weidenschilling 1977) where τ s is the dimensionless stopping time (Stokes number). The drift speed is determined by the aerodynamical size of the pebble τ s and η . Based on Lambrechts & Johansen (2014) and Schoonenberg et al. (2018), the pebble mass flux ( ˙ M peb = 2 πrv 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 ˙ M peb = 100 M ⊕ / Myr (consistent with Eq. (14) of Lambrechts & Johansen (2014)) and neglect the time dependence for simplicity. A lower ˙ M 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 (Youdin & Lithwick 2007)  \nH peb = √ α t α t + τ s H gas . (8)  \nThe pebble flux therefore needs to be smaller than 2 πrv r Σ gas H peb /H gas . Our default value is well below this limit. We also explore two additional pebble fluxes ˙ M peb = 200 M ⊕ Myr -1 and 50 M ⊕ Myr -1 in Sect. 4.  \nThe 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 M glyph[circledot] yr -1 , the dust mass is ∼ 100 M ⊕ (Fig.1 of Manara et al. (2016)), consistent with our adopted fiducial value.  \nIn Eq. (8) α 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 (Johansen & Klahr 2005; Zhu et al. 2015; Yang et al. 2018). The value of α 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 α t = 10 -3 and test the influence of a lower turbulent disk ( α t = 10 -4 ) in Sect. 4.  \nFor pebbels we adopt a fiducial aerodynamical size of τ 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 τ s = 0 . 03 is also explored in Sect. 4.  \nA fraction of pebbles will be accreted onto a planetesimal when pebbles drift through its orbit. The pebble accretion effi-  \nciency ( ε PA = ˙ M PA / ˙ M peb ) is taken from Liu & Ormel (2018) and Ormel & Liu (2018). This quantity depends on the disk properties ( τ s , η , α t ) and the planet properties ( m p , a, e, i ). In the limit of 2D and 3D pebble accretion these are given by  \nε 2D = 0 . 32 η √ m p ∆ v M glyph[star] v K τ s f set (9)  \nand  \nε 3D = 0 . 39 m p ηh peb M glyph[star] f 2 set (10)  \nrespectively, 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 Eq. (10) we have already assume that i < h peb . The above expressions include a modulation factor,  \nf set = exp [ -0 . 5(∆ v/v ∗ ) 2 ] (11)  \nwhere v ∗ = ( M p /τ s M glyph[star] ) 1 / 3 v K 1 . Physically, when the pebbleplanet 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 glyph[lessmuch] 1 and pebble accretion fails (Visser & Ormel 2016).  \nIn 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  \n˙ M i , peb =        ˙ M peb i = 1 ˙ M peb i -1 ∏ k =1 (1 -ε k , PA ) , i ≥ 2 (12)  \nwhere 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.  \nWe focus on the formation of protoplanets with m p glyph[greaterorsimilar] 1 M ⊕ (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 ∆ r glyph[similarequal] ηr ice (Yang & Johansen 2016; Li et al. 2018). The threshold condition to trigger gravitational collapse requires ρ peb glyph[similarequal] ρ 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 . 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 ⊕ .  \nFigure 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 ⊕ , 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 ⊕ . Even though the accretion timescale in this  \n<!-- image -->  \nFig. 2. Mass growth of planetesimals for planets on coplanar orbits (left, run md test in Table 1) and non-coplanar orbits (right, run md 1 ). Grey lines trace the mass of each planetesimal. Red dots represents the merging of two bodies, marked at the position of the less massive body. The early sudden growth is caused by planetesimal accretion and the later smooth growth by pebble accretion.  \n<!-- image -->', '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; Schafer 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.  \nIn this section we investigate the question that provided in total 100 M ⊕ of pebbles ( 100 M glyph[circledot] yr -1 of ˙ M 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 ⊕ planetesimals at r ice . In the mono-dispersed initial conditions, we assume that all planetesimals are 400 km in radius. ( 6 . 7 × 10 -5 M ⊕ 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 / 2 e 2 0 ] .  \nWe 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 (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.  \nIn hydrodynamic simulations (Simon et al. 2016; Schafer 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 glyph[lessmuch] η ). 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).  \nTable 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 Myr.  \nconfiguration 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 accretiondominated phase (phase B) when a massive body of ∼ 10 -3 M ⊕ to 10 -2 M ⊕ 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 M ⊕ , 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 ( δv ) of planetesimals is smaller than their Hill velocity ( v H = R H Ω K ) (Ida 1990), where the mutual Hill radius R H = (2 m p / 3 M glyph[star] ) 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 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 rms glyph[similarequal] e rms / 2 (Ida et al. 1993).  \nBased 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),  \nThis 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)  \nτ tidal glyph[similarequal] ( M glyph[star] m p )( M glyph[star] Σ gas a 2 p ) ( h 4 gas Ω K ) glyph[similarequal] 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 AU ) 3 / 2 yr . (15)  \nThe 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 M ⊕ ). From the above equations we obtain τ vs glyph[lessmuch] τ 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.  \nThe 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  \nτ vs glyph[similarequal] 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 AU ) -1 / 2 yr .  \n(13)  \nThe damping timescale of eccentricity from the gas drag is (Adachi et al. 1976),  \nτ drag glyph[similarequal] 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 AU ) 11 / 4 yr . (14)  \n∆ v glyph[lessorsimilar] v ∗ = ( M p /M glyph[star] τ s ) 1 / 3 v K . (16)  \nInserting ∆ v glyph[similarequal] 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 glyph[star] ≈ 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 ) (Levison et al. 2015). In Fig. 3 we find that this happens at t glyph[similarequal] 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'.", '3.1.3. Spreading of the planetesimal belt': "In Fig. 4 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  \nFig. 3. Root-mean-square (rms) eccentricity (orange) and inclination (cyan) of planetesimals as functions of time in run md 1 . The dashed line indicates the separation between the sheardominated and dispersion-dominated regime, e = v H /v K . The dotted line indicates the inclination equal to the aspect ratio of the pebble disk.  \n<!-- image -->  \nFig. 4. Orbital expansion of planetesimals as a function of time. The grey and black lines are the semimajor axis of each planetesimal and the mean value of the planetesimal population, respectively. The red line denotes the analytical expression for the width of the planetesimal belt from Equation (17).  \n<!-- image -->  \n⋂]⌉ ˜ {glyph[arrowbt]√}  \nis η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 Fig. 4 we find that the width increases steadily with time and becomes 0 . 05 AU after 10 5 yr.  \nThis 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 ( ν in their Eq. (18)) into ∆ W = √ νt , we have  \n∆ W = C w 3 √ ( N Ω K t/e 2 ) R 2 H /a, (17)  \nwhere the prefactor C 3 w = C fit [24 I ( β ) /π 2 ln(Λ 2 +1)] , I ( β ) = 0 . 2 , t is the time, Λ = ¯ i rms (¯ e 2 rms + ¯ i 2 rms ) / 3 , and ¯ i rms = ai rms /R H and ¯ e rms = ae rms /R H . The numerical factor of C fit = 30 is calibrated with our numerical simulations and is glyph[similarequal] 3 times larger than Ohtsuki & Tanaka (2003)'s result. This  \nFig. 5. Sum of planetesimal mass in each bin for all mass planetesimals (red + blue) and for the simulated mass branch (blue) vs the initial planetesimal mass (size) from Eq. (18). The initial mass function of planetesimals was based on streaming instability simulations with a minimum mass of 10 -6 M ⊕ ( 100 km) and a characteristic mass of 1 . 6 × 10 -5 M ⊕ ( 250 km). We discrete the this distribution into different mass bins.  \n<!-- image -->  \nformula 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 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 ).  \nSince 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).  \nIt 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.  \nIn 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 M ⊕ 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 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.  \nFig. 6. Panel (a): mass growth of the planetesimals for a polydispersed size population of planetesimals based on the blue branch of Fig. 5. The black line represent the growth of the final most massive body and grey lines are for other planetesimals. Red dots represents the merging of two bodies, marked at the position of the less massive body. The sudden growth is caused by planetesimal accretion and the smooth growth by pebble accretion. Note that f set approaches 1 when the planet is 0 . 01 M ⊕ . Panel (b): semi-major axis evolution of the final most massive body. Panel (c): Eccentricity (solid) and Inclination (dashed) evolution where black lines represent the final most massive body and the red lines indicate rms values of the planetesimals. The simulation is from run pd 1 .  \n<!-- image -->  \n6  \nThe 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 Schafer 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 Schafer et al. (2017)). From Schafer et al. (2017) simulations, the number fraction of planetesimals with m p > m is  \ngiven by their Eq.(19),  \nN > ( m ) N tot = ( m m min ) -p exp [ ( m min m e ) β -( m m e ) β ] . (18)  \nFollowing 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 M ⊕ ( 100 km), m e = 1 . 6 × 10 -5 M ⊕ ( 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 M ⊕ 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.  \nOther 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 semimajor 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.  \nFor 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 M ⊕ , 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.  \nRather 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 M ⊕ . 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.  \nTo conclude, the presence of ∼ 10 3 kmplanetesimals 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 glyph[similarequal] 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  \nwould 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 singlesite 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)  \nτ rg = C rg ρ · R 0 Σ p Ω K = 10 5 ( C rg 0 . 1 )( R 0 100 km ) ( Σ p 10 g cm -2 ) -1 ( a 2 . 7 AU ) 3 / 2 yr . (19)  \nwhere C rg glyph[similarequal] 0 . 1 is a numerical-corrected prefactor from Ormel et al. (2010), 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 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 ∝ R 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.  \nAccounting for these effects, we assume that the planetesimals are all 100 km in size ( glyph[similarequal] 10 -6 M ⊕ 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 p ∝ m -2 . 5 p (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.  \nThe 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  \nin Fig. 6a). Altogether, the situation resembles the classical runaway/oligarchy transition, except that, in our case, there is only one single 'oligarch'.  \nWhy 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 (Eq. (13)) equals the runaway timescale τ rg (Eq. (19)) 2 . 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 Ormel et al. (2010),  \nR rg / oli glyph[similarequal] 850 ( C rg 0 . 1 ) 3 / 7 ( R 0 100 km ) 3 / 7 ( a 2 . 7 AU ) 5 / 7 ( Σ p 1 g cm -2 ) 2 / (20)  \nSubstituting Eq. (19) into Eq. (17), we obtain that when the massive body grows into a transition radius, the planetesimal belt width is  \n∆ W = √ C rg C 3 w R 3 H 2 R 0 . (21)  \nWe 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 Eq. (21), 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 Eq. (20), we obtain,  \nR rg / oli glyph[similarequal] 1000 ( C rg 0 . 1 ) 3 / 7 ( R 0 100 km ) 1 / 7 ( M tot , pl 0 . 039 M ⊕ ) 2 / 7 km . (22)  \nWe therefore conclude that the feeding zone of this embryo ( glyph[similarequal] 10 R 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 M ⊕ and small planetesimals dominate the mass, there are in total 35000 100-kmsized planetesimals ( glyph[similarequal] 10 -6 M ⊕ ) 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 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 glyph[lessmuch] m sp glyph[lessmuch] 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  \nartificially stochastic for a too large N sp (low number of superparticles).  \nWe 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 glyph[similarequal] v esc = √ 2 Gm 0 /R 0 ), i.e., e glyph[similarequal] 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 (16) is far from satisfied), it is reasonable to neglect pebble accretion of the superparticles. The embryo is initially placed at r ice .", '3.3.3. 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 Sects. 3.3.3 and 3.3.4 are based on the N sp = 50 run.  \nFig. 7. Panel (a): Mass growth of the embryo for the twocomponent initial conditions in run sp . Panel (b): Semi-major axis evolution of the embryo. The red thick line is for the averaged value and the light red is for the spreading among three different runs. The black dashed line is for the growth of a single embryo on a circular and coplanar orbit without any planetesimal accretion.  \n<!-- image -->  \nThe 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 p = 10 -3 M ⊕ when R p = R rg / oli in Eq. (22)). Three individual simulations ( run sp 1 to run sp 3 ) are conducted, varying only the initial  \npositions 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.  \nFigure 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., ¯ m ( t ) = ∑ 3 i =1 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.  \nWhen the mass of the embryo is beyond 0 . 1 M ⊕ , 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 M ⊕ 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  \nτ PA , 2D = m p ε PA , 2D ˙ M peb glyph[similarequal] 9 × 10 4 f -1 set ( m p 0 . 05 M ⊕ ) 1 / 3 ( τ s 0 . 1 ) 1 / 3 ( η 2 . 5 × 10 -3 ) ( ˙ M peb 100 M ⊕ / Myr ) -1 yr . (23)  \nwhere we use Eq. (9) and ∆ 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 set evaluates as  \nf set = exp [ -0 . 07 ( η 2 . 5 × 10 -3 ) 2 ( m p 0 . 01 M ⊕ ) -2 / 3 ( τ s 0 . 1 ) 2 / 3 ] . (24)  \nWhen the planet is glyph[greaterorsimilar] 10 -2 M ⊕ , f set glyph[similarequal] 1 . From Eq. (24), we clearly see that f set becomes smaller when a planet is less massive or their pebble-planet relative velocity is higher.  \nSimilarly, in the 3D limit (Eq. (10)), the pebble accretion timescale becomes  \nτ PA , 3D = m p ε PA , 3D ˙ M peb glyph[similarequal] 9 × 10 4 f -2 set ( h peb 4 . 2 × 10 -3 ) ( η 2 . 5 × 10 -3 ) ( ˙ M peb 100 M ⊕ / Myr ) -1 yr . (25)  \nFrom 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 Dand 3 Dis glyph[similarequal] 0 . 05 M ⊕ for the adopted parameters.  \nIn 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  \nFig. 8. Growth timescale as a function of the embryo's mass. The grey line represents the simulation whereas the blue and magenta lines are timescales for the analytical 3 Dand 2 Dpebble accretion. The dark (light) blue line corresponds to τ 3D including (or excluding) the f set factor (Eq. (25)).  \n<!-- image -->  \nsimulations starting with a less massive embryo ( m p = 3 × 10 -4 M ⊕ ) 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.  \nThese 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 set ) and in 2 D, respectively. We can identify the following different growth stages.  \nWhen the mass of the embryo is lower than ∼ 10 -3 M ⊕ , 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 p in Eq. (25)) and hence the weak setting accretion (small f 2 set 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 M ⊕ . With increasing mass f set increases. When f 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 M ⊕ . In that case growth timescale starts to increase, τ PA , 2D ∝ m 1 / 3 p (magenta line in Fig. 8), since ε 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 rg / oli ∝ R 1 / 7 0 , Eq. (22)), and 2) since  \nthe orbital expansion is less significant for smaller planetesimals ( ∆ W ∝ R 0 from Equation (21)) and Σ p ∝ ∆ W -1 , the runaway accretion timescale is strongly correlated with the initial size ( t rg ∝ R 2 0 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.  \nThe 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 set becomes smaller when η 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.  \nOur 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 M ⊕ . As illustrated in Fig. 7, the pebble accretion dominates the growth and the planetesimals' contribution can be neglected after m p glyph[greaterorsimilar] 10 -3 M ⊕ . 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.  \nIn 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/a = 0 . In this  \nTable 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 AU.  \ncase, 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.  \nWe 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 M ⊕ and 6 . 8 M ⊕ , 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 M ⊕ in run fid while in run nmig the planet always stays at the ice line and attains 3 . 8 M ⊕ . 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 glyph[greaterorsimilar] 1 M ⊕ does not form within 1 Myr from these two configurations. In run tau the planet attains 4 . 8 M ⊕ 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 M ⊕ planet at t = 2 . 1 Myr.  \nThe 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 M ⊕ and migrates into the inner disk cavity at t = 0 . 7 Myr in run hpeb , the planet mass is only 0 . 02 M ⊕ in run lpeb . We find a super linear correlation between the planet mass ( m p ) and the integrated pebble flux ( ˙ M peb t ). A factor of 4 change in ˙ M peb t 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 2 Dand 3 Dregime. 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.  \nWhen 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 superEarth planet in a lower turbulent disk.  \nThe 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 Eqs. (23) and (25) transition mass from 3 D to 2 D in run tau is 0 . 15 M ⊕ , three times higher than the fiducial run. In this case starting with  \n10 -3 M ⊕ , 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.  \nThe 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.  \nFig. 9. Mass growth of one single embryo in Sect. 3.3 for different runs given in Table 2. The red, cyan, orange and magenta corresponds to run fid , run nmig , run tau and run alpha , respectively. The thick and thin red lines represent run hpeb and run lpeb .  \n<!-- image -->", '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.  \nThis 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  \nof 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 glyph[similarequal] 10 -3 M ⊕ .  \nWe 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 ( τ s , η , α and ˙ M peb ) and planet properties ( m p , a, e, i ). Simulations with different initial planetesimal sizes and disk parameters are investigated.  \nThe key findings of this study are the following:  \n- 1. Protoplanets cannot emerge from a mono-dispersed population of 400 km size planetesimals, fuelled by a 100 M ⊕ 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).\n- 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 M ⊕ , 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).\n- 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)\n- 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).  \nPlanetesimal 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.  \nIn 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  \nthen emerge sequentially and a chain of multiple planets forms, as envisioned by Ormel et al. (2017).  \nTo 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': 'Weshow the simulations with different mass of the superparticle ( m sp = N sp m 0 ) for run sp starting with m p = 3 × 10 -4 M ⊕ in Sect. 3.3.4. Three simulations with randomlized initial conditions are performed for each set of N sp . The results are shown in Fig. A.1 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 ) = ∑ 3 i =1 m i ( t ) / 3 whereas the light region marks the range between the minimum and maximum values from the three simulations.  \nFig. A.1 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 Fig. A.1a the difference between the mean mass for the above N 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. A.1c.  \nTo summarize, results from the above three tested N sp are in general agreement with each other.  \nAcknowledgements. 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.', 'References': "Abod, C. P., Simon, J. B., Li, R., et al. 2018, ArXiv e-prints  \nAdachi, I., Hayashi, C., & Nakazawa, K. 1976, Progress of Theoretical Physics, 56, 1756 Andrews, S. M., Rosenfeld, K. A., Kraus, A. L., & Wilner, D. J. 2013, ApJ, 771, 129 Andrews, S. M., Wilner, D. J., Hughes, A. M., Qi, C., & Dullemond, C. P. 2009, ApJ, 700, 1502 Ansdell, M., Williams, J. P., Manara, C. F., et al. 2017, AJ, 153, 240 Artymowicz, P. 1993, ApJ, 419, 155 Ataiee, S., Baruteau, C., Alibert, Y., & Benz, W. 2018, A&A, 615, A110 Bai, X.-N. & Stone, J. M. 2010, ApJ, 722, 1437 Baruteau, C., Crida, A., Paardekooper, S.-J., et al. 2014, Protostars and Planets VI, 667  \nBen'ıtez-Llambay, P., Masset, F., Koenigsberger, G., & Szul'agyi, J. 2015, Nature, 520, 63  \nBirnstiel, T., Klahr, H., & Ercolano, B. 2012, A&A, 539, A148  \nFig. A.1. Convergence test for N sp = 25 (red) and 50 (blue). Mass, semi-major axis and eccentricity evolution of the embryo are shown. The line represents the mean value of three individual runs while the area indicates the scattering from these runs. In the bottom panel the RMS eccentricity of the planetesimals are shown in the dashed line.  \n<!-- image -->  \nBitsch, B., Johansen, A., Lambrechts, M., & Morbidelli, A. 2015, A&A, 575, A28  \nBitsch, B., Morbidelli, A., Johansen, A., et al. 2018, A&A, 612, A30  \nBrasser, R., Bitsch, B., & Matsumura, S. 2017, AJ, 153, 222  \nCarrera, D., Gorti, U., Johansen, A., & Davies, M. B. 2017, ApJ, 839, 16  \nCarrera, D., Johansen, A., & Davies, M. B. 2015, A&A, 579, A43  \nChambers, J. E. 1999, MNRAS, 304, 793  \nChatterjee, S. & Tan, J. C. 2014, ApJ, 780, 53  \nChrenko, O., Broˇz, M., & Lambrechts, M. 2017, A&A, 606, A114  \nCresswell, P. & Nelson, R. P. 2008, A&A, 482, 677  \nDominik, C. & Tielens, A. G. G. M. 1997, ApJ, 480, 647  \nDra¸˙zkowska, J. & Alibert, Y. 2017, A&A, 608, A92  \nFlaherty, K. M., Hughes, A. M., Rose, S. C., et al. 2017, ApJ, 843, 150  \nFlaherty, K. M., Hughes, A. M., Rosenfeld, K. A., et al. 2015, ApJ, 813, 99  \nGaraud, P. & Lin, D. N. C. 2007, ApJ, 654, 606  \nGillon, M., Triaud, A. H. M. J., Demory, B.-O., et al. 2017, Nature, 542, 456  \nGoldreich, P., Lithwick, Y., & Sari, R. 2004, ARA&A, 42, 549  \nGoldreich, P. & Tremaine, S. 1979, ApJ, 233, 857  \n- Greenberg, R., Wacker, J. F., Hartmann, W. K., & Chapman, C. R. 1978, Icarus, 35, 1  \nGressel, O., Nelson, R. P., & Turner, N. J. 2011, MNRAS, 415, 3291  \nGuillot, T., Ida, S., & Ormel, C. W. 2014, A&A, 572, A72  \nGuttler, C., Blum, J., Zsom, A., Ormel, C. W., & Dullemond, C. P. 2010, A&A,  \n513, A56  \nHansen, B. M. S. 2009, ApJ, 703, 1131  \nHayashi, C. 1981, Progress of Theoretical Physics Supplement, 70, 35  \n- Hu, X., Tan, J. C., Zhu, Z., et al. 2018, ApJ, 857, 20\n- Ida, S. 1990, Icarus, 88, 129  \nIda, S., Guillot, T., & Morbidelli, A. 2016, A&A, 591, A72  \nIda, S., Kokubo, E., & Makino, J. 1993, MNRAS, 263, 875  \nIda, S. & Lin, D. N. C. 2004, ApJ, 604, 388  \nIda, S. & Lin, D. N. C. 2008, ApJ, 673, 487  \nIda, S. & Makino, J. 1993, Icarus, 106, 210  \nIzidoro, A. & Raymond, S. N. 2018, ArXiv e-prints:1803.08830  \nJohansen, A., Blum, J., Tanaka, H., et al. 2014, Protostars and Planets VI, 547  \n- Johansen, A. & Klahr, H. 2005, ApJ, 634, 1353  \nJohansen, A. & Lambrechts, M. 2017, Annual Review of Earth and Planetary Sciences, 45, 359  \nJohansen, A., Mac Low, M.-M., Lacerda, P., & Bizzarro, M. 2015, Science Advances, 1, 1500109  \nJohansen, A., Oishi, J. S., Mac Low, M.-M., et al. 2007, Nature, 448, 1022  \nJohansen, A., Youdin, A., & Mac Low, M.-M. 2009, ApJ, 704, L75  \n- Johansen, A., Youdin, A. N., & Lithwick, Y. 2012, A&A, 537, A125  \nKataoka, A., Tanaka, H., Okuzumi, S., & Wada, K. 2013, A&A, 557, L4  \nKenyon, S. J. & Bromley, B. C. 2010, ApJS, 188, 242  \nKley, W. & Nelson, R. P. 2012, ARA&A, 50, 211  \nKokubo, E. & Ida, S. 1996, Icarus, 123, 180  \nKokubo, E. & Ida, S. 1998, Icarus, 131, 171  \nKokubo, E. & Ida, S. 2000, Icarus, 143, 15  \nKrijt, S., Ormel, C. W., Dominik, C., & Tielens, A. G. G. M. 2016, A&A, 586, A20  \nLambrechts, M. & Johansen, A. 2012, A&A, 544, A32  \nLambrechts, M. & Johansen, A. 2014, A&A, 572, A107  \nLambrechts, M., Johansen, A., & Morbidelli, A. 2014, A&A, 572, A35  \nLeinhardt, Z. M. & Stewart, S. T. 2012, ApJ, 745, 79  \nLevison, H. F., Kretke, K. A., & Duncan, M. J. 2015, Nature, 524, 322  \nLi, R., Youdin, A. N., & Simon, J. B. 2018, ApJ, 862, 14  \nLissauer, J. J. 1987, Icarus, 69, 249  \nLiu, B. & Ormel, C. W. 2018, A&A, 615, A138  \nManara, C. F., Rosotti, G., Testi, L., et al. 2016, A&A, 591, L3  \nMasset, F. S. 2017, MNRAS, 472, 4204  \nMcNally, C. P., Nelson, R. P., & Paardekooper, S.-J. 2018, MNRAS, 477, 4596  \nMorbidelli, A., Bottke, W. F., Nesvorn'y, D., & Levison, H. F. 2009, Icarus, 204,  \n558  \nMorbidelli, A., Lambrechts, M., Jacobson, S., & Bitsch, B. 2015, Icarus, 258, 418  \nMordasini, C., Alibert, Y., & Benz, W. 2009, A&A, 501, 1139  \nMorishima, R., Schmidt, M. W., Stadel, J., & Moore, B. 2008, ApJ, 685, 1247  \nMustill, A. J., Davies, M. B., & Johansen, A. 2018, MNRAS, 478, 2896  \nOhtsuki, K. & Tanaka, H. 2003, Icarus, 162, 47  \n- Okuzumi, S. & Ormel, C. W. 2013, ApJ, 771, 43  \nOkuzumi, S., Tanaka, H., Kobayashi, H., & Wada, K. 2012, ApJ, 752, 106  \nOrmel, C. W. 2017, in Astrophysics and Space Science Library, Vol. 445,  \nAstrophysics and Space Science Library, ed. M. Pessah & O. Gressel, 197  \n- Ormel, C. W., Dullemond, C. P., & Spaans, M. 2010, ApJ, 714, L103\n- Ormel, C. W. & Klahr, H. H. 2010, A&A, 520, A43\n- Ormel, C. W. & Liu, B. 2018, A&A, 615, A178  \nOrmel, C. W., Liu, B., & Schoonenberg, D. 2017, A&A, 604, A1  \n- Paardekooper, S.-J. 2014, MNRAS, 444, 2031  \nPaardekooper, S.-J., Baruteau, C., Crida, A., & Kley, W. 2010, MNRAS, 401, 1950  \nPaardekooper, S.-J., Baruteau, C., & Kley, W. 2011, MNRAS, 410, 293  \nP'erez, L. M., Chandler, C. J., Isella, A., et al. 2015, ApJ, 813, 41  \n- Pinte, C., Dent, W. R. F., M'enard, F., et al. 2016, ApJ, 816, 25  \nRafikov, R. R. 2004, AJ, 128, 1348  \n- Raymond, S. N., Izidoro, A., Bitsch, B., & Jacobson, S. A. 2016, MNRAS, 458, 2962\n- Raymond, S. N., Kokubo, E., Morbidelli, A., Morishima, R., & Walsh, K. J. 2014, Protostars and Planets VI, 595  \nRicci, L., Testi, L., Natta, A., et al. 2010, A&A, 512, A15  \n- Ros, K. & Johansen, A. 2013, A&A, 552, A137  \nSchafer, U., Yang, C.-C., & Johansen, A. 2017, A&A, 597, A69  \nSchoonenberg, D. & Ormel, C. W. 2017, A&A, 602, A21  \nSchoonenberg, D., Ormel, C. W., & Krijt, S. 2018, submitted  \nSimon, J. B., Armitage, P. J., Li, R., & Youdin, A. N. 2016, ApJ, 822, 55  \nSimon, J. B., Armitage, P. J., Youdin, A. N., & Li, R. 2017, ApJ, 847, L12  \nTazzari, M., Testi, L., Ercolano, B., et al. 2016, A&A, 588, A53  \nThommes, E. W., Duncan, M. J., & Levison, H. F. 2003, Icarus, 161, 431  \nVisser, R. G. & Ormel, C. W. 2016, A&A, 586, A66  \nWeidenschilling, S. J. 1977, MNRAS, 180, 57  \nWeidenschilling, S. J. 2011, Icarus, 214, 671  \nWetherill, G. W. & Stewart, G. R. 1989, Icarus, 77, 330  \nWetherill, G. W. & Stewart, G. R. 1993, Icarus, 106, 190  \n- Yang, C.-C. & Johansen, A. 2014, ApJ, 792, 86\n- Yang, C.-C. & Johansen, A. 2016, ApJS, 224, 39  \nYang, C.-C., Johansen, A., & Carrera, D. 2017, A&A, 606, A80  \n- Yang, C.-C., Mac Low, M.-M., & Johansen, A. 2018, ApJ, 868, 27  \nYang, C.-C., Mac Low, M.-M., & Menou, K. 2012, ApJ, 748, 79  \nYoudin, A. N. & Goodman, J. 2005, ApJ, 620, 459 Youdin, A. N. & Lithwick, Y. 2007, Icarus, 192, 588 Zhu, Z., Stone, J. M., & Bai, X.-N. 2015, ApJ, 801, 81 Zsom, A., Ormel, C. W., Guttler, C., Blum, J., & Dullemond, C. P. 2010, A&A,  \n513, A57"}