Context. Coagulation theory predicts that micronsized dust grains grow into pebbles which drift inward towards the star when they reach sizes of mmcm. When they cross the orbit of a planet a fraction of these drifting pebbles will be accreted. In the pebble accretion mechanism the combined effects of the planets gravitational attraction and gas drag greatly increase the accretion rate. BR  Aims We calculate the pebble accretion efficiency SUB2DSUB  the probability that a pebble is accreted by the planet  in the 2D limit pebbles reside in the midplane. In particular we investigate the dependence of SUB2DSUB on the planet eccentricity and its implications for planet formation models. BR  Methods We conduct Nbody simulations to calculate the pebble accretion efficiency in both the local frame and the global frame. With the global method we investigate the pebble accretion efficiency when the planet is on an eccentric orbit. BR  Results We find that the local and the global methods generally give consistent results. However the global method becomes more accurate when the planet is more massive than a few Earth masses or when the aerodynamic size Stokes number of the pebble is larger than 1. The efficiency increases with the planets eccentricity once the relative velocity between the pebble and the planet is determined by the planets eccentric velocity. At high eccentricities however the relative velocity becomes too high for pebble accretion. The efficiency then drops significantly and the accretion enters the ballistic regime. We present general expressions for SUB2DSUB. Applying the obtained formula to the formation of a secondary planet in resonance with an alreadyformed giant planet we find that the embryo grows quickly due to its higher eccentricity. BR  Conclusions The maximum SUB2DSUB for a planet on an eccentric orbit is several times higher than for a planet on a circular orbit but this increase gives the planet an important headstart and boosts its following mass growth. The recipe for SUB2DSUB that we have obtained is designed to be implemented into Nbody codes to simulate the growth and evolution of planetary systems.
['methods: numerical', 'planets and satellites: formation', 'Astrophysics - Earth and Planetary Astrophysics']
Catching drifting pebbles. I. Enhanced pebble accretion efficiencies for eccentric planets
{'I. Enhanced pebble accretion efficiencies for eccentric planets': 'Beibei Liu 1 , Chris W. Ormel 1  \nAnton Pannekoek Institute (API), University of Amsterdam, Science Park 904,1090GE Amsterdam, The Netherlands e-mail: b.liu@uva.nl,c.w.ormel@uva.nl  \nApril 18, 2018', 'ABSTRACT': "Context. Coagulation theory predicts that micron-sized dust grains grow into pebbles, which drift inward towards the star when they reach sizes of mm -cm. When they cross the orbit of a planet, a fraction of these drifting pebbles will be accreted. In the pebble accretion mechanism, the combined effects of the planet's gravitational attraction and gas drag greatly increase the accretion rate.  \nAims. We calculate the pebble accretion efficiency ε 2D - the probability that a pebble is accreted by the planet - in the 2D limit (pebbles reside in the midplane). In particular, we investigate the dependence of ε 2D on the planet eccentricity and its implications for planet formation models.  \nMethods. We conduct N-body simulations to calculate the pebble accretion efficiency in both the local frame and the global frame. With the global method we investigate the pebble accretion efficiency when the planet is on an eccentric orbit.  \nResults. We find that the local and the global methods generally give consistent results. However, the global method becomes more accurate when the planet is more massive than a few Earth masses or when the aerodynamic size (Stokes number) of the pebble is larger than 1 . The efficiency increases with the planet's eccentricity once the relative velocity between the pebble and the planet is determined by the planet's eccentric velocity. At high eccentricities, however, the relative velocity becomes too high for pebble accretion. The efficiency then drops significantly and the accretion enters the ballistic regime. We present general expressions for ε 2D . Applying the obtained formula to the formation of a secondary planet, in resonance with an already-formed giant planet, we find that the embryo grows quickly due to its higher eccentricity.  \nConclusions. The maximum ε 2D for a planet on an eccentric orbit is several times higher than for a planet on a circular orbit, but this increase gives the planet an important headstart and boosts its following mass growth. The recipe for ε 2D that we have obtained is designed to be implemented into N-body codes to simulate the growth and evolution of planetary systems.  \nKey words. methods: numerical - planets and satellites: formation", '1. Introduction': "In protoplanetary disks, micron-sized dust grains grow by coagulation (Weidenschilling & Cuzzi 1993; Dominik et al. 2007; Johansen et al. 2014). Typically, dust coagulation can be divided into two phases: an initial growth phase, where particles grow from micron-sized dust grains into much larger pebbles; and a drift phase, where particles are transported to the inner disk regions. Aerodynamically, a particle is defined in terms of a dimensionless quantity τ s (often referred to as the Stokes number). Dust grains have τ s glyph[lessmuch] 10 -3 and are strongly coupled to the gas. Particles that have become aerodynamically larger ( τ s ∼ 10 -2 -1 ), on the other hand, decouple from the gas and drift inward (Weidenschilling 1977). The radial drift effectively limits the size of the particles. Other processes that limit the size occur once the particles' relative velocities reach the bouncing or fragmentation thresholds (Güttler et al. 2010). The drifting particles are defined as pebbles, whose size therefore depends on the material strength and the disk properties but typically lies in the mm -cm range (Brauer et al. 2008; Birnstiel et al. 2010). Because the growth timescale of pebbles increases with the disk radius, the dust evolution proceeds in an inside-out manner. The interior, planet-forming regions are therefore constantly supplied by pebbles that drift from the outer regions of disks (Birnstiel  \net al. 2012; Lambrechts & Johansen 2014; Sato et al. 2016; Krijt et al. 2016).  \nObservationally, the existence of pebble-sized particles in disks is inferred from optically thin emission at sub-mm and radio wavelengths. Particles that have grown to sizes larger than the observing wavelength tend to emit as a gray body. Therefore, a millimeter spectral index that is reduced compared to the interstellar medium (ISM) value is a signature of the existence of pebble-sized particles (Natta & Testi 2004; Draine 2006; Pérez et al. 2015; Tazzari et al. 2016). Assuming a temperature (profile) and a millimeter opacity (quantities that are rather uncertain), the total mass of the pebble reservoir can be calculated (Williams & Cieza 2011). Typical values for the pebble disk mass lie from a few Earth masses to hundreds of Earth masses (Ricci et al. 2010; Andrews et al. 2013; Ansdell et al. 2016; Barenfeld et al. 2016; Pascucci et al. 2016). However, the observed pebble mass may not represent the full reservoir of solid material in the disk. For instance, depending on the rapidity of the planet formation process, some of the erstwhile pebbles may already have been locked up in larger bodies, invisible to millimeter disk observations. Indeed, Ansdell et al. (2017) find lower average dust masses for the older σ Orionis cluster.  \nWhen the drifting pebbles cross the orbit of big bodies (planetesimals, planetary embryos, or planets), a fraction of them can be accreted by the combined effects of the planet's grav-  \nitational attraction and gas drag (Ormel & Klahr 2010; Lambrechts & Johansen 2012). This process is known as pebble accretion (see Ormel (2017); Johansen & Lambrechts (2017) for recent reviews). Depending on the importance of gas drag during the pebble-planet encounter, the accretion can be classified into two regimes: settling and ballistic. In the classical, planetesimal-driven planet formation paradigm (Safronov 1972; Lissauer 1987; Ida & Lin 2004), accretion takes place entirely through ballistic interactions. In the ballistic regime, gas drag is unimportant during the encounters and accretion relies on hitting the surface. On the other hand, in the settling regime, pebbles settle onto the planet at their terminal velocity (the velocity where the planet's gravity balances gas drag). The accretion rate does not depend on the planet's physical radius; only the mass of the planet matters. For large planets accreting τ s ∼ 1 particles, the accretion cross section can be as large as the planet's Hill sphere.  \nThe large accretion cross sections and continuous supply of pebbles from the outer disk may, at first sight, offer ideal conditions to produce planets. However, there is one catch: the planet may not accrete all pebbles, simply because they drift too fast or they are stirred to a large height. The viability of pebble accretion as a planet formation process therefore depends not only on the (large) accretion cross sections, but also on the disk conditions, i.e., whether pebbles drift fast and how settled they are. Specifically, we define the pebble accretion efficiency ( ε ) as the number of pebbles that are accreted over the total number of pebbles that the disk supplies. In terms of mass fluxes, the definition reads 1  \nε ≡ ˙ M PA ˙ M disk , (1)  \nwhere ˙ M PA is the pebble accretion rate on the planet and ˙ M disk is the mass flux of pebbles through the disk. Put simply, ε is the probability that a pebble is accreted by the planet. A high value of ε indicates that pebble accretion is efficient. For example, if ε = 1 (the highest value possible) and the disk contains a total 10 M ⊕ in pebbles, an initial 1 M ⊕ planet will grow to 11 M ⊕ , which is large enough to trigger giant planet formation (Pollack et al. 1996). On the other hand, when only one in a thousand pebbles is accreted ( ε = 10 -3 ), the gained planetary mass will be 10 -2 M ⊕ , meaning that the planet growth has significantly stalled. The efficiency ε is a local quantity: it can be computed entirely from the conditions at the location of the planet. Nevertheless, as illustrated by the example, ε carries global significance since it directly states the total amount of pebbles a disk needs to contain in order to grow planets, i.e., how efficiently the pebble mass can be converted into planet mass.  \nThe goal of our paper series is to obtain a general recipe for pebble accretion efficiency ε , i.e., to quantify how ε depends on planet properties (e.g., the planet mass, the eccentricity and the inclination), pebble properties (the pebble's aerodynamical size τ s ), and disk properties (temperature and turbulence strength). In this work (hereafter Paper I) only the planar limit is considered, where we assume that the pebbles have settled to the disk midplane, or more generally, that the pebble accretion radius exceeds the scaleheight of the pebble layer. We elucidate the role of the planet eccentricity on the pebble accretion efficiency in this 2D limit ( ε 2D ). In the subsequent paper (Paper II, Ormel & Liu (2018)), we calculate the 3D pebble accretion efficiency ( ε 3D ) by investigating the roles of the planet inclination and the disk turbulence. The prescriptions that we obtain in these studies can  \nbe implemented in N-body codes, in order to study the formation and long-term dynamical evolution of planetary systems.  \nIn the literature, most theoretical and numerical work consider pebble accretion for planets on circular orbits (Ormel & Klahr 2010; Lambrechts & Johansen 2012, 2014; Guillot et al. 2014; Morbidelli et al. 2015; Bitsch et al. 2015; Levison et al. 2015a; Ida et al. 2016; Matsumura et al. 2017; Picogna et al. 2018). The relatively velocity between the pebble and the planet changes when the planet is on an eccentric orbit. Therefore the accretion efficiency can change as well. During planet formation, planets can acquire moderate eccentricities through mechanisms as planet-planet scattering, secular and mean motion resonances (Lee & Peale 2002; Raymond et al. 2006; Zhou & Lin 2007; Zheng et al. 2017). To investigate how the accretion efficiency depends on the planet eccentricity is therefore relevant. Johansen et al. (2015) already considered the eccentric situation in which they focus on planetesimals with relatively low eccentricities in the local (co-moving) frame. In this paper, we will also conduct numerical orbital integrations in the global (heliocentric) frame, which is especially appropriate for planets with high eccentricities.  \nThe paper is structured as follows. In Sect. 2, we outline two approaches to calculate ε 2D by considering the equation of motion in the local frame and the global frame, respectively. Results are presented in Sect. 3. We compare the results from the above two approaches (Sect. 3.1), and investigate how ε 2D depends on properties of the planet, the pebble and the disk (Sect. 3.2). We also provide analytical fit expressions for ε 2D (Sect. 3.3). In Sect. 4, we apply our results by assessing how fast a secondary planet can grow from a planetary embryo, in the presence of an already-formed giant planet. We summarize our key findings in Sect. 5. The derivation of the pebble accretion efficiency expressions and list of notations are given in the Appendix A.", '2. Method': "In this section, we present two ways to calculate the pebble accretion efficiency. One approach is to treat the pebble's motion with respect to the planet in a non-inertial, local frame (Sect. 2.1). The alternatively approach is to consider it in a global frame centred on the star (Sect. 2.2). Three numerical methods based on the above two approaches are described in Sect. 2.3.", '2.1. Local frame': "In literature, Ormel & Klahr (2010) and Lambrechts & Johansen (2012) adopt a local frame approach to investigate the pebbleplanet interaction. As shown in Fig. 1a, the coordinate is in a comoving frame centred on and rotating with the planet. In the shearing box approximation, the pebble's motion is linearized with respect to the planet (Eq. (7) in Ormel & Klahr (2010)).  \nIn a 2 D limit, the pebble accretion efficiency in the local frame is given by (Guillot et al. 2014; Ormel 2017)  \nε = ˙ M PA ˙ M disk = ∫ | v y ( b ) | ¯ H ( b )Σ P d b 2 πr p v r Σ P , (2)  \nwhere b is the pebble's impact distance measured from the x -axis and v y ( b ) is the pebble's velocity perpendicular to b (Fig. 1a), Ω p is the angular velocity of the planet, and ¯ H ( b ) is the Heaviside step function; ¯ H ( b ) = 1 means the pebble with a impact distance b hits the planet, otherwise, ¯ H ( b ) = 0 . The total pebble mass flux is given by ˙ M disk = 2 πr p v r Σ P , where v r is the radial drift  \nFig. 1. The sketch illustrating pebble accretion in two different frames. a) Local frame: the co-moving frame is centred on the planet where the x-axis is pointing radially away from the star and the y-axis follows the orbital direction of the planet. Grey arrows purely indicate the gas velocity of the shearing flow, and the blue arrow shows the trajectory of the pebble, with the x-axis impact distance b and the y-axis velocity v y ( b ) . b) Global frame: the planet orbits around the central star with a semi-major axis a p . The black arrow depicts the planet motion and the blue arrows illustrate the trajectories of pebbles that cross the orbit of the planet.  \n<!-- image -->  \nvelocity of the pebble and Σ P is the pebble surface density. All above quantities refer to the values at the planet's location ( r p ). In Eq. (2), the efficiency is obtained by the mass flux ratio of the accreted pebbles to total pebbles that cross the orbit of the planet.", '2.2. Global frame': "Here we introduce a new method to calculate the pebble accretion efficiency ε . In the global frame centred on the central star, we integrate the equations of motion to find the orbital evolution of the planet as well as pebbles. As pebbles drift in from regions exterior to the planet and bypass its orbit, the accretion efficiency is simply obtained by counting the fraction of pebbles that hit the planet.  \nThe equation of motion for the pebble in heliocentric coordinates is  \nd 2 r d t 2 = F glyph[star] + F planet + F drag = -GM glyph[star] r r 3 + GM p ( r p -r | r p -r | 3 -r p r 3 p ) -v -v gas t s , (3)  \nwhere G is the gravitational constant, r = ( x, y, z ) , r p = ( x p , y p , z p ) are the positions of the pebble and the planet with respect to the central star, M glyph[star] and M p are the masses of the central star and the planet, v = ( v x , v y , v z ) and v gas = ( v K -v hw ) e φ are the velocity vector of the pebble and the gas's azimuthal velocity at the pebble's location r , and t s is the pebble's stopping time. It should be note that v hw = ηv K is the headwind velocity which measures the deviation between the Keplerian velocity 2  \n( v K ) and the gas azimuthal velocity, η = -1 2 H gas ∂ log P r 2 ∂ log r is the headwind prefactor, H gas is the disk scale height and P is the gas pressure at the planet's location r p .  \nIn Eq. (3), the total (per unit mass) forces of the pebble include the gravitational forces from the central star and the planet, and the drag force from the disk gas.  \nThe pebble is treated as a body with zero gravitational mass. The planet only feels the gravity from the central star, and its equation of motion is expressed as  \nd 2 r p d t 2 = -GM glyph[star] r p r 3 p . (4)  \nThus, we can obtain the orbital motions of the planet and the pebble from Eqs. (3) and (4).  \nIntegrating Eq. (3) in the global frame is computationally expensive, in particular when the pebble is strongly coupled to the gas (small t s ), resulting in a tiny relative velocity ( ∆ v = v -v gas ). In such cases, the pebble tends to move at its terminal velocity, where gas drag balances the gravitational forces. Therefore, the relative acceleration d 2 d t 2 ( r -r gas ) vanishes and the pebble's velocity can be approximated as  \nv ≈ v gas + t s ( F glyph[star] + F planet -Ω 2 gas r ) , (5)  \nWe refer to this approximation as the strong coupling approximation (SCA) (Johansen & Klahr 2005). This condition is satisfied when the pebble and the gas are well coupled ( v ≈ v gas ). It is invalid when the perturbation from the planet is significant or τ s becomes large. Clearly, in Eq. (5) the velocity can be explicitly calculated whereas in Eq. (3) the velocity needs to solve from a differential equation. Therefore, the SCA greatly simplifies the calculation of velocity compared to directly integrating the equation of motion.  \nThe efficiency in the global frame is simply  \nε 2D = N hit N tot , (6)  \nwhere N tot is total number of pebbles across the orbit of the planet, and N hit is the number of pebbles accreted by the planet.", '3. Results': "We employ three methods to calculate the pebble accretion efficiency:  \n- 1. Local - direct integration of the equation of motion (Ormel & Klahr 2010) in the local frame. The pebble accretion efficiency ε 2D is obtained from Eq. (2);\n- 2. Global direct - direct integration of the equation of motion in the global frame (Eq. (3)), and ε 2D is calculated directly from the fraction of pebbles that hit the planet (Eq. (6));\n- 3. Global hybrid - The SCA (Eq. (5)) is applied for | r p -r | > 2 R H ; otherwise switch to direct integration (Eq. (3)). The efficiency ε 2D is also obtained the same as global direct method. The Hill radius is defined as R H ≡ ( M p / 3 M glyph[star] ) 1 / 3 a p and a p is the planet's semi-major axis.  \nWe assume that the inclination of the planet i p is zero. More generally, the inclination can be neglected when it is much smaller than the the scale height of the pebble disk. In this work we only consider the 2 D planar limit where the planet and pebbles are all settled into the disk midplane. 3D effects (planet inclination and disk turbulence) will be modeled in detail in paper II. The local method can be used for a planet on a circular orbit ( e p = 0 ). To model a planet on an eccentric orbit in the local frame, the elliptic motion of the planet needs to be considered, see e.g., Eqs. (15) and (16) in the supplementary material of Johansen et al. (2015). This is nevertheless only a first order approximation, valid for a relatively low eccentricity. The local frame is not a good approach in case that the eccentricity is not so small. In contrast, the global approach is conceptually straightforward to conduct simulations with the planet moving on an eccentric orbit. In this paper, the local method is restricted to planets on circular orbits, while the global method will be applied to planets on both circular and eccentric orbits.  \nIn the local simulation, the pebbles' initial x locations ( x 0 ) range from -5 R H to 5 R H , with an interval of 10 -4 R H , and y 0 is initialized at either -40 R H or 40 R H . The initial radial and azimuthal velocity of the pebble are: v x = -2 τ s v hw / (1 + τ 2 s ) and v y = 1 . 5 x 0 -v hw / (1 + τ 2 s ) (Weidenschilling 1977) where τ s ≡ t s Ω p is the dimensionless stopping time (Stokes number). The simulation terminates when the pebble leaves the domain of the shearing box ( | y | > 1 . 05 y 0 ), or when the separation between the pebble and the planet is smaller than the planet radius R p , indicating a collision.  \nIn the global simulation, the planet is initiated on a circular ( e p = 0 ) orbit (Sect. 3.1) or on an eccentric orbit ( e p > 0 ) with a random true anomaly θ (Sect. 3.2). For circular cases, N tot pebbles are uniformly distributed along a ring exterior to the planet's orbit at an initial distance r 0 = a p (1 + e ) + 5 R H . However, the orbital phase of the planet affects the accretion efficiency when the planet is on an eccentric orbit. In such a situation ( e p > 0 ), pebbles are not only distributed azimuthally, but also radially, from r 0 to r 0 + 2 πv r / Ω K . There are √ N tot grid points in both the vertical and radial directions and N tot in total. The initial radial and azimuthal velocity of the pebble are: v r = -2 τ s v hw / (1 + τ 2 s ) and v φ = v K -v hw / (1 + τ 2 s ) . The simulation terminates when the pebble drifts interior to the orbit of the planet where r < a p (1 -e ) -R H , or collides with the planet where | r -r p | < R p . We set a p = 1 AU for all simulations in this paper.  \nFor both methods, the equations are numerically integrated with a Runge-Kutta-Fehlberg variable timestep scheme (Fehlberg 1969). The relative error tolerance is adopted to be 10 -8 to ensure the numerical convergence.  \nResults comparing the different methods for planets on circular orbits are presented in Sect. 3.1, and global simulations for planets on eccentric orbits are conducted in Sect. 3.2. The analytical fit is given in Sect. 3.3.", '3.1. Comparison between different approaches': "We conduct simulations with the three different methods (Sect. 2.3) for planets on circular orbits ( e p = 0 ) and compare the obtained pebble accretion efficiencies. We assume that the planet density is 3 gcm -3 and central star has 1 M glyph[circledot] . Three planet masses with varying τ s have been explored. Figure 2 shows the pebble accretion efficiency as functions of τ s and M p . The symbols and colors depict results from three different methods, and black lines represent our analytical fit expression (to be discussed in Sect. 3.3). For the global method, the dot and the errorbar represent the mean value from the simulations and the Poisson counting error ( ∆ ε = √ N hit /N tot ), respectively. Hereafter the efficiency refers to its mean value from the simulations. We focus on the pebble accretion regime (settling regime) when 10 -3 glyph[lessorsimilar] τ s < 1 . But for the additional method comparison, the cases of τ s glyph[greaterorsimilar] 1 have also been investigated with the local and the global direct methods, in which the gravity of the central star also becomes important during the pebble-planet encounter (referred to as the three-body regime in Ormel & Klahr (2010)).  \nThe strong coupling approximation (SCA) works well for the small τ s pebble where the gas-pebble coupling time is short. We do not apply the global hybrid method for τ s > 1 cases. It can be seen in Fig. 2 that the hybrid method (orange circle) well matches with the direct method (blue triangle) when τ s glyph[lessorsimilar] 0 . 1 . Moreover, the hybrid algorithm significantly reduces the computational time. For a 0 . 1 M ⊕ planet accreting pebbles of τ s = 10 -2 , the hybrid integration is 40 times faster than the direct integration. It becomes 200 times faster when τ s is 10 -3 .  \nWe also see in Fig. 2 that the results from the local method (red square) and the global methods are in good agreement with each other for the low and intermediate mass planet ( 10 -3 M ⊕ and 0 . 1 M ⊕ ). Only for the massive planet case ( 10 M ⊕ ) the local efficiencies are clearly higher than the global ones.  \nThere are two reasons for the difference in ε 2D between the two methods in the high mass planet case. First, there is a risk of multiple counting accreting pebbles for the local method. In general pebbles can be accreted by the planet through two ways: (i) immediate accretion upon the first penetration of the planet's Hill sphere; and (ii) secondary accretion after a synodical time of the first close-encounter where t syn ∼ Ω -1 p r p / ∆ r co and ∆ r co is the width of the co-orbital region of the planet. In the local frame, we initially place pebbles at both sides of the planet (Sects. 2.1 and 2.3). The problem with this setup, however, is that it may amount to simulate the same pebble twice. The risk of double counting occurs in particular when the radial drift of the pebble within a synodical time is small compared to the accretion radius, v r t syn < b set . It follows that the pebble accretion efficiency of the planet whose mass is higher than a few Earth masses is overestimated due to this multiple counting effect.  \nThe other reason is due to the linearization of the equations of motions. In the local method, we ignore the high order terms of (∆ r p /r p ) n ( n ≥ 2 ) in the pebble's equation of motion by virtue of the approximation that ∆ r p glyph[lessmuch] r p , where ∆ r p is the distance between the planet and the pebble. The integration domain (the maximum ∆ r p we consider in the local frame) is larger for a more massive planet. Thus, the condition ∆ r p glyph[lessmuch] r p  \nFig. 2. Pebble accretion efficiency ε 2D vs τ s and M p . Three planet masses ( M p = 10 -3 M ⊕ , 10 -1 M ⊕ and 10 M ⊕ ) and nine dimensionless stopping times ( τ s ranging from 10 -3 to 10 ) are considered. Three different simulations are conducted for each case, the local method (magenta square), global direct method (blue triangle) and global hybrid method (orange circle), respectively. The errorbar represents the Poission error. Different total number of pebbles are used for different planet mass cases, N tot = 3 × 10 4 for M p = 10 -3 M ⊕ , N tot = 3000 for M p = 10 -1 M ⊕ , and N tot = 300 for M p = 10 M ⊕ , respectively. The solid line corresponds to the analytical fit described in Sect. 3.3. The vertical dashed line gives the boundary where our fit can be applied. For the global approach, the hybrid method matches the direct method well when τ s glyph[lessorsimilar] 0 . 1 . Local simulations are consistent with global simulations in the settling regime ( τ s < 1 ) except when the planet is very massive ( e.g., M p = 10 M ⊕ ).  \n<!-- image -->  \nbreaks down and the first order linear approximation becomes inappropriate for very massive planets. To conclude, both effects are biased towards high planet mass, and therefore the local results overestimate ε in the high mass planet case.  \nWe also see in Fig. 2 that the local and the global simulations are inconsistent when τ s becomes larger than 1 . The two methods give very different results and this mismatch increases with the planet mass. The pebble's eccentricity can be excited by the planet during the encounter. For pebbles with large stopping time, their eccentricity damping from the gas drag is much longer than the encounter time ( glyph[similarequal] orbital time). These pebbles should remain at moderate eccentricities during pebbleplanet interactions. However, in the local method pebbles' are initialized on unperturbed ( e p = 0 ) trajectories, which is not valid when τ s is large. On the other hand, the global method does not suffer from this effect, since it integrates these pebbles for many complete orbits. These differences become more extreme for larger stopping times and higher planet masses. For the M p = 10 M ⊕ and τ s = 10 case, ε 2D evaluates to zero in the global simulation, due to the fact that the pebbles get trapped into resonance with the planet. Of course, resonance trapping cannot be incorporated in a local framework.  \nIn addition, in Fig. 2 the analytical fit expression derived in Sect. 3.3 matches well with both methods when the pebble accretion is in the settling regime ( 10 -3 glyph[lessorsimilar] τ s < 1 ). We find that the efficiency increases with the mass of the planet and decreases with the size of the accreted pebbles. The efficiency is increased, because massive planets have a larger impact cross section Eq. (A.1), while small pebbles, because of their slow radial drift, are more likely to encounter the planet.  \nIn summary, the results of the global methods and the local method are mostly in good agreement with each other in the settling regime. Although computationally more expensive, the global methods are intuitive, and more accurate for higher planet mass ( M p large than a few Earth mass) and large pebbles ( τ s glyph[greaterorsimilar] 1 ). The global method is arguably the better approach to use, in particular for configurations such as resonance trapping,  \nand, more generally, for eccentric orbits. In the following subsection, we only use the global method to calculate the accretion efficiency.", '3.2. Eccentric pebble accretion': 'In Sect. 3.2.1, we describe the default run for eccentric pebble accretion. In Sect. 3.2.2, a parameters study is conducted to explore the dependence of pebble accretion efficiency on the planet mass ( M p ), the stellar mass ( M glyph[star] ), the headwind velocity ( v hw ) and the dimensionless stopping time ( τ s ) with respect to the planet eccentricity.', '3.2.1. Default run': "The parameters of the default model are: M p = 0 . 1 M ⊕ , M glyph[star] = 1 M glyph[circledot] , v hw = 30 m s -1 and τ s = 10 -2 . Figure 3 shows that the pebble accretion efficiency varies with the planet eccentricity. We conduct simulations with 18 different eccentricities, logarithmically spaced between 6 × 10 -5 and 0 . 2 . For each eccentricity, we perform ten simulations randomly varying with different initial planet phase angle θ . In this case we set the total number of pebbles N tot = 200 , N hit is counted from the simulations , and the pebble accretion efficiency is N hit /N tot . The black line is the analytical fit, which will be described in Sect. 3.3.  \nWe find that the pebble accretion efficiency varies with eccentricity. When the planet is on a nearly circular orbit ( e p glyph[lessorsimilar] 10 -3 ), it remains at 0 . 05 . When the eccentricity keeps increasing, the efficiency first increases and then decreases. It attains a maximum value of ε 2D glyph[similarequal] 0 . 18 when the eccentricity approaches 0 . 03 . After that, the efficiency quickly reduces to 0 . 07 when e p glyph[similarequal] 0 . 1 . Clearly, the planet's eccentricity plays an important role in determining the pebble accretion efficiency. At its peak, the efficiency of an eccentric orbit planet ( ε 2D ( e p ) ) is higher than a circular orbit planet by a factor of 4 .  \nFig. 3. Pebble accretion efficiency vs planet eccentricity for the default run: M p = 0 . 1 M ⊕ , M glyph[star] = 1 M glyph[circledot] , v hw = 30 m s -1 , and τ s = 10 -2 . Eccentricities have been explored from e p = 6 × 10 -5 to 0 . 2 . The dot gives the mean value, and the errorbar indicates the Poisson counting error. The solid line represents the analytical fit. The efficiency peaks around an eccentricity of 0 . 03 , and this peak efficiency is four times of the circular efficiency when e p glyph[similarequal] 0 .  \n<!-- image -->  \nBy balancing the settling time and the encounter time, it can be shown that the pebble accretion radius is (Ormel & Klahr 2010; Lambrechts & Johansen 2012)  \nb set ∼ √ GM p t s ∆ v , (7)  \nwhere ∆ v is the relative velocity between the planet and the pebble (see Sect. A).  \nWhen a planet is on a circular orbit, the relative velocity between the pebble and the planet is either dominated by Keplerian shear or the headwind velocity. When a planet is on an eccentric orbit, the relative velocity in addition includes an eccentricity contribution due to the elliptic motion of the planet, which increases with the planet eccentricity. Therefore, when e p is tiny, ∆ v is very close to the circular case and ε 2D still remains constant. The rise of ε 2D is due to an increase in ∆ v when the eccentric velocity ∼ ev K starts to dominate the relative velocity. The flux of pebbles increases. Physically, in the 2D limit, a planet on a more eccentric orbit sweeps up more pebbles due to a larger pebble feeding zone, resulting in a higher efficiency.  \nThe rapid drop of ε 2 D at higher e p glyph[similarequal] 0 . 04 , however, indicates the settling interactions fail. The reason is that the time for the planet-pebble encounter decreases as ∆ v (eccentricity) increases. When the planet-pebble encounter is short, the change of the pebble's velocity due to the gas drag is modest. Therefore, the accretion transitions from the settling regime to the ballistic regime, where the gas drag effect is negligible and the accretion radius reduces to the planet's physical radius (gravitational focusing is unimportant at these high eccentricities). Thus, the accretion efficiency declines significantly during this transition.  \nIn order to analyse this issue in detail, we investigate at which location (planet's orbital phase angle) pebbles are captured by the planet. The capture condition is defined as the separation between the pebble and the planet is smaller than the accretion radius ( | r -r p | < b set ). In order to clarify the effect we focus on the low mass planet when it is initially in the headwind regime, therefore we adopt: M p = 0 . 01 M ⊕ , M glyph[star] = 1 M glyph[circledot] ,  \nτ s = 10 -2 , v hw = 60 m / s . We artificially set the planet radius 10 times smaller than its normal value to eliminate any ballistic accretion. Figure 4 plots the relative velocity between the planet and pebble (left y axis) and the number of accreted pebbles (right y axis) as functions of the planet's phase angle for e p = 7 × 10 -3 (Fig. 4a) and e p = 3 × 10 -2 (Fig. 4b). The radial and azimuthal velocity of a planet on eccentric orbit are v K ( r p ) e p sin θ/ √ 1 -e 2 p and v K ( r p )(1 + e p cos θ ) / √ 1 -e 2 p , where θ is the true anomaly of the planet (Murray & Dermott 1999). Assuming that the orbit of the pebble is not affected by the planet, its radial and azimuthal velocity are -2 v hw τ s v K / (1+ τ 2 s ) and v K -v hw / (1 + τ 2 s ) . Therefore, the relative velocity between the planet and the pebble can be analytically calculated, which is shown by the blue lines of Fig. 4.  \nWe find that in Fig. 4 the number of accreted pebbles correlates with ∆ v . In the low eccentricity case (Fig. 4a), more pebbles are accreted when ∆ v is higher. For instance, when the pebble moves on a sub-Keplerian orbit it has the lowest relative velocity with the planet at aphelion ( 180 · ). We find that at this location the planet accretes the least number of pebbles. Since in Eq. (7) b set ∝ ∆ v -1 / 2 , the accretion radius is largest at aphelion. However, the accretion rate correlates with b set ∆ v in the 2D limit and with b 2 set ∆ v in the 3D limit (Johansen et al. 2015; Morbidelli et al. 2015). Therefore, ε 2D increases as ∆ v 1 / 2 . Less pebbles are accreted at aphelion than elsewhere.  \nOn the other hand, when the eccentricity is high, the pebble accretion rate displays the opposite dependence on the phase angle (Fig. 4b), i.e., more pebbles are accreted at lowest ∆ v (aphelion). This is because in this high eccentricity case, ∆ v becomes too high for the settling. Thus, only pebbles with relatively low ∆ v are able to be accreted by the planet. The transition velocity from the settling to the ballistic accretion is marked by the black dashed line in Fig. 4 (Eq. (13)).", '3.2.2. Parameter study': 'To investigate the effects of different parameters on the pebble accretion efficiency, we conduct a parameter study by varying the planet mass ( M p ), the stellar mass ( M glyph[star] ), the headwind velocity ( v hw ) and the dimensionless stopping time ( τ s ). We explore four planet masses ( 10 -3 M ⊕ , 10 -2 M ⊕ , 10 -1 M ⊕ , 1 M ⊕ ), three stellar mass ( 3 M glyph[circledot] , 1 M glyph[circledot] , 0 . 3 M glyph[circledot] ), three headwind velocities ( 15 m s -1 , 30 m s -1 , 60 m s -1 ) and three dimensionless stopping times ( 10 -3 , 10 -2 , 10 -1 ). Nine of these simulations are illustrated in Table 1 and Fig. 5.', 'Planet mass': "Fig. 5a shows the efficiency as a function of the eccentricity for three different planet masses. The other parameters are the same as the default values. The default case ( M p = 0 . 1 M ⊕ ), the massive planet case ( M p = 1 M ⊕ ), and the less massive planet case ( M p = 0 . 01 M ⊕ ) are shown in black, dark red and light red, respectively.  \nThe three cases exhibit a similar trend with eccentricity. The efficiency curves are initially flat, gradually increase with e p later and drop to lower values. For the massive planet case, the efficiency is 0 . 22 when the planet is on a circular orbit. The maximum efficiency attains 0 . 59 when e p = 0 . 05 . For the less massive planet case, the efficiency of a circular orbit planet is 0 . 013 and the maximum efficiency is 0 . 04 when e p is close to 0 . 01 . Comparing the three cases, we find that (i) the efficiency increases with the planet mass for all eccentricities; more massive  \nFig. 5c plots three cases with different headwind velocities. The disk headwind now varies with low speed ( v hw = 15 m s -1 ) and high speed ( v hw = 60 m s -1 ) shown in dark blue and light blue, respectively. Efficiencies obtained from the default run ( v hw = 30 m s -1 ) are shown in black.  \n<!-- image -->  \n0  \n90  \n180  \n270  \n360  \nPhase angle (degree)  \nFig. 4. Relative velocity between the planet and pebbles (dark blue line) and number of accreted pebbles (light blue bars) as function of the planet's orbital phase angle. Perihelion is at 0 · and aphelion is at 180 · . The planet eccentricity is e p = 0 . 007 in the left and e p = 0 . 03 in the right panel, respectively. Other parameters are the same for the two panels: M p = 0 . 01 M ⊕ , M glyph[star] = 1 M glyph[circledot] , τ s = 10 -2 , v hw = 60 m / s and N tot = 300 × 300 . The black dashed line represents the transition velocity from the settling to the ballistic accretion, v = v ∗ (Eq. (13)).  \nplanets accrete pebbles more efficiently (the same as the circular case in Sect. 3.1). (ii) The transition eccentricity where the efficiency attains its maximum increases with the planet mass. This is because the transition occurs when the relative velocity ∆ v is so high that the settling accretion enters the ballistic accretion (see Sect. A.1). Since a more massive planet has a stronger gravitational attraction that enables to accrete more eccentric pebbles, the transition velocity increases with the planet mass (Eq. (13)). Consequently, the transition eccentricity also increases with the planet mass.", 'Stellar mass': 'Fig. 5b illustrates the dependence of the stellar mass on the efficiency where the rest of parameters are identical as the default case. The default case ( M glyph[star] = 1 M glyph[circledot] ) is shown in black, and the massive giant star case ( M glyph[star] = 3 M glyph[circledot] ) and the low-mass M dwarf star case ( M glyph[star] = 0 . 3 M glyph[circledot] ) are given in light purple and dark purple, respectively.  \nFor the massive star, the efficiency is 0 . 026 when the planet is on a circular orbit. The maximum efficiency attains 0 . 09 when e p = 0 . 03 . For the less massive star, the efficiency for a circular orbit planet is 0 . 1 and the maximum efficiency is 0 . 34 when e p = 0 . 04 . Comparing the three cases, we find that (i) the efficiency decreases with the stellar mass. (ii) The transition eccentricity decreases with the stellar mass. Both effects can be explained from the fact that a planet around a less massive star has a larger Hill sphere and therefore can accrete more pebbles.', 'Headwind velocity': 'For the low headwind case, the efficiency of pebble accretion is 0 . 10 when the planet is on a circular orbit and the maximum efficiency is 0 . 33 . For higher headwind velocity case, the efficiency when the planet is on a circular orbit is 0 . 03 , and the  \nmaximum efficiency is 0 . 08 . We conclude that (i) the efficiency decreases with the headwind velocity. Because the headwind is related with the radial drift velocity of the pebble, a lower headwind speed results in a slower radial drift, and the planet is therefore able to accrete more pebbles. Furthermore, (ii) the eccentricity for which the maximum efficiency is attained, is independent of the headwind velocity. This is because as the eccentricity increase, ∆ v is dominated by the eccentricity velocity but not the headwind velocity. Therefore, the transition velocity between the settling and the ballistic regime is independent of the headwind velocity (Eq. (13)).', 'Dimensionless stopping time': 'Finally, three cases with different τ s are shown in Fig. 5d. The default case ( τ s = 10 -2 ) is shown in black while small pebbles ( τ s = 10 -3 ) and large pebbles ( τ s = 10 -1 ) are shown in dark green and light green, respectively.  \nFor τ s = 10 -3 the efficiency of pebble accretion is 0 . 13 when the planet is on a circular orbit. When e p glyph[similarequal] 0 . 06 , the efficiency attains a maximum value of 0 . 60 . For τ s = 10 -1 the efficiency of a circular orbit planet is 0 . 02 . The maximum efficiency is close to 0 . 04 when e p = 0 . 01 . It is clear that the planet is more efficient at accreting small size pebbles due to their slow radial drift, which is consistent with the results of circular pebble accretion in Sect. 3.1.  \nIn addition, we find that the ratio of the maximum ε 2D to ε 2D when e p glyph[similarequal] 0 is high for τ s = 10 -3 and low for τ s = 10 -1 . The reason is that the settling interactions extend to higher e p when τ s is small. Therefore, the transition to the ballistic regime occurs at higher e p and the peak efficiency has a stronger boost for small τ s .  \nTo summarize, the pebble accretion efficiency ε 2D is an increasing function of the planet-to-star mass ratio ( q p = M p /M glyph[star] ), a decreasing function of the headwind velocity ( v hw ) and the dimensionless stopping time of the pebble ( τ s ). The efficiency is independent of eccentricity when it is relatively small ( e p glyph[lessorsimilar] 10 -3 ), gradually increases with eccentricity when it is moderate ( e p glyph[similarequal] 10 -2 ), and drops quickly when the eccentricity becomes relatively high ( e p glyph[greaterorsimilar] 0 . 1 ).  \nb  \n)  \ne  \np = 3  \n×  \n10  \n2  \n-  \n)  \nK  \nv  \n(  \ny  \nt  \ni  \nc  \no  \nl  \ne  \nv  \ne  \nv  \ni  \nt  \na  \nl  \ne  \nR  \n0  \n.  \n03  \n0  \n.  \n02  \n0  \n.  \n01  \n0  \n.  \n005  \n120  \n80  \n40  \ns  \nt  \ni  \nh  \nf  \no  \nr  \ne  \nb  \nm  \nu  \nN  \nD  \n2  \nε  \ny  \nc  \nn  \ne  \ni  \nc  \ni  \nf  \nf  \ne  \nn  \no  \ni  \nt  \ne  \nr  \nc  \nc  \na  \ne  \nl  \nb  \nb  \ne  \nP  \n0  \n.  \n01  \n0  \n.  \n5  \n0  \n.  \n1  \n(a)  \n10  \n4  \n-  \nM  \np = 1  \nM  \nM  \np = 0  \n.  \n1  \nM  \nM  \np = 0  \n.  \n01  \nM  \n10  \n⊕  \n⊕  \n3  \n-  \n10  \n-  \nPlanet eccentricity  \ne  \n2  \np  \nFig. 5. Pebble accretion efficiency vs planet eccentricity. The default run is shown in black, with M p = 0 . 1 M ⊕ , M glyph[star] = 1 M glyph[circledot] , v hw = 30 m s -1 and τ s = 10 -2 . In each panel one parameter varies and the other three keep the same. The upper left panel: three different planet masses, M p = 1 M ⊕ , M p = 0 . 1 M ⊕ and M p = 0 . 01 M ⊕ ; The upper right panel: three different stellar masses, M glyph[star] = 3 M glyph[circledot] , M glyph[star] = 1 M glyph[circledot] and M glyph[star] = 0 . 3 M glyph[circledot] ; The lower left: three different headwind velocities, v hw = 15 m s -1 , v hw = 30 m s -1 and v hw = 60 m s -1 ; The lower right: three dimensionless stopping times, τ s = 10 -3 , τ s = 10 -2 and τ s = 10 -1 . The dot gives the mean value, and the errorbar indicates the Poisson counting error. The solid line is the analytical fit of Eq. (16). The efficiency increases with M p , and decreases with M glyph[star] , v hw and τ s .  \n<!-- image -->  \nM  \n=0  \n.  \n3  \nM  \nM  \n=1  \nM  \nM  \n=3  \nM  \nτ  \ns = 10  \nτ  \ns = 10  \nτ  \ns = 10  \n-  \n3  \n2  \n-  \n1  \n-  \n10  \n/circledot  \n/circledot  \n10  \n-  \n-  \n/circledot  \n3  \n10  \n-  \nPlanet eccentricity  \ne  \n3  \n10  \n-  \nPlanet eccentricity  \ne  \n2  \np', '3.3. Analytical fit expression': "In this subsection, we present analytical fitting formulas, the detailed derivations of these formulas are described at length in the appendix.  \nWe find that the obtained accretion efficiency in the settling regime can be well fitted by  \nε set , 2D = 0 . 32 √ q p τ s η 2 ( ∆ v v K ) f set , (8)  \nwhere q p is the mass ratio of the planet to the central star. The relative velocity between the planet and the pebble ( ∆ v ) is  \n∆ v = max( v cir , v ecc ) , (9)  \nHere v cir is the relative velocity between the circular orbit planet and the pebble. We adopt an expression of v cir that combines the  \nArticle number, page 8 of 15  \nheadwind and the shear regimes,  \nv cir = [ 1 + 5 . 7 ( q p q hw / sh )] -1 v hw + v sh , (10)  \nwhere v hw = ηv K , v sh = 0 . 52( q p τ s ) 1 / 3 v K and the transition mass ratio for two regimes q hw / sh = η 3 /τ s . The eccentric velocity of the planet relative to its circular Keplerian value is v ecc , where we fit  \nv ecc = 0 . 76 e p v K . (11)  \nIn the epicycle approximation, it can be shown that the velocity relative to a circular orbit ranges from a minimum of e p v K / 2 to maximum of e p v K (e.g. Johansen et al. 2015). Our numerical coefficient of 0 . 76 falls between these limits. Once η > e p , ∆ v will be dominated by the headwind, also consistent with the analytical expression by (Johansen et al. 2015).  \nWhen ∆ v becomes higher, the accretion is no longer in the settling regime. We adopt an exponential function f set to express  \n⊕  \n0  \n.  \n5  \n0  \n.  \n1  \n10  \n1  \n-  \nD  \n2  \nε  \ny  \nc  \nn  \ne  \ni  \nc  \ni  \nf  \nf  \ne  \nn  \no  \ni  \nt  \ne  \nr  \nc  \nc  \na  \ne  \nl  \nb  \nb  \ne  \nP  \nD  \n2  \nε  \ny  \nc  \nn  \ne  \ni  \nc  \ni  \nf  \nf  \ne  \nn  \no  \ni  \nt  \ne  \nr  \nc  \nc  \na  \ne  \nl  \nb  \nb  \ne  \nP  \n0  \n.  \n01  \n0  \n.  \n01  \n0  \n.  \n5  \n0  \n.  \n1  \n(b)  \n10  \n(d)  \n10  \n-  \n4  \n-  \n4  \n2  \np  \n10  \n10  \n-  \n1  \n-  \n1  \nTable 1. Model set-up in a parameter study (planet density is 3 gcm -3 )  \nthe decay of ε set , 2D ,  \nf set = exp [ -0 . 5 ( ∆ v v ∗ ) 2 ] , (12)  \nand  \nv ∗ = ( q p /τ s ) 1 / 3 v K (13)  \nis the transition velocity from the settling to the ballistic regimes (see derivations in Sect. A.1).  \nThe accretion efficiency in the ballistic regime is  \nε bal , 2D = R p 2 πτ s ηr p √ 2 q p r p R p + ( ∆ v v K ) 2 (1 -f set ) . (14)  \nIt is important to note that ε set , 2D is independent of r p and R p , whereas ε bal , 2D is related with the ratio of above two quantities, R p /r p .  \nThe total accretion efficiency is  \nε 2D = ε set , 2D + ε bal , 2D (15)  \nThe above recipe for ε 2D is calculated based on rate expression. It is therefore not guaranteed that the accretion probability ε 2D remains no larger than unity ( ε 2D glyph[lessorequalslant] 1 ). The situation is analogous to radioactive decay, where, for a decay rate of λ , the probability of decaying after a time is P = 1 -exp( -λt ) . We can therefore correct for this effect, simply by redefining ˜ ε 2D = 1 -exp( -ε 2D ) Clearly, in such a expression, ˜ ε 2D = ε 2D when ε 2D glyph[lessmuch] 1 and ˜ ε 2D glyph[similarequal] 1 when ε 2D glyph[greaterorequalslant] 1 . Accounting for these probabilistic nature of ε , we give the efficiency expression as,  \nε 2D → 1 -exp( -ε 2D ) . (16)  \nto ensure that ε ≤ 1 . In Fig. 2 and Fig. 5, we find that this analytical fit expression (Eq. (16)) agrees with the simulations quite well for planets on circular orbits as well as on eccentric orbits.  \nLambrechts & Johansen (2014) also calculate the efficiency when the planet is on a circular orbit. In the shear-dominated regime, Eq. (8) can be written as  \nε sh = 0 . 1 ( M p M ⊕ ) 2 / 3 ( η 10 -3 ) -1 ( τ s 0 . 1 ) -1 / 3 glyph[similarequal] 0 . 022 ( M p M ⊕ ) 2 / 3 ( τ s 0 . 1 ) -1 / 3 ( r 10 AU ) -1 / 2 (17)  \nThe latter expression of the above formula adopts the same disk model as Lambrechts & Johansen (2014). Comparing Eq. (17)  \nwith Lambrechts & Johansen (2014)'s Eq.(33), we find that the scaling relations are identical, and our prefactor is 35% lower than theirs.  \nIn the headwind-dominated regime, Eq. (8) can be written as  \nε hw = 0 . 055 ( M p M ⊕ ) 1 / 2 ( η 10 -3 ) -1 / 2 ( τ s 0 . 1 ) -1 / 2 . (18)  \nThe transition mass from the headwind to the shear regimes is  \nM hw / sh = 0 . 03 ( η 10 -3 ) 3 ( τ s 0 . 1 ) -1 M ⊕ . (19)  \nIn the eccentricity-dominated regime, Eq. (8) can be written as (neglecting f set )  \nε ecc = 0 . 055 ( e p 10 -3 ) 1 / 2 ( M p M ⊕ ) 1 / 2 ( η 10 -3 ) -1 ( τ s 0 . 1 ) -1 / 2 . (20)  \nWe find that in this regime ε set , 2D increases as e 1 / 2 p .", '3.4. Neglected effects': "In deriving the expression above, we have neglected several effects:  \n- 1. Evaporation/molten of pebbles for planetesimals travelling on supersonic velocity (relevant for e p > h gas ). One caveat is that planetesimals on eccentric orbits can produce bow shocks as they move supersonically through the disk gas (i.e., e p > h gas , Morris et al. (2012)). The surrounding gas temperature is raised due to energy deposition at the shock front. Therefore, pebbles are likely to be sublimated/melted before they settle onto the planetesimal, depending on the gas density and pressure (e.g., Fig. 6 of Brouwers et al. (2017)), the planetesimal's mass, and the pebble's size and composition (Hood 1998; Miura & Nakamoto 2005). The detail treatment is beyond the scope of this paper; here we do not solve for the thermal balance of the pebble. Nevertheless, stopping the mass growth due to the pebble evaporation is most relevant for low-mass planetesimals. For high mass planets capable of accreting the primordial disk gas, evaporating pebbles may not reach the surfaces of the cores, but still enrich their envelopes (Venturini et al. 2016; Alibert 2017; Brouwers et al. 2017).  \n- 2. Aerodynamic deflection (Guillot et al. 2014; Johansen et al. 2015; Visser & Ormel 2016), relevant for very low-mass planetesimals and small τ s pebbles. When the planet mass is very low, its impact radius reduces to the physical radius. However, when the stopping time of the pebble is short compared to the crossing time of the planetesimal ( ∼ R p /v hw ), the pebble will follow gas streamlines and avoid accretion.This occurs for typical planetesimal size R p glyph[lessorsimilar] 100 km and τ s glyph[lessorsimilar] 10 -3 .\n- 3. Pre-planetary atmosphere formation, relevant for b set < R Bondi . Once the planet's surface escape velocity becomes larger than the local sound speed, a planet will start to accrete the disk gas, creating a pre-planetary atmosphere. The gas density in planetary atmosphere is higher compared to the surrounding disk gas. When the pebble enters the planetary atmosphere, gas drag becomes larger and the accretion cross section also increases due to this enhanced drag. However, for pebble accretion, the accretion cross section ( b set glyph[similarequal] τ 1 / 3 s R H in the shear regime) can be as large as the planet's Hill sphere. The radius of the planetary atmosphere cannot exceed the Bondi radius, R Bondi = GM p /c 2 s , where c s is the sound speed at the planet's location. As long as R Bondi < b set (or q p < 2 . 3 τ 1 / 2 s h 3 gas in dimensionless units), it is therefore appropriate to neglect the atmosphere enhancement on pebble accretion. For a planet accreting τ s = 10 -2 pebbles in the disk with a typical scale height h gas = 0 . 05 , this condition is justified as long as the planet is less than 10 M ⊕ .  \nIn addition, we assume that the gas moves on an unperturbed sub-Keplerian orbit. In reality the planet will perturb the gas flow. A natural scale on which these effects appear is again the planet's Bondi radius (Ormel et al. 2015). Therefore, the flow pattern effect on pebble accretion is also minor when R Bondi < b set .  \n- 4. Pebble isolation mass. When the planet is massive enough to strongly perturb the disk gas, a gap can be formed in the vicinity of the planet and inverses the locally gas pressure gradient. This process truncates pebbles' inward drift, and therefore terminates pebble accretion. The planet mass is defined as the pebble isolation mass ( M p , iso ). Approximating the gap opening mass mass (Lin & Papaloizou 1993), M p , iso is around 20 M ⊕ for a solar mass star ( M p , iso ∼ h 3 gas M glyph[star] , Lambrechts et al. (2014); Bitsch et al. (2018)). In our study, we consider the planet mass in a range of 10 -3 M ⊕ glyph[lessorsimilar] M p glyph[lessorsimilar] 10 M ⊕ < M p , iso .", '4. Application: formation of a secondary planet': "In this section we envision an already-formed primary planet (a gas giant) and a planetary embryo for the future second planet. We discuss how secondary planet formation is aided by a higher pebble accretion efficiency at resonance locations, where eccentricities are excited by the primary giant planet.  \nBased on the observational statistics of the occurrence rate, 10% -20% of solar-type systems contain gas giants (Cumming et al. 2008; Mayor et al. 2011). The formation of a gas giant requires a solid core mass that exceeds 10 M ⊕ (Pollack et al. 1996) before the gas disk is depleted. Disk migration is an efficient way to transport embryos to construct such cores, but it also depends on detailed disk models (Cossou et al. 2014; Coleman & Nelson 2014; Liu et al. 2015). On the other hand, (sub)millimeter dust observations suggest the existence of massive pebble reservoirs ( ∼ 100 M ⊕ ) in the outer region of disks  \nduring the gas-rich phase (Ricci et al. 2010; Andrews et al. 2013; Ansdell et al. 2016). In the pebble accretion scenario the inward drift of these pebbles provide the building blocks to form massive planets (Lambrechts & Johansen 2014; Bitsch et al. 2015). Here we only focus on the pebble accretion scenario model.  \nOnce the core mass reaches the 10 M ⊕ , the gas accretion enters a runaway mode (Pollack et al. 1996). As a result, the surrounding embryos and planetesimals are strongly perturbed by the sudden increase of planet's gravity (Zhou & Lin 2007; Raymond & Izidoro 2017). These bodies could be scattered outward during close encounters. They subsequently migrate inward either due to the type I migration (embryos) or aerodynamic gas drag (planetesimal), and could be capture into the 2 : 1 mean motion resonance with this giant planet (Zhou & Lin 2007). Here, we only focus on the largest embryo among these bodies because it has the highest chance to grow into a secondary planet. The embryo's eccentricity will be excited at the resonant location. Meanwhile, it can accrete pebbles that drift from the outer part of the disk. With the pebble accretion prescription of Sect. 3.3, we conduct N-body simulations to simulate the embryo's orbital and mass evolution. We explore how massive the pebble disk must be to produce a secondary planet.  \nWe restrict our problem to a quiescent (low turbulent) disk such that pebbles are settled into the disk midplane. The accretion is thus in the 2 D regime (pebble scale height H P is smaller than impact radius b set ). The influence of disk turbulence and the prescription of 3 Dpebble accretion will be presented in Paper II. We use the pebble generation model derived by Lambrechts & Johansen (2014). The adopted gas surface density, the disk aspect ratio are (Lambrechts & Johansen 2014)  \nΣ gas = 500 r -1 AU gcm -2 , h gas = 0 . 033 r 1 / 4 AU , (21)  \nwhere r AU = r/ (1 AU) , the headwind prefactor η = 1 . 5 × 10 -3 r 1 / 2 AU and the dust-to-gas ratio is 1% .  \nThe formation of the first gas giant may already take a few Myr, depending on the opacity in its envelope (Movshovitz et al. 2010). Thus, ˙ M disk = 5 × 10 -5 M ⊕ / yr is the pebble mass flux estimated from Lambrechts & Johansen (2014) (assuming t = 5 × 10 6 yr in their Eq. (14)). The pebble aerodynamal size is τ s = 0 . 05 based on their Eqs. (20) and (25). The above τ s is consistent with advanced dust coagulation model (Birnstiel et al. 2012) as well as disk observations (Pérez et al. 2015; Tazzari et al. 2016).  \nA gas giant of 1 M J is initialized at 1 AU . We ignore the migration of the gas giant, but do consider the eccentricity damping of the disk gas, which operates on a timescale of τ e = 10 3 yr (Bitsch et al. 2013). Since the Jupiter mass planet is much more massive than the embryo, the eccentricity of the gas giant remains low even without any gas damping. We find that the specific choice of τ e for the giant planet does not affect our simulation results. A 0 . 1 M ⊕ embryo is initialized at a period ratio of 2 . 1 (just exterior to the 2 : 1 resonance) with the inner giant planet.  \nWe consider two scenarios: (i) a resonant case where the embryo experiences type I migration and is captured in resonance; and (ii) a non-resonant case without type I migration but including eccentricity damping for a comparison. For the resonant case, we implement the additional accelerations of the planetdisk interaction and eccentricity damping into the N-body code based on Papaloizou & Larwood (2000) and Cresswell & Nelson (2008): a m = -v p /τ m and a e = -2( v p · r p ) r p /r 2 p τ e . The embryo's migration timescale τ m and the eccentricity damping timescale τ e are adopted from Eqs. (11) and (13) in Cresswell  \nFig. 6 shows a) the eccentricity, b) the mass growth, and c) the pebble accretion efficiency of the embryo as functions of the total mass of drifting pebbles (lower x-axis) and time (upper xaxis). For the resonant case, the embryo gets trapped into the 2 : 1 resonance where the eccentricity is gradually excited (black line in Fig. 6a). Meanwhile, the gas damps the eccentricity and circularizes the orbit. An equilibrium eccentricity is achieved by balancing gas damping and resonant excitation. The equilib-  \n<!-- image -->  \n⊕  \nFig. 6. Eccentricity (a), pebble accretion efficiency (b) and mass of the embryo (c) vs total mass of drifting pebbles (lower x axis) and time (upper x axis). The black line corresponds to the resonant case with type I migration and resonance trapping, while the red line represents the non-resonant, comparison case that the embryo without type I migration moves on a nearly circular orbit.  \n&Nelson (2008). These expressions take into account the supersonic regime when e p glyph[greaterorsimilar] h gas . Thus, for the resonant case we expect that the embryo will migrate and get trapped into the 2 : 1 resonance while for the comparison case the embryo will remain at the original orbit. It should be noted that we do not consider the inclinations of both planets. Since the inclination damping timescale is much shorter compared to the migration timescale (Bitsch & Kley 2011), their inclinations would be soon damped down, restricting both planets to coplanar orbits. Thus, the above simplification is appropriate.  \nThe gas giant also opens a gap in the disk gas with a typical width of Hill radius. The depletion of the gas near the planet also changes the gas pressure gradient, resulting in the truncation of drifting pebbles. Therefore, a wider dust gap can be formed exterior to the gas gap. Hydrodynamic simulations show that the width of the dust gap produced by a Jupiter mass planet is less than the distance between the giant planet and its 2 : 1 resonance (Zhu et al. 2012). Therefore, embryos at the 2 : 1 resonance can still accrete pebbles that drift from the outer region of the disk. 2  \nrium eccentricity ( e eq ) can be analytically derived by solving Lagrange's planetary equations with migration and eccentricity damping (see details in Teyssandier & Terquem (2014)). Due to the high mass ratio of the giant planet and the embryo, from their Eq. (46) e eq can be simplified as e eq glyph[similarequal] √ τ e /τ m glyph[similarequal] h gas . In Fig. 6a we also find e eq glyph[similarequal] 0 . 035 , consistent with the above analysis. However, the embryo in the comparison case without migration has a nearly circular orbit ( e p < 4 × 10 -3 ; red line in Fig. 6a). It is non-zero due to weak, secular perturbations of the inner gas giant.  \nWe find in Fig. 6b that after 3 × 10 4 yr the pebble accretion efficiency in the resonant case is higher than the non-resonant, comparison case. As was already seen from Fig. 3, this is again because for the resonant case the moderate eccentricity in the 2 : 1 resonance facilitates a higher pebble accretion rate. After that, a higher eccentricity in the resonant case increases the pebble accretion efficiency, leading to a more massive planet. Because of its higher mass, it accretes pebbles more efficiently, promoting an even higher pebble accretion efficiency. That is why in Fig. 6b both efficiencies increase with time, but the efficiency in the resonant case (black) is higher than the comparison case (red). As a result, in Fig. 6c the resonant embryo grows more rapidly than the non-resonant planet. And this mass difference increases with time.  \nIn Fig. 6c we see that the final mass of the embryo depends on the total amount of pebbles available in the outer disk. Since the formation of the giant planet already consumes large amounts of pebbles, the amount of pebbles left behind that feed the secondary planet should be limited. For instance, if the total mass of drifting pebbles beyond the planet's orbit is limited to 20 M ⊕ , the embryo at resonance can grow to a 3 M ⊕ super-Earth, whereas the embryo in a non-resonant orbit will only reach 1 M ⊕ . Similarly, if the total pebble mass is limited to 40 M ⊕ , the embryo in resonance can grow into a planet with 11 M ⊕ , whereas the non-resonant embryo only attains 3 M ⊕ . For the resonant case, then, the massive core is able to initiate the rapid gas accretion and grow into a giant planet. However, for the non-resonant case, the slow growth of the embryo results in a final super-Earth mass planet.  \nTo conclude, an embryo at resonance accretes pebbles more efficiently and grows faster than neighboring embryos that move on nearly circular orbits. A secondary planet therefore preferentially forms at the resonance location. Radial velocity surveys show that multiple gas giants pile up the 2 : 1 resonance (Wright et al. 2011), supporting our hypothesis.", '5. Conclusions': "Pebble accretion is an important mechanism that drives planet growth by accreting millimeter to centimeter size particles in gas-rich disks. Formed in the outer regions of disks, pebbles drift inward due to the aerodynamic drag from the disk gas. When pebbles cross the orbit of a planet, a fraction of them will be accreted by the planet. In this paper, we have calculated the efficiency of pebble accretion in the 2 Dlimit by conducting a series of numerical integrations of the pebble's equation of motion in both the local (co-moving) frame and the global (heliocentric) frame.  \nThe key findings of this study are the following:  \n- 1. The results of the local and global methods are generally consistent. However, the global method more accurately simulates the pebble-planet interaction, due to the fact that it accounts for curvature effects and models the dynamics of  \nthe pebble properly. Since the equations of motions are linearized in the local frame, the local method tends to overestimate ε when the planet is more massive than a few Earth masses, or when the aerodynamic size of the pebble is larger than 1 (Sect. 3.1).  \n- 2. We find that the 2D efficiency ( ε 2D ) is a function of the planet eccentricity. A planet will accrete pebbles at a higher efficiency once the eccentricity velocity is higher than the relative velocity obtained for a circular orbit planet ( v ecc > v cir ). The pebble accretion efficiency then increases with eccentricity. However, ε 2D drops quickly when the eccentricity becomes too large for encounters to satisfy the settling conditions. The accretion therefore transitions to the ballistic regime. Planets with moderate eccentricities ( e p ∼ 10 -2 -0 . 1 ) accrete pebbles at rates a factor of 3 -5 higher than planets on circular orbits (Sect. 3.2).\n- 3. We have obtained a recipe for the pebble accretion efficiency ε 2D as functions of the planet eccentricity ( e p ), the mass ratio of the planet to the star ( q p ), the disk headwind prefactor ( η ), and the aerodynamic size of the pebble ( τ s ). Consistent with previous work, the efficiency increases with the planetto-star mass ratio, and decreases with both the headwind velocity and the pebble size. An analytical fit expression of ε 2D ( e p , q p , η, τ s ) is derived from our simulations (Sect. 3.3). Such a recipe can be readily implemented into N-body codes to study the long-term growth and evolution of planetary systems.\n- 4. In the 2D limit embryos trapped in resonance and on eccentric orbits grow faster than those on circular orbits. Therefore, the secondary planet formation occurs preferentially at resonances with the first giant planet.  \nIn this work we have focused on the 2D limit where all pebbles are in the midplane. However, disk turbulence may lift small pebbles from the midplane. In addition, the inclinations of the embryos/planetesimals can be excited by their mutual gravitational interactions. Under these circumstances, the pebble accretion efficiency is not ε 2D , but rather involves effects related with the planet's inclination, the pebble accretion radius, and the scale height of the pebble layer (Ormel & Kobayashi 2012; Morbidelli et al. 2015; Levison et al. 2015a,b; Xu et al. 2017). These 3D effects will be addressed in the subsequent paper (Paper II).", 'Appendix A: Derivation of the accretion efficiencies': "Appendix A.1: Settling regime - planar limit (2D)  \nAppendix A.1.1: General expression  \nWe use physical, order-of-magnitude arguments to derive the scaling relations of the pebble accretion efficiency (see also Ormel & Klahr (2010); Ormel & Kobayashi (2012); Lambrechts & Johansen (2012); Guillot et al. (2014); Morbidelli et al. (2015); Ida et al. (2016)). These expressions are complemented by prefactors which we obtain from the numerical simulations.  \nThere are two key requirements for pebble accretion (see Ormel (2017) for a review):  \n- 1. During a pebble-planet encounter, the time a pebble settles onto the planet is shorter than its encounter time, t set < t enc such that the pebble is able to hit the planet.\n- 2. The stopping time is shorter than the encounter time, t s < t enc . Otherwise, the gas drag is not important during this encounter.  \nFig. A.1. Sketch illustrating the pebble-planet encounter in a local frame where the pebble is co-moving with the central planet. The approach velocity is given by the relative velocity between these two objects ( ∆ v ). The perturbed and unperturbed trajectories are represented in the blue solid and dashed lines. The important timescales are: the pebble-planet encounter time t enc = b/ ∆ v , and the settling time t set = b/v set , where v set is the sedimentation velocity when the planet's gravity equals to the gas drag.  \n<!-- image -->  \nThe failure of either criterion implies that accretion is not in the settling regime. In Fig. A.1, the settling velocity of a pebble is obtained from balancing the planets' gravitational force with the gas drag force, v set = ( GM p /b 2 ) t s where b is the impact parameter. The settling time is defined as t set = b/v set while the encounter time approximates as t enc glyph[similarequal] b/ ∆ v , where ∆ v is the relative velocity between the planet and the pebble. Therefore, from criterion (1) the accretion radius of the planet in the settling regime, b set can be written as  \nb set ∼ √ GM p t s ∆ v . (A.1)  \nEquivalently, criterion (1) can also be obtained from the gravitational deflection time t g (Perets & Murray-Clay 2011; Lambrechts & Johansen 2012; Baruteau et al. 2016). For settling t g must be shorter than the stopping time, t g < t s . Here t g = ∆ v/ ( GM p /b 2 ) is the time for deflecting the approaching velocity ∆ v by the gravity of the planet.  \nIn the coplanar (2D) accretion case, the mass flux of pebbles accreted by the planet is  \n˙ M PA , 2D ∼ 2 b set ∆ v Σ P ∼ 2 √ GM p t s ∆ v Σ P . (A.2)  \nAccording to the definition, the pebble accretion efficiency for a circular orbit planet is given by  \nε 0 , 2D = ˙ M PA ˙ M disk ∼ √ GM p t s ∆ v (1 + τ 2 s ) 2 πτ s ηv K r . (A.3)  \nRewritten in terms of dimensionless quantities, this expression becomes  \nε 0 , 2D = A 2D √ q p τ s η 2 ( ∆ v v K ) , (A.4)  \nwhere A 2D = 0 . 32 is a fitting factor. In the above expression, the factor (1 + τ 2 s ) has been omitted since we focus on pebbles with τ s < 1 . In addition, dimensionless quantities are used  \nq p ≡ M p M glyph[star] = 3 × 10 -6 ( M p 1 M ⊕ ) ; η ≡ v hw v K = 10 -3 ( v hw 30 m / s ) (A.5)  \nSince the pebble accretion criterion (2) breaks down when the encounter time is shorter than the stopping time, a transition velocity v ∗ can be calculated by t enc = t s  \nv ∗ = ( q p /τ s ) 1 / 3 v K (A.6)  \nWhen ∆ v approaches v ∗ , the accretion transitions from the setting regime to the ballistic regime. We adopt an exponential decay function to fit such a transition 3  \nf set = exp[ -0 . 5(∆ v/v ∗ ) 2 ] . (A.7)  \nWhen ∆ v glyph[lessmuch] v ∗ , accretion operates in the settling regime and f set = 1 ; on the other hand, when ∆ v glyph[greatermuch] v ∗ , accretion is in the ballistic regime and f set = 0 . Therefore, the general expression of the accretion efficiency in the settling regime reads  \nε set , 2D = A 2D √ q p τ s η 2 ( ∆ v v K ) f set , (A.8)", 'Appendix A.1.2: Circular & Eccentric cases': 'For a planet on a circular orbit, its relative velocity ∆ v is in principle the sum of the headwind velocity, v hw glyph[similarequal] ηv K , and the Keplerian shear velocity between the planet and the pebble, v sh ∼ Ω K b set (Fig. 1). In the literature, the pebble accretion in the circular ( e p = 0 ) limit is classified into two regimes depending on which contribution dominates ∆ v : the headwind regime (Bondi regime) or the shear regime (Hill regime) (Lambrechts & Johansen 2012; Guillot et al. 2014). Note that the shear velocity ( v sh ∝ b set ) increases with τ s and M p . Therefore, for given pebble size and disk properties (same τ s and v hw ), planetesimals and low-mass planets tend to be in the headwind regime, whereas massive planets would be in the shear regime. Equivalently, expressed in terms of τ s , accretion of small τ s pebbles takes place in the headwind regime whereas accretion of large τ s pebbles will be in the shear regime (as seen from Fig. 2).  \nWhen v hw glyph[similarequal] v sh (or ε 0 , hw glyph[similarequal] ε 0 , sh ), a transition mass ratio for the above two regimes is defined as q hw / sh = η 3 /τ s .  \nWe construct an expression for the relative velocity ∆ v to be used in Eq. (A.3), which combines both the headwind and the shear regimes. For a circular orbit planet, we define ∆ v = v cir with  \nv cir = [ 1 + a hw / sh ( q p q hw / sh )] -1 v hw + v sh , (A.9)  \nwhere v hw = ηv K and v sh = 0 . 52( q p τ s ) 1 / 3 v K . Numerically, we fit a hw / sh = 5 . 7 . Equation (A.9) provides a smooth transition at the boundary, ensuring that v cir = v hw for q p glyph[lessmuch] q hw / sh and v cir = v sh for q p glyph[greatermuch] q hw / sh .  \nWhen a planet is on a eccentric orbit, the relative velocity in addition includes an eccentricity contribution due to the elliptic  \norbit of the planet. v ecc is the eccentric velocity of the planet relative to its circular Keplerian value. Therefore, in Eq. (A.8)  \n∆ v = max(v cir , v ecc ) . (A.10)  \nWe fit v ecc = 0 . 76 e p v K from simulations, consistent with the analysis of Guillot et al. (2014). When the eccentricity is large that v ecc > v cir , the accretion is in the eccentricity regime instead of the shear/headwind regime.', 'Appendix A.2: Settling regime - 3D limit': 'For completeness, we also include the expression of 3D accretion efficiency here. The mass flux of pebbles accreted by the planet is  \n˙ M set , 3D = b 2 set ∆ vρ P ∼ GM p t s ρ P , (A.11)  \nwhere ρ P = Σ P / ( √ 2 πr p h P ) is the volume density of the pebbles and h P is the aspect ratio of the pebble layer. The expression of the pebble accretion efficiency in the 3D limit reads  \nε set , 3D = A 3D q p ηh P f 2 set (A.12)  \nwhere A 3D will be calculated in Paper II. Note that in the expression for ε set , 3D ∆ v only appears in f set .', 'Appendix A.3: Ballistic regime': 'In the ballistic regime, the accretion radius is significantly reduced due to the lack of gas drag, resulting in a drop of the accretion efficiency (see Fig. 3). The ballistic accretion radius now reads (Safronov 1972)  \nb bal = R p √ ( v esc ∆ v ) 2 +1 . (A.13)  \nThe accretion efficiency in the 2 D ballistic regime becomes  \nε bal , 2D = ∆ vb bal 2 πr p τ s ηv K = R p √ v esc 2 +∆ v 2 2 πr p τ s ηv K . (A.14)  \nWe note that when v ∗ glyph[lessorsimilar] ∆ v glyph[lessorsimilar] v esc , the accretion is in the gravitational focusing regime and ε bal , 2D is independent of the eccentricity. On the other hand, when ∆ v glyph[greaterorsimilar] v esc , the accretion is geometric regime. The efficiency increases with the eccentricity and b bal is reduced to the physical radius of the planet (Guillot et al. 2014).  \nWe rewrite the above formula including the ( 1 -f set ) term  \nε bal , 2D = R p 2 πτ s ηr p √ 2 q p r p R p + ( ∆ v v K ) 2 (1 -f set ) , (A.15)  \nIn the 3D limit, we have  \nε bal , 3D = 1 4 √ 2 πητ s h P ( 2 q p v K ∆ v R p r p + R 2 p r 2 p ∆ v v K ) ( 1 -f 2 set ) (A.16)  \nFor notations of masses, velocities, timescales, etc., see Table A.1.  \nAcknowledgements. Wethank Anders Johansen, Carsten Dominik for useful discussions, and Ramon Brasser for proofreading the manuscript. We also thank the anonymous referee for their insightful suggestions and comments. B.L. and C.W.O are supported by the Netherlands Organization for Scientific Research (NWO; VIDI project 639.042.422).', 'References': 'Alibert, Y. 2017, A&A, 606, A69  \nAndrews, S. M., Rosenfeld, K. A., Kraus, A. L., & Wilner, D. J. 2013, ApJ, 771,  \n129  \nAnsdell, M., Williams, J. P., Manara, C. F., et al. 2017, AJ, 153, 240  \nAnsdell, M., Williams, J. P., van der Marel, N., et al. 2016, ApJ, 828, 46  \nBarenfeld, S. A., Carpenter, J. M., Ricci, L., & Isella, A. 2016, ApJ, 827, 142 Baruteau, C., Bai, X., Mordasini, C., & Mollière, P. 2016, Space Sci. Rev., 205, 77  \nBirnstiel, T., Dullemond, C. P., & Brauer, F. 2010, A&A, 513, A79  \nBirnstiel, T., Klahr, H., & Ercolano, B. 2012, A&A, 539, A148  \n- Bitsch, B., Crida, A., Morbidelli, A., Kley, W., & Dobbs-Dixon, I. 2013, A&A, 549, A124  \n- Bitsch, B. & Kley, W. 2011, A&A, 530, A41\n- Bitsch, B., Lambrechts, M., & Johansen, A. 2015, A&A, 582, A112  \nBitsch, B., Morbidelli, A., Johansen, A., et al. 2018, ArXiv e-prints  \nBrauer, F., Dullemond, C. P., & Henning, T. 2008, A&A, 480, 859  \n- Brouwers, M. G., Vazan, A., & Ormel, C. W. 2017, ArXiv e-prints  \nColeman, G. A. L. & Nelson, R. P. 2014, MNRAS, 445, 479  \n- Cossou, C., Raymond, S. N., Hersant, F., & Pierens, A. 2014, A&A, 569, A56\n- Cresswell, P. & Nelson, R. P. 2008, A&A, 482, 677\n- Cumming, A., Butler, R. P., Marcy, G. W., et al. 2008, PASP, 120, 531  \nDominik, C., Blum, J., Cuzzi, J. N., & Wurm, G. 2007, in Protostars and Planets  \nV, ed. B. Reipurth, D. Jewitt, & K. Keil (Univ. of Arizona Press, Tucson),  \n783-800  \nDraine, B. T. 2006, ApJ, 636, 1114  \nFehlberg, E. 1969, NASA Technical Report  \nGuillot, T., Ida, S., & Ormel, C. W. 2014, A&A, 572, A72  \n- Güttler, C., Blum, J., Zsom, A., Ormel, C. W., & Dullemond, C. P. 2010, A&A,', '513, A56': 'Hood, L. L. 1998, Meteoritics and Planetary Science, 33  \nIda, S., Guillot, T., & Morbidelli, A. 2016, A&A, 591, A72  \nIda, S. & Lin, D. N. C. 2004, ApJ, 604, 388  \nJohansen, A., Blum, J., Tanaka, H., et al. 2014, Protostars and Planets VI, 547  \nJohansen, A. & Klahr, H. 2005, ApJ, 634, 1353  \nJohansen, A. & Lambrechts, M. 2017, Annual Review of Earth and Planetary  \nSciences, 45, 359  \n- Johansen, A., Mac Low, M.-M., Lacerda, P., & Bizzarro, M. 2015, Science Advances, 1, 1500109  \nKrijt, S., Ormel, C. W., Dominik, C., & Tielens, A. G. G. M. 2016, A&A, 586,  \nA20  \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  \nLee, M. H. & Peale, S. J. 2002, ApJ, 567, 596  \nLevison, H. F., Kretke, K. A., & Duncan, M. J. 2015a, Nature, 524, 322  \nLevison, H. F., Kretke, K. A., Walsh, K. J., & Bottke, W. F. 2015b, Proceedings  \nof the National Academy of Science, 112, 14180  \nLin, D. N. C. & Papaloizou, J. C. B. 1993, in Protostars and Planets III, ed. E. H. Levy & J. I. Lunine, 749-835  \nLissauer, J. J. 1987, Icarus, 69, 249  \nLiu, B., Zhang, X., Lin, D. N. C., & Aarseth, S. J. 2015, ApJ, 798, 62  \nMatsumura, S., Brasser, R., & Ida, S. 2017, A&A, 607, A67  \nMayor, M., Marmier, M., Lovis, C., et al. 2011, ArXiv e-prints: 1109.2497  \nMiura, H. & Nakamoto, T. 2005, Icarus, 175, 289  \n- Morbidelli, A., Lambrechts, M., Jacobson, S., & Bitsch, B. 2015, Icarus, 258, 418\n- Morris, M. A., Boley, A. C., Desch, S. J., & Athanassiadou, T. 2012, ApJ, 752, 27\n- Movshovitz, N., Bodenheimer, P., Podolak, M., & Lissauer, J. J. 2010, Icarus, 209, 616\n- Murray, C. D. & Dermott, S. F. 1999, Solar system dynamics  \nNatta, A. & Testi, L. 2004, in Astronomical Society of the Pacific Conference Series, Vol. 323, Star Formation in the Interstellar Medium: In Honor of David Hollenbach, ed. D. Johnstone, F. C. Adams, D. N. C. Lin, D. A. Neufeeld, & E. C. Ostriker, 279  \nOrmel, C. W. 2017, in Astrophysics and Space Science Library, Vol. 445, Astro-  \nphysics and Space Science Library, ed. M. Pessah & O. Gressel, 197  \nOrmel, C. W. & Klahr, H. H. 2010, A&A, 520, A43  \nOrmel, C. W. & Kobayashi, H. 2012, ApJ, 747, 115  \n- Ormel, C. W. & Liu, B. 2018, ArXiv e-prints  \nOrmel, C. W., Shi, J.-M., & Kuiper, R. 2015, MNRAS, 447, 3512  \nPapaloizou, J. C. B. & Larwood, J. D. 2000, MNRAS, 315, 823  \nPascucci, I., Testi, L., Herczeg, G. J., et al. 2016, ApJ, 831, 125  \nPerets, H. B. & Murray-Clay, R. A. 2011, ApJ, 733, 56  \nPérez, L. M., Chandler, C. J., Isella, A., et al. 2015, ApJ, 813, 41  \nPicogna, G., Stoll, M. H. R., & Kley, W. 2018, ArXiv e-prints  \nPollack, J. B., Hubickyj, O., Bodenheimer, P., et al. 1996, Icarus, 124, 62  \nRaymond, S. N. & Izidoro, A. 2017, Icarus, 297, 134  \nRaymond, S. N., Quinn, T., & Lunine, J. I. 2006, Icarus, 183, 265  \nRicci, L., Testi, L., Natta, A., et al. 2010, A&A, 512, A15  \nSafronov, V. S. 1972, Evolution of the protoplanetary cloud and formation of the  \nearth and planets.  \n- Sato, T., Okuzumi, S., & Ida, S. 2016, A&A, 589, A15  \nTazzari, M., Testi, L., Ercolano, B., et al. 2016, A&A, 588, A53  \nTeyssandier, J. & Terquem, C. 2014, MNRAS, 443, 568  \nVenturini, J., Alibert, Y., & Benz, W. 2016, A&A, 596, A90  \n- Visser, R. G. & Ormel, C. W. 2016, A&A, 586, A66  \nWeidenschilling, S. J. 1977, MNRAS, 180, 57  \nWeidenschilling, S. J. & Cuzzi, J. N. 1993, in Protostars and Planets III, ed. E. H.  \nLevy & J. I. Lunine, 1031-1060  \nWilliams, J. P. & Cieza, L. A. 2011, ARA&A, 49, 67  \nWright, J. T., Veras, D., Ford, E. B., et al. 2011, ApJ, 730, 93  \nXu, Z., Bai, X.-N., & Murray-Clay, R. A. 2017, ApJ, 847, 52  \nZheng, X., Lin, D. N. C., & Kouwenhoven, M. B. N. 2017, ApJ, 836, 207  \nZhou, J.-L. & Lin, D. N. C. 2007, ApJ, 666, 447  \n- Zhu, Z., Nelson, R. P., Dong, R., Espaillat, C., & Hartmann, L. 2012, ApJ, 755,'}