Issue 
A&A
Volume 615, July 2018



Article Number  A138  
Number of page(s)  15  
Section  Planets and planetary systems  
DOI  https://doi.org/10.1051/00046361/201732307  
Published online  27 July 2018 
Catching drifting pebbles
I. Enhanced pebble accretion efficiencies for eccentric planets
Anton Pannekoek Institute (API), University of Amsterdam,
Science Park 904,
1090GE
Amsterdam, The Netherlands
email: b.liu@uva.nl,c.w.ormel@uva.nl
Received:
16
November
2017
Accepted:
27
March
2018
Context. Coagulation theory predicts that micronsized 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.
Aims. 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.
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.
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 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 alreadyformed giant planet, we find that the embryo grows quickly due to its higher eccentricity.
Conclusions. 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 Nbody codes to simulate the growth and evolution of planetary systems.
Key words: methods: numerical / planets and satellites: formation
© ESO 2018
1 Introduction
In protoplanetary disks, micronsized 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 micronsized 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} ≪ 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 insideout manner. The interior, planetforming regions are therefore constantly supplied by pebbles that drift from the outer regions of disks (Birnstiel et al. 2012; Lambrechts & Johansen 2014; Sato et al. 2016; Krijt et al. 2016).
Observationally, the existence of pebblesized particles in disks is inferred from optically thin emission at submm 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 pebblesized 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 inlarger bodies, invisible to millimeter disk observations. Indeed, Ansdell et al. (2017) find lower average dust masses for the older σ Orionis cluster.
When 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 gravitational attraction and gas drag (Ormel & Klahr 2010; Lambrechts & Johansen 2012). This process is known as pebble accretion (see Ormel 2017 and Johansen & Lambrechts 2017 for recent reviews). Depending on the importance of gas drag during the pebbleplanet encounter, the accretion can be classified into two regimes: settling and ballistic. In the classical, planetesimaldriven 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.
The 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} (1)
where Ṁ_{PA} is the pebble accretion rate on the planet and Ṁ_{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.
The 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), 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 be implemented in Nbody codes, in order to study the formation and longterm dynamical evolution of planetary systems.
In 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 planetplanet 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 (comoving) 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.
The 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 alreadyformed 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 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 noninertial, 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).
In a 2D limit, the pebble accretion efficiency in the local frame is given by (Guillot et al. 2014; Ormel 2017): (2)
where b is the pebble’s impact distance measured from the xaxis and v_{y}(b) is the pebble’s velocity perpendicular to b (Fig. 1a), Ω_{p} is the angular velocity of the planet, and is the Heaviside step function; means the pebble with a impact distance b hits the planet, otherwise, . The total pebble mass flux is given by Ṁ_{disk} = 2πr_{p}v_{r}Σ_{P}, where v_{r} is the radial drift velocity 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.
Fig. 1 Sketch illustrating pebble accretion in two different frames. a) Local frame: the comoving frame is centred on the planet where the xaxis is pointing radially away from the star and the yaxis 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 xaxis impact distance b and the yaxis velocity v_{y}(b). b) Global frame: the planet orbits around the central star with a semimajor 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. 
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.
The equation of motion for the pebble in heliocentric coordinates is (3)
where 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_{⋆} 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 (v_{K}) and the gas azimuthal velocity, is the headwind pre factor, H_{gas} is the disk scale height and P is the gas pressure at the planet’s location r_{p}.
In 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.
The 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 (4)
Thus, we can obtain the orbital motions of the planet and the pebble from Eqs. (3) and (4).
Integrating 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 vanishes and the pebble’s velocity can be approximated as (5)
We 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.
The efficiency in the global frame is simply (6)
where 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.
Fig. 2 Pebble accretion efficiency ε_{2D} vs τ_{s} and M_{p}. Three planetmasses (M_{p} = 10^{−3} M_{⊕}, 10^{−1} M_{⊕} and 10 M_{⊕}) and nine dimensionless stopping times (τ_{s} ranging from10^{−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} ≲ 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_{⊕}). 
2.3 Three numerical methods
We employ three methods to calculate the pebble accretion efficiency:

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);

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));

Global hybrid – the SCA (Eq. (5)) is applied for r_{p} − r > 2R_{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 and a_{p} is the planet’s semimajor axis.
We 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 scale height of the pebble disk. In this work we only consider the 2D 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.
In the local simulation, the pebbles’ initial x locations (x_{0}) range from − 5R_{H} to 5R_{H}, with an interval of 10^{−4}R_{H}, and y_{0} is initialized at either − 40R_{H} or 40R_{H}. The initial radial and azimuthal velocity of the pebble are and (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.
In 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) + 5R_{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 grid points in both the vertical and radial directions and N_{tot} in total. The initial radial and azimuthal velocity of the pebble are and . 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.
For 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.
3 Results
Results 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_{⊙}. 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 (), respectively. Hereafter the efficiency refers to its mean value from the simulations. We focus on the pebble accretion regime (settling regime) when 10^{−3} ≲ τ_{s} < 1. But for the additional method comparison, the cases of τ_{s} ≳ 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 pebbleplanet encounter (referred to as the threebody regime in Ormel & Klahr 2010).
The strong coupling approximation (SCA) works well for the small τ_{s} pebble where the gaspebble 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} ≲ 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}.
We alsosee 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.
There are two reasons for the difference in ε_{2D} between thetwo methods in the high mass planet case. First, there is a risk of multiple counting accreting pebbles for the localmethod. 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 closeencounter where and Δr_{co} is the widthof the coorbital 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.
The other reason is due to the linearization of the equations of motions. In the local method, we ignore the highorder terms of (n ≥ 2) in the pebble’s equation of motion by virtue of the approximation that Δr_{p} ≪ 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} ≪ r_{p} breaks down and the firstorder 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.
We 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 (≃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.
In 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} ≲ τ_{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.
In 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} ≳ 1). The global method is arguably the better approach to use, in particular for configurations such as resonance trapping, and, 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_{⋆}), the headwind velocity (v_{hw}) and the dimensionless stopping time (τ_{s}) with respectto the planet eccentricity.
3.2.1 Default run
The parameters of the default model are: M_{p} = 0.1 M_{⊕}, M_{⋆} = 1 M_{⊙}, 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.
We find that the pebble accretion efficiency varies with eccentricity. When the planet is on a nearly circular orbit (e_{p} ≲ 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}≃ 0.18 when the eccentricity approaches 0.03. After that, the efficiency quickly reduces to 0.07 when e_{p} ≃ 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.
By balancing the settling time and the encounter time, it can be shown that the pebble accretion radius is (Ormel & Klahr 2010; Lambrechts & Johansen 2012) (7)
where Δv is the relative velocity between the planet and the pebble (see Sect. A).
When 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.
The rapid drop of ε_{2D} at higher e_{p} ≃ 0.04, however, indicates the settling interactions fail. The reason is that the time for the planetpebble encounter decreases as Δv (eccentricity) increases. When the planetpebble 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.
In orderto 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_{⋆} = 1 M_{⊙}, τ_{s} = 10^{−2}, v_{hw} = 60 m s^{−1}. 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 yaxis) and the number of accreted pebbles (right yaxis) 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 and , 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 and . Therefore, the relative velocity between the planet and the pebble can be analytically calculated, which is shown by the blue lines of Fig. 4.
We 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 subKeplerian 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 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.
On 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)).
Fig. 3 Pebble accretion efficiency vs planet eccentricity for the default run: M_{p} = 0.1 M_{⊕}, M_{⋆} = 1 M_{⊙}, 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} ≃ 0. 
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_{⋆}), 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_{⊙}, 1 M_{⊙}, 0.3 M_{⊙}), 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. Figure 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.
The 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 lessmassive 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 planets 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. Figure 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_{⋆} = 1 M_{⊙}) is shown in black, and the massive giant star case (M_{⋆} = 3 M_{⊙}) and the lowmass M dwarf star case (M_{⋆} = 0.3 M_{⊙}) are given in light purple and dark purple, respectively.
For 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 fora 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. Figure 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.
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 maximum 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.
For τ_{s} = 10^{−3} the efficiency of pebble accretion is 0.13 when the planet is on a circular orbit. When e_{p} ≃ 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.
In addition, we find that the ratio of the maximum ε_{2D} to ε_{2D} when e_{p} ≃ 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}.
To summarize, the pebble accretion efficiency ε_{2D} is an increasing function of the planettostar mass ratio (q_{p} = M_{p}∕M_{⋆}), 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} ≲ 10^{−3}), gradually increases with eccentricity when it is moderate (e_{p} ≃ 10^{−2}), and drops quickly when the eccentricity becomes relatively high (e_{p} ≳ 0.1).
Fig. 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 aphelionis at 180°. The planet eccentricity is e_{p} = 0.007 in the left panel and e_{p} = 0.03 in the right panel. Other parameters are the same for the two panels: M_{p} = 0.01 M_{⊕}, M_{⋆} = 1 M_{⊙}, τ_{s} = 10^{−2}, v_{hw} = 60 m s^{−1} and N_{tot} = 300 × 300. The black dashed line represents the transition velocity from the settling to the ballistic accretion, v = v_{*} (Eq. (13)). 
Model setup in a parameter study (planet density is 3 gcm^{−3}).
Fig. 5 Pebble accretion efficiency vs planet eccentricity. The default run is shown in black, with M_{p} = 0.1 M_{⊕}, M_{⋆} = 1 M_{⊙}, v_{hw} = 30 m s^{−1} and τ_{s} = 10^{−2}. In each panel one parameter varies and the other three keep the same. Upper left panel: three different planet masses, M_{p} = 1 M_{⊕}, M_{p} = 0.1 M_{⊕} and M_{p} = 0.01 M_{⊕}; upper right panel: three different stellar masses, M_{⋆} = 3 M_{⊙}, M_{⋆} = 1 M_{⊙} and M_{⋆} = 0.3 M_{⊙}; lower left panel: three different headwind velocities, v_{hw} = 15 m s^{−1}, v_{hw} = 30 m s^{−1} and v_{hw} = 60 m s^{−1}; lower right panel: 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_{⋆}, v_{hw} and τ_{s}. 
3.3 Analytical fit expression
In this subsection, we present analytical fitting formulas, the detailed derivations of these formulas are described at length in Appendix A.
We find that the obtained accretion efficiency in the settling regime can be well fitted by (8)
where q_{p} is the mass ratio of the planet to the central star. The relative velocity between the planet and the pebble (Δv) is (9)
Here v_{cir} is the relative velocity between the circular orbit planet and the pebble. We adopt an expression of v_{cir} that combines the headwind and the shear regimes: (10)
where v_{hw} = ηv_{K}, and the transition mass ratio for tworegimes q_{hw/sh} = η^{3}∕τ_{s}. The eccentric velocity of the planet relative to its circular Keplerian value is v_{ecc}, where we fit (11)
In the epicycle approximation, it can be shown that the velo city 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).
When Δv becomes higher, the accretion is no longer in the settling regime. We adopt an exponential function f_{set} to express the decay of ε_{set,2D}, (12)
is the transition velocity from the settling to the ballistic regimes (see derivations in Sect. A.1).
The accretion efficiency in the ballistic regime is (14)
It 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 twoquantities, R_{p}∕r_{p}.
The total accretion efficiency is (15)
The 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} ≤ 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 . Clearly, in such a expression, when ε_{2D} ≪ 1 and when ε_{2D} ≥ 1. Accounting for these probabilistic nature of ε, we give the efficiency expression as (16)
to ensure that ε ≤ 1. In Figs. 2 and 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.
Lambrechts & Johansen (2014) also calculate the efficiency when the planet is on a circular orbit. In the sheardominated regime, Eq. (8) can be written as (17)
The latter expression of the above formula adopts the same disk model as Lambrechts & Johansen (2014). Comparing Eq. (17) with Lambrechts & Johansen (2014)’s Eq. (33), we find that the scaling relations are identical, and our prefactor is 35% lower than theirs.
In the headwinddominated regime, Eq. (8) can be written as (18)
The transition mass from the headwind to the shear regimes is (19)
In the eccentricitydominated regime, Eq. (8) can be written as (neglecting f_{set}): (20)
We find that in this regime ε_{set,2D} increases as .
3.4 Neglected effects
In deriving the expression above, we have neglected several effects: 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 lowmass 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).
2. Aerodynamic deflection (Guillot et al. 2014; Johansen et al. 2015; Visser & Ormel 2016), relevant for very lowmass 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} ≲ 100 km and τ_{s} ≲ 10^{−3}.
3. Preplanetary 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 preplanetary 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 ( 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, , where c_{s} is the sound speed at the planet’s location. As long as R_{Bondi} < b_{set} (or 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_{⊕}.
In addition, we assume that the gas moves on an unperturbed subKeplerian 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}.
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 (Lin & Papaloizou 1993), M_{p, iso} is around 20 M_{⊕} for a solar mass star (; Lambrechts et al. 2014; Bitsch et al. 2018). In our study, we consider the planet mass in a range of 10^{−3} M_{⊕}≲ M_{p} ≲ 10 M_{⊕} < M_{p,iso}.
4 Application: formation of a secondary planet
In this section we envision an alreadyformed 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.
Based on the observational statistics of the occurrence rate, 10%–20% of solartype 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 during the gasrich 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.
Once 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 Nbody simulations to simulate the embryo’s orbital and mass evolution. We explore how massive the pebble disk must be to produce a secondary planet.
We restrict our problem to a quiescent (low turbulent) disk such that pebbles are settled into the disk midplane. The accretion is thus in the 2D regime (pebble scale height H_{P} is smaller than impact radius b_{set}). The influence of disk turbulence and the prescription of 3D pebble 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) (21)
where r_{AU} = r∕(1 AU), the headwind prefactor and the dusttogas ratio is 1%.
The formation of the first gas giant may already take a few Myr, depending on the opacity in its envelope (Movshovitz et al. 2010). Thus, Ṁ_{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 consistentwith advanced dust coagulation model (Birnstiel et al. 2012) as well as disk observations (Pérez et al. 2015; Tazzari et al. 2016).
A gas giant of 1M_{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.
We consider two scenarios: (i) a resonant case where the embryo experiences type I migration and is captured in resonance; and (ii) a nonresonant 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 Nbody code based on Papaloizou & Larwood (2000) and Cresswell & Nelson (2008): a_{m} = −v_{p}∕τ_{m} and . The embryo’s migration timescale τ_{m} and the eccentricity damping timescale τ_{e} are adopted from Eqs. (11) and (13) in Cresswell & Nelson (2008). These expressions take into account the supersonic regime when e_{p} ≳ 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.
The 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}.
Figure 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 xaxis) 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 equilibrium 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 . In Fig. 6a we also find e_{eq} ≃ 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 nonzero due to weak, secular perturbations of the inner gas giant.
We find in Fig. 6b that after 3 × 10^{4} yr the pebble accretion efficiency in the resonant case is higher than the nonresonant, 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 nonresonant planet. And this mass difference increases with time.
In 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_{⊕} superEarth, whereas the embryo in a nonresonant 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 nonresonant 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 nonresonant case, the slow growth of the embryo results in a final superEarth mass planet.
To 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.
Fig. 6 Eccentricity (panel a), pebble accretion efficiency (panel b) and mass of the embryo (panel c) vs total mass of drifting pebbles (lower xaxis) and time (upper xaxis). The black line corresponds to the resonant case with type I migration and resonance trapping, while the red line represents the nonresonant, comparison case that the embryo without type I migration moves on a nearly circular orbit. 
5 Conclusions
Pebble accretion is an important mechanism that drives planet growth by accreting millimeter to centimeter size particles in gasrich 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 2D limit by conducting a series of numerical integrations of the pebble’s equation of motion in both the local (comoving) frame and the global (heliocentric) frame.
The key findings of this study are the following:

The results of the local and global methods are generally consistent. However, the global method more accurately simulates the pebbleplanet interaction, due to the fact that it accounts for curvature effects and models the dynamics of the 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).

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 a factor of 3−5 higher than planets on circular orbits (Sect. 3.2).

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 planettostar 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 Nbody codes to study the longterm growth and evolution of planetary systems.

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.
In 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).
Acknowledgements
We thank Anders Johansen and 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).
Appendix A Derivation of the accretion efficiencies
A.1 Settling regime – planar limit (2D)
A.1.1 General expression
We use physical, orderofmagnitude 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.
Fig. A.1 Sketch illustrating the pebbleplanet encounter in a local frame where the pebble is comoving 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 pebbleplanet 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. 
There are two key requirements for pebble accretion (see Ormel 2017 for a review):

During a pebbleplanet 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.

The stopping time is shorter than the encounter time, t_{s} < t_{enc}. Otherwise, the gas drag is not important during this encounter.
The 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} ≃ 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 (A.1)
Equivalently, criterion (1) can also be obtained from the gravitational deflection time t_{g} (Perets & MurrayClay 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.
In the coplanar (2D) accretion case, the mass flux of pebbles accreted by the planet is (A.2)
According to the definition, the pebble accretion efficiency for a circular orbit planet is given by (A.3)
Rewritten in terms of dimensionless quantities, this expression becomes (A.4)
where A_{2D} = 0.32 is a fitting factor. In the above expression, the factor has been omitted since we focus on pebbles with τ_{s} < 1. In addition, dimensionless quantities are used (A.5)
Since 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}: (A.6)
When Δ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}: (A.7)
When Δv ≪ v_{*}, accretion operates in the settling regime and f_{set} = 1; on the other hand, when Δv ≫ v_{*}, accretion is in the ballistic regime and f_{set} = 0. Therefore, the general expression of the accretion efficiency in the settling regime reads (A.8)
A.1.2 Circular and eccentric cases
For a planet on a circular orbit, its relative velocity Δv is in principle the sum of the headwind velocity, v_{hw} ≃ ηv_{K}, and the Keplerian shear velocity between the planet and the pebble, v_{sh} ~ Ωkb_{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 diskproperties (same τ_{s} and v_{hw}), planetesimals and lowmass 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).
When v_{hw} ≃ v_{sh} (or ε_{0,hw} ≃ ε_{0,sh}), a transition mass ratio for the above two regimes is defined as q_{hw/sh} = η^{3}∕τ_{s}.
We 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 (A.9)
where v_{hw} = ηv_{K} and . 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} ≪ q_{hw/sh} and v_{cir} = v_{sh} for q_{p} ≫ q_{hw/sh}.
When a planet is on a eccentric orbit, the relative velocity in addition includes an eccentricity contribution due to the elliptic orbit of the planet. v_{ecc} is the eccentric velocity of the planet relative to its circular Keplerian value. Therefore, in Eq. (A.8) (A.10)
We fit v_{ecc} = 0.76e_{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.
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 (A.11)
where 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 (A.12)
where A_{3D} will be calculated in Paper II. Note thatin the expression for ε_{set,3D} Δv only appears in f_{set}.
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) (A.13)
The accretion efficiency in the 2D ballistic regime becomes (A.14)
We note that when v_{*}≲ Δv ≲ v_{esc}, the accretion is in the gravitational focusing regime and ε_{bal,2D} is independent of the eccentricity. On the other hand, when Δv ≳ 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).
We rewrite the above formula including the (1 − f_{set}) term (A.15)
In the 3D limit, we have (A.16)
For notations of masses, velocities, timescales, etc., see Table A.1.
List of notations.
References
 Alibert, Y. 2017, A&A, 606, A69 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Andrews, S. M., Rosenfeld, K. A., Kraus, A. L., & Wilner, D. J. 2013, ApJ, 771, 129 [NASA ADS] [CrossRef] [Google Scholar]
 Ansdell, M., Williams, J. P., van der Marel, N., et al. 2016, ApJ, 828, 46 [Google Scholar]
 Ansdell, M., Williams, J. P., Manara, C. F., et al. 2017, AJ, 153, 240 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Barenfeld, S. A., Carpenter, J. M., Ricci, L., & Isella, A. 2016, ApJ, 827, 142 [NASA ADS] [CrossRef] [Google Scholar]
 Baruteau, C., Bai, X., Mordasini, C., & Mollière, P. 2016, Space Sci. Rev., 205, 77 [NASA ADS] [CrossRef] [Google Scholar]
 Birnstiel, T., Dullemond, C. P., & Brauer, F. 2010, A&A, 513, A79 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Birnstiel, T., Klahr, H., & Ercolano, B. 2012, A&A, 539, A148 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bitsch, B., & Kley, W. 2011, A&A, 530, A41 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bitsch, B., Crida, A., Morbidelli, A., Kley, W., & DobbsDixon, I. 2013, A&A, 549, A124 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bitsch, B., Lambrechts, M., & Johansen, A. 2015, A&A, 582, A112 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bitsch, B., Morbidelli, A., Johansen, A., et al. 2018, A&A, 612, A30 [Google Scholar]
 Brauer, F., Dullemond, C. P., & Henning, T. 2008, A&A, 480, 859 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Brouwers, M. G., Vazan, A., & Ormel, C. W. 2018, A&A, 611, A65 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Coleman, G. A. L., & Nelson, R. P. 2014, MNRAS, 445, 479 [NASA ADS] [CrossRef] [Google Scholar]
 Cossou, C., Raymond, S. N., Hersant, F., & Pierens, A. 2014, A&A, 569, A56 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Cresswell, P., & Nelson, R. P. 2008, A&A, 482, 677 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Cumming, A., Butler, R. P., Marcy, G. W., et al. 2008, PASP, 120, 531 [CrossRef] [Google Scholar]
 Dominik, C., Blum, J., Cuzzi, J. N., & Wurm, G. 2007, in Protostars and Planets V, eds. B. Reipurth, D. Jewitt, & K. Keil (Tucson: University of Arizona Press), 783 [Google Scholar]
 Draine, B. T. 2006, ApJ, 636, 1114 [NASA ADS] [CrossRef] [Google Scholar]
 Fehlberg, E. 1969, NASA Technical Report [Google Scholar]
 Guillot, T., Ida, S., & Ormel, C. W. 2014, A&A, 572, A72 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Güttler, C., Blum, J., Zsom, A., Ormel, C. W., & Dullemond, C. P. 2010, A&A, 513, A56 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hood, L. L. 1998, Meteorit. Planet. Sci., 33, 97 [Google Scholar]
 Ida, S., & Lin, D. N. C. 2004, ApJ, 604, 388 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Ida, S., Guillot, T., & Morbidelli, A. 2016, A&A, 591, A72 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Johansen, A., & Klahr, H. 2005, ApJ, 634, 1353 [NASA ADS] [CrossRef] [Google Scholar]
 Johansen, A., & Lambrechts, M. 2017, Annu. Rev. Earth Planet. Sci., 45, 359 [Google Scholar]
 Johansen, A., Blum, J., Tanaka, H., et al. 2014, Protostars and Planets VI, 547 [Google Scholar]
 Johansen, A., Mac Low, M.M., Lacerda, P., & Bizzarro, M. 2015, Sci. Adv., 1, 1500109 [NASA ADS] [CrossRef] [Google Scholar]
 Krijt, S., Ormel, C. W., Dominik, C., & Tielens, A. G. G. M. 2016, A&A, 586, A20 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lambrechts, M., & Johansen, A. 2012, A&A, 544, A32 [Google Scholar]
 Lambrechts, M., & Johansen, A. 2014, A&A, 572, A107 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lambrechts, M., Johansen, A., & Morbidelli, A. 2014, A&A, 572, A35 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lee, M. H., & Peale, S. J. 2002, ApJ, 567, 596 [NASA ADS] [CrossRef] [Google Scholar]
 Levison, H. F., Kretke, K. A., & Duncan, M. J. 2015a, Nature, 524, 322 [CrossRef] [Google Scholar]
 Levison, H. F., Kretke, K. A., Walsh, K. J., & Bottke, W. F. 2015b, Proc. Nat. Acad. Sci., 112, 14180 [Google Scholar]
 Lin, D. N. C., & Papaloizou, J. C. B. 1993, in Protostars and Planets III, eds. E. H. Levy, & J. I. Lunine, 749 [Google Scholar]
 Lissauer, J. J. 1987, Icarus, 69, 249 [NASA ADS] [CrossRef] [Google Scholar]
 Liu, B., Zhang, X., Lin, D. N. C., & Aarseth, S. J. 2015, ApJ, 798, 62 [NASA ADS] [CrossRef] [Google Scholar]
 Matsumura, S., Brasser, R., & Ida, S. 2017, A&A, 607, A67 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mayor, M., Marmier, M., Lovis, C., et al. 2011, ArXiv eprints [arXiv:1109.2497] [Google Scholar]
 Miura, H., & Nakamoto, T. 2005, Icarus, 175, 289 [NASA ADS] [CrossRef] [Google Scholar]
 Morbidelli, A., Lambrechts, M., Jacobson, S., & Bitsch, B. 2015, Icarus, 258, 418 [NASA ADS] [CrossRef] [Google Scholar]
 Morris, M. A., Boley, A. C., Desch, S. J., & Athanassiadou, T. 2012, ApJ, 752, 27 [NASA ADS] [CrossRef] [Google Scholar]
 Movshovitz, N., Bodenheimer, P., Podolak, M., & Lissauer, J. J. 2010, Icarus, 209, 616 [NASA ADS] [CrossRef] [Google Scholar]
 Murray, C. D., & Dermott, S. F. 1999, Solar System Dynamics [Google Scholar]
 Natta, A., & Testi, L. 2004, Star Formation in the Interstellar Medium: In Honor of David Hollenbach, eds. D. Johnstone, F. C. Adams, D. N. C. Lin, D. A. Neufeeld, & E. C. Ostriker, ASP Conf. Ser., 323, 279 [NASA ADS] [Google Scholar]
 Ormel, C. W. 2017, in Formation, Evolution, and Dynamics of Young Solar Systems, eds. M. Pessah & O. Gressel, Astrophys. Space Sci. Lib., 445, 197 [NASA ADS] [CrossRef] [Google Scholar]
 Ormel, C. W., & Klahr, H. H. 2010, A&A, 520, A43 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ormel, C. W., & Kobayashi, H. 2012, ApJ, 747, 115 [NASA ADS] [CrossRef] [Google Scholar]
 Ormel, C. W., Shi, J.M., & Kuiper, R. 2015, MNRAS, 447, 3512 [NASA ADS] [CrossRef] [Google Scholar]
 Papaloizou, J. C. B., & Larwood, J. D. 2000, MNRAS, 315, 823 [NASA ADS] [CrossRef] [Google Scholar]
 Pascucci, I., Testi, L., Herczeg, G. J., et al. 2016, ApJ, 831, 125 [NASA ADS] [CrossRef] [Google Scholar]
 Perets, H. B., & MurrayClay, R. A. 2011, ApJ, 733, 56 [NASA ADS] [CrossRef] [Google Scholar]
 Pérez, L. M., Chandler, C. J., Isella, A., et al. 2015, ApJ, 813, 41 [NASA ADS] [CrossRef] [Google Scholar]
 Picogna, G., Stoll, M. H. R., & Kley, W. 2018, A&A, in press, DOI: 10.1051/00046361/201732523 [Google Scholar]
 Pollack, J. B., Hubickyj, O., Bodenheimer, P., et al. 1996, Icarus, 124, 62 [NASA ADS] [CrossRef] [Google Scholar]
 Raymond, S. N., & Izidoro, A. 2017, Icarus, 297, 134 [NASA ADS] [CrossRef] [Google Scholar]
 Raymond, S. N., Quinn, T., & Lunine, J. I. 2006, Icarus, 183, 265 [NASA ADS] [CrossRef] [Google Scholar]
 Ricci, L., Testi, L., Natta, A., et al. 2010, A&A, 512, A15 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Safronov, V. S. 1972, Evolution of the Protoplanetary Cloud and Formation of the Earth and Planets [Google Scholar]
 Sato, T., Okuzumi, S., & Ida, S. 2016, A&A, 589, A15 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Tazzari, M., Testi, L., Ercolano, B., et al. 2016, A&A, 588, A53 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Teyssandier, J., & Terquem, C. 2014, MNRAS, 443, 568 [NASA ADS] [CrossRef] [Google Scholar]
 Venturini, J., Alibert, Y., & Benz, W. 2016, A&A, 596, A90 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Visser, R. G., & Ormel, C. W. 2016, A&A, 586, A66 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Weidenschilling, S. J. 1977, MNRAS, 180, 57 [NASA ADS] [CrossRef] [Google Scholar]
 Weidenschilling, S. J., & Cuzzi, J. N. 1993, Protostars and Planets III, eds. E. H. Levy, & J. I. Lunine, 1031 [Google Scholar]
 Williams, J. P., & Cieza, L. A. 2011, ARA&A, 49, 67 [NASA ADS] [CrossRef] [Google Scholar]
 Wright, J. T., Veras, D., Ford, E. B., et al. 2011, ApJ, 730, 93 [NASA ADS] [CrossRef] [Google Scholar]
 Xu, Z., Bai, X.N., & MurrayClay, R. A. 2017, ApJ, 847, 52 [NASA ADS] [CrossRef] [Google Scholar]
 Zheng, X., Lin, D. N. C., & Kouwenhoven, M. B. N. 2017, ApJ, 836, 207 [NASA ADS] [CrossRef] [Google Scholar]
 Zhou, J.L., & Lin, D. N. C. 2007, ApJ, 666, 447 [NASA ADS] [CrossRef] [Google Scholar]
 Zhu, Z., Nelson, R. P., Dong, R., Espaillat, C., & Hartmann, L. 2012, ApJ, 755, 6 [NASA ADS] [CrossRef] [Google Scholar]
This definition is the same as the filtering factor/efficiency of Guillot et al. (2014) and Lambrechts & Johansen (2014).
The form of the modulation factor is very similar to Ormel & Klahr (2010), Ormel & Kobayashi (2012) and Visser & Ormel (2016), but the numerical factors are slightly different from these works.
All Tables
All Figures
Fig. 1 Sketch illustrating pebble accretion in two different frames. a) Local frame: the comoving frame is centred on the planet where the xaxis is pointing radially away from the star and the yaxis 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 xaxis impact distance b and the yaxis velocity v_{y}(b). b) Global frame: the planet orbits around the central star with a semimajor 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. 

In the text 
Fig. 2 Pebble accretion efficiency ε_{2D} vs τ_{s} and M_{p}. Three planetmasses (M_{p} = 10^{−3} M_{⊕}, 10^{−1} M_{⊕} and 10 M_{⊕}) and nine dimensionless stopping times (τ_{s} ranging from10^{−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} ≲ 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_{⊕}). 

In the text 
Fig. 3 Pebble accretion efficiency vs planet eccentricity for the default run: M_{p} = 0.1 M_{⊕}, M_{⋆} = 1 M_{⊙}, 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} ≃ 0. 

In the text 
Fig. 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 aphelionis at 180°. The planet eccentricity is e_{p} = 0.007 in the left panel and e_{p} = 0.03 in the right panel. Other parameters are the same for the two panels: M_{p} = 0.01 M_{⊕}, M_{⋆} = 1 M_{⊙}, τ_{s} = 10^{−2}, v_{hw} = 60 m s^{−1} and N_{tot} = 300 × 300. The black dashed line represents the transition velocity from the settling to the ballistic accretion, v = v_{*} (Eq. (13)). 

In the text 
Fig. 5 Pebble accretion efficiency vs planet eccentricity. The default run is shown in black, with M_{p} = 0.1 M_{⊕}, M_{⋆} = 1 M_{⊙}, v_{hw} = 30 m s^{−1} and τ_{s} = 10^{−2}. In each panel one parameter varies and the other three keep the same. Upper left panel: three different planet masses, M_{p} = 1 M_{⊕}, M_{p} = 0.1 M_{⊕} and M_{p} = 0.01 M_{⊕}; upper right panel: three different stellar masses, M_{⋆} = 3 M_{⊙}, M_{⋆} = 1 M_{⊙} and M_{⋆} = 0.3 M_{⊙}; lower left panel: three different headwind velocities, v_{hw} = 15 m s^{−1}, v_{hw} = 30 m s^{−1} and v_{hw} = 60 m s^{−1}; lower right panel: 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_{⋆}, v_{hw} and τ_{s}. 

In the text 
Fig. 6 Eccentricity (panel a), pebble accretion efficiency (panel b) and mass of the embryo (panel c) vs total mass of drifting pebbles (lower xaxis) and time (upper xaxis). The black line corresponds to the resonant case with type I migration and resonance trapping, while the red line represents the nonresonant, comparison case that the embryo without type I migration moves on a nearly circular orbit. 

In the text 
Fig. A.1 Sketch illustrating the pebbleplanet encounter in a local frame where the pebble is comoving 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 pebbleplanet 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. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.