© ESO 2018

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 ≪ 10⁻³ and are strongly coupled to the gas. Particles that have become aerodynamically larger (τ_s ~ 10⁻²−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 et al. 2012; Lambrechts & Johansen 2014; Sato et al. 2016; Krijt et al. 2016).

Observationally, 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 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 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.

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)

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⁻³), the gained planetary mass will be 10⁻² 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 N-body codes, in order to study the formation and long-term 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 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.

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 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 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 pebble-planet 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 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 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_pv_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 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 vy(b). b) Global frame: the planet orbits around the central star with a semi-major axis ap. 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 Mp. Three planetmasses (Mp = 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, Ntot = 3 × 104 for Mp = 10−3 M⊕, Ntot = 3000 for Mp = 10−1 M⊕ and Ntot = 300 for Mp = 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., Mp = 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 semi-major 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₀) range from − 5R_H to 5R_H, with an interval of 10⁻⁴R_H, and y₀ 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₀), 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₀ = 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₀ to r₀ + 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⁻⁸ 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⁻³ 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⁻³ ≲ τ_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 pebble-planet encounter (referred to as the three-body regime in Ormel & Klahr 2010).

The 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 ≲ 0.1. Moreover, the hybrid algorithm significantly reduces the computational time. For a 0.1 M_⊕ planet accreting pebbles of τ_s = 10⁻², the hybrid integration is 40 times faster than the direct integration. It becomes 200 times faster when τ_s is 10⁻³.

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⁻³ 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 close-encounter where and Δr_co is the widthof 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.

The other reason is due to the linearization of the equations of motions. In the local method, we ignore the high-order 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 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.

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 pebble-planet 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⁻³ ≲ τ_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⁻¹ and τ_s = 10⁻². 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⁻⁵ 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⁻³), 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 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.

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⁻², 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⁻³ (Fig. 4a) and e_p = 3 × 10⁻² (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 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 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: Mp = 0.1 M⊕, M⋆ = 1 M⊙, vhw = 30 m s−1, and τs = 10−2. Eccentricities have been explored from ep = 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 ep ≃ 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⁻³ M_⊕, 10⁻² M_⊕, 10⁻¹ M_⊕, 1 M_⊕), three stellar mass (3 M_⊙, 1 M_⊙, 0.3 M_⊙), three headwind velocities (15 m s⁻¹, 30 m s⁻¹, 60 m s⁻¹) and three dimensionless stopping times (10⁻³, 10⁻², 10⁻¹). 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 low-mass 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⁻¹) and high speed (v_hw = 60 m s⁻¹) shown in dark blue and light blue, respectively. Efficiencies obtained from the default run (v_hw = 30 m s⁻¹) 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⁻²) is shown in black while small pebbles (τ_s = 10⁻³) and large pebbles (τ_s = 10⁻¹) are shown in dark green and light green, respectively.

For τ_s = 10⁻³ 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⁻¹ 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⁻³ and low for τ_s = 10⁻¹. 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 planet-to-star 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⁻³), gradually increases with eccentricity when it is moderate (e_p ≃ 10⁻²), 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 ep = 0.007 in the left panel and ep = 0.03 in the right panel. Other parameters are the same for the two panels: Mp = 0.01 M⊕, M⋆ = 1 M⊙, τs = 10−2, vhw = 60 m s−1 and Ntot = 300 × 300. The black dashed line represents the transition velocity from the settling to the ballistic accretion, v = v* (Eq. (13)).

Fig. 5 Pebble accretion efficiency vs planet eccentricity. The default run is shown in black, with Mp = 0.1 M⊕, M⋆ = 1 M⊙, vhw = 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, Mp = 1 M⊕, Mp = 0.1 M⊕ and Mp = 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, vhw = 15 m s−1, vhw = 30 m s−1 and vhw = 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 Mp, and decreases with M⋆, vhw 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 = η³∕τ_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)

and (13)

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 shear-dominated 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 headwind-dominated regime, Eq. (8) can be written as (18)

The transition mass from the headwind to the shear regimes is (19)

In the eccentricity-dominated 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 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).

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 ≲ 100 km and τ_s ≲ 10⁻³.

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 ( 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⁻² 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 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.

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⁻³ M_⊕≲ M_p ≲ 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.

Based 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 during 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.

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

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 dust-to-gas 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⁻⁵ M_⊕∕yr is the pebble mass flux estimated from Lambrechts & Johansen (2014; assuming t = 5 × 10⁶ 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³ 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 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 planet-disk interaction and eccentricity damping into the N-body 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².

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 x-axis) and time (upper x-axis). 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⁻³; red line in Fig. 6a). It is non-zero due to weak, secular perturbations of the inner gas giant.

We find in Fig. 6b that after 3 × 10⁴ 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.

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_⊕ 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.

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

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 2D limit 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.

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 pebble-planet 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⁻²−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 planet-to-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.

• 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, 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.

Fig. 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 tenc = b∕Δv, and the settling time tset = b∕vset, where vset 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 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.

• 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²)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 & 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²) 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³: (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 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).

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 = η³∕τ_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_pv_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.

References

All Tables

All Figures

Fig. 1 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 vy(b). b) Global frame: the planet orbits around the central star with a semi-major axis ap. 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 Mp. Three planetmasses (Mp = 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, Ntot = 3 × 104 for Mp = 10−3 M⊕, Ntot = 3000 for Mp = 10−1 M⊕ and Ntot = 300 for Mp = 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., Mp = 10 M⊕).

In the text

Fig. 3 Pebble accretion efficiency vs planet eccentricity for the default run: Mp = 0.1 M⊕, M⋆ = 1 M⊙, vhw = 30 m s−1, and τs = 10−2. Eccentricities have been explored from ep = 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 ep ≃ 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 ep = 0.007 in the left panel and ep = 0.03 in the right panel. Other parameters are the same for the two panels: Mp = 0.01 M⊕, M⋆ = 1 M⊙, τs = 10−2, vhw = 60 m s−1 and Ntot = 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 Mp = 0.1 M⊕, M⋆ = 1 M⊙, vhw = 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, Mp = 1 M⊕, Mp = 0.1 M⊕ and Mp = 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, vhw = 15 m s−1, vhw = 30 m s−1 and vhw = 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 Mp, and decreases with M⋆, vhw 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 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.

In the text

Fig. 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 tenc = b∕Δv, and the settling time tset = b∕vset, where vset is the sedimentation velocity when the planet’s gravity equals to the gas drag.

In the text