© ESO, 2012

1. Introduction

The gas giants (Jupiter and Saturn) and ice giants (Uranus and Neptune) in our solar system consist of a dense rocky and/or icy core surrounded by a varying degree of hydrogen and helium atmosphere (Guillot 2005). The strong positive correlation between stellar metallicity and exoplanet occurrence (Santos et al. 2004; Fischer & Valenti 2005) is also accompanied by a correlation between stellar metallicity and the amount of heavy elements present in the exoplanetary interior (Guillot et al. 2006; Miller & Fortney 2011), for objects in the gas giant mass range between about 0.3 and 10 Jupiter masses (M_J). Additionally, a careful statistical inspection of the planet candidates from the Kepler transit survey reveals the evaporation and sublimation of the smaller ice and gas giant planets to their naked cores as they get close to their host star (Youdin 2011). After their formation, the migration of these massive planets in the later stages of the protoplanetary disc shapes the final planetary system (Walsh et al. 2011). However, reconstructing how ice and gas giants form in the first place has proven to be challenging.

In the disc instability scenario, gravitational instabilities in the protoplanetary disc excite dense spiral arms which fragment directly into gas giant planets (Boss 1997). The core accretion (or nucleated instability) scenario requires the formation of a 10-Earth-mass (M_⊕) solid core, capable of holding on to a gaseous atmosphere. When the envelope reaches a mass comparable to the core mass, a run-away accretion of the surrounding gas is triggered (Mizuno 1980; Pollack et al. 1996).

Using the solar system as a template for the end result of planet formation is challenged by direct imaging of planetary companions to A-stars at wide stellar separations. The system HR 8799, for example, contains at least 4 planets separated from their host star by 14.5, 24, 38 and 68 AU, confined by an inner and outer debris disc (Marois et al. 2010). Best estimates of the planetary masses place them all in the gas-giant range. The presence of the debris discs reveals that growth to planetesimals occurs at wide orbital distances as well. Another example of a directly imaged gas-giant planet, β Pictoris b (Lagrange et al. 2010), orbits the host star at approximately 10 AU. Fomalhaut b detected in reflected visble light (Kalas et al. 2008), with an upper mass below 1 M_J (Janson et al. 2012), is located far from the central start at approximately 120 AU. LkCa 15b is a newly discovered gas-giant planet of about 6 M_J, likely caught in the epoch of formation, orbiting at approximately 20 AU around a young solar-like star, with an estimated age of only 2 Myr (Kraus & Ireland 2012).

Formation of gas giants by direct gravitational collapse has been shown to be problematic. At large distances from the host star, Kratter et al. (2010) point out that it becomes increasingly difficult to clump gas with masses below the deuterium burning limit. Additionally, at smaller stellar separations, gas cools too slowly for the spiral arms to fragment into bound clumps (Matzner & Levin 2005; Rafikov 2005).

On the other hand, gas-giant formation by core accretion suffers from exceedingly long timescales at wide stellar separations (Dodson-Robinson et al. 2009; Rafikov 2011). Observations of dust infra-red emission (Haisch et al. 2001; Currie et al. 2009) and disc accretion (Jayawardhana et al. 2006) limit the lifetime of the gaseous component of the protoplanetary disc to 10^6...7 yr. Classical core formation by runaway planetesimal accretion is believed to take more than 10⁷ yr beyond 5 AU, where the planetesimal number densities are low (Goldreich et al. 2004). Planetesimals (>km) get gravitationally focused on to the core, but this effect can be significantly reduced when scattering events drive up the random velocity component of the planetesimals.

The formation of planetesimals, larger-than-km-sized solid bodies bound by self-gravity, is problematic in its own way. While classically considered to be the building blocks of both rocky planets and gas-giant cores, the formation of solids this size remains difficult to explain both theoretically and experimentally. Particle growth beyond cm-sizes by coagulation is inefficient (Blum & Wurm 2008; Brauer et al. 2008; Windmark et al. 2012) and radial drift timescales for m-sized boulders are as short as a hundred orbital timescales (Adachi et al. 1976; Weidenschilling 1977). On the other hand, one can circumvent this so-called meter barrier with turbulence induced by the magnetorotational instability (MRI, Balbus & Hawley 1991), which excites local pressure bumps, ideal regions for dust particle trapping and growth (Whipple 1972; Johansen et al. 2009a). In dead zones where the MRI does not operate, streaming instabilities can destabilize the relative motion between gas and particles (Youdin & Goodman 2005; Johansen & Youdin 2007; Bai & Stone 2010) and lead to the formation of dense filaments. When the particle density becomes sufficiently high, large Ceres-sized planetesimals are formed through gravitational collapse (Johansen et al. 2007). The streaming instability benefits strongly from increased disc metallicities (Johansen et al. 2009b; Bai & Stone 2010), explaining partly the higher occurrence rate of exoplanets around higher metallicity stars.

Instead of building up cores of ice and gas giants with planetesimals, we investigate in this paper the accretion of smaller particles, coupled to the gas on approximately orbital timescales. Dust continuum observations of young circumstellar discs around low-mass pre-main-sequence stars show growth of the dominant particle size to mm and cm sizes within less than 1 Myr (Testi et al. 2003; Wilner et al. 2005; Rodmann et al. 2006). The dynamics of these small particles is influenced by the presence of the surrounding protoplanetary gas (Weidenschilling 1977), through Epstein drag (Epstein 1924).

While drag helps reducing the random velocities of large planetesimals (≥1 km), Rafikov (2004) carefully investigated analytically the effect of drag on smaller bodies (≤1 km) as assumed products of a collisional cascade. However, he did not consider particles coupled to the gas on shorter-than-orbital timescales, excluding the pebble-sized objects seen in T Tauri discs. He finds that nearly all fragments settle to the midplane of the nebula and that gas drag is efficient enough to prevent dynamical excitation, making core formation possible within nearly 10⁶ yr, as was later confirmed in coagulation models by Kenyon & Bromley (2009). Accretion of smaller, pebble-sized particles onto protoplanets was first investigated by Johansen & Lacerda (2010), who numerically found that pebbles are accreted from the entire Hill sphere, the region roughly corresponding to the maximal gravitational reach of the core. They identify a prograde particle disc, which could explain the spin periods of asteroids and preferential prograde spin of large asteroids. The influence of gas drag on the interaction of single small bodies and low-mass planets was explored by Ormel & Klahr (2010). Analytically, they calculated that protoplanets starting from  ~10³ km can efficiently accrete  ~cm-sized particles with impact parameters comparable to the radius of the Hill sphere. Ormel & Kobayashi (2012) further investigated the protoplanet growth stage with a thorough toy model including fragments, planetesimals and embryos and stressed the importance of the gas disc properties, such as a reduced local headwind and turbulence for fast growth. Perets & Murray-Clay (2011) analytically investigated the coalescence of binary planetesimals due to drag forces and commented on the possibility of growth through this mechanism. Lyra et al. (2008a) had already earlier ran full disc models of pressure bumps formed near the edges of the dead zone. After merely 200 orbits, they observed bound embryos with masses similar to the planet Mars, consisting of pebble-sized particles.

In this paper, we investigate core growth from a seed mass by gas-drag-aided capture of cm-sized pebbles. In Sect. 2, we describe the physics included in the shearing coordinate frame used to numerically model the growth of the core. In Sect. 3 we present the results from our simulations and analyse the accretion rates for various core masses. We compare our results to analytic expressions capturing the essential physics underlying the phenomena at hand, namely the sub-Keplerian gas velocity, the particle size, the Keplerian shear and the gravitational pull from the seed core. The effect of local changes in the pressure gradient are analysed and we present the effect of including the backreaction of the particles on the gas flow. By extrapolating the measured accretion rates, we discuss the formation of gas and ice giant cores and derive a characteristic timescale for core formation by pebble accretion in Sect. 4. We discuss the approximations made in this paper and the limitations of our model in Sect. 5. Finally, in Sect. 6, we conclude that pebble accretion can explain rapid gas and ice giant formation in the core accretion scenario, even at wide stellar separations.

2. Physical model

The growth of a gas giant’s core occurs in a protoplanetary disc, a gaseous disc in the process of accreting onto the young star. Based on the mass distribution in the solar system and assuming a mean gas-to-dust ratio or metallicity of $\mathit{Z}\mathrm{=}\frac{{\mathrm{\Sigma }}_{\mathrm{p}}}{\mathrm{\Sigma }}\mathrm{=}\mathrm{0.01}\mathit{,}$(1)with Σ_p and Σ denoting the solid (dust+ice) and gas column densities, Hayashi (1981) constructed the minimum mass solar nebula (MMSN). He found the radial dependence of the gas column density to be $\mathrm{\Sigma }\mathrm{=}\mathrm{1700}{\left(\frac{\mathit{r}}{\mathrm{AU}}\right)}^{\mathrm{-}\mathrm{3}\mathit{/}\mathrm{2}}{\mathrm{g}\hspace{0.17em}\mathrm{cm}}^{-2}\mathit{,}$(2)with the orbital radius r within 0.35–36 AU. The thin disc is characterised by the ratio of the gas scale height H to the orbital distance $\frac{\mathit{H}}{\mathit{r}}\mathrm{=}\frac{{\mathit{c}}_{\mathrm{s}}}{{\mathit{v}}_{\mathrm{K}}}\mathrm{\approx }\mathrm{0.033}{\left(\frac{\mathit{r}}{\mathrm{AU}}\right)}^{\mathrm{1}\mathit{/}\mathrm{4}}\mathit{,}$(3)with c_s the sound speed of the gas and v_K the Keplerian velocity, ${\mathit{v}}_{\mathrm{K}}\mathrm{=}\mathit{r}{\mathit{\Omega }}_{\mathrm{K}}\mathrm{=}{\left(\frac{\mathit{GM}}{\mathit{r}}\right)}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{·}$(4)Here Ω_K is the Keplerian frequency.

Solids with radii smaller than the local mean free path of the gas, R ≤ (9/4)λ, are in the Epstein regime of gas-particle coupling (Epstein 1924). They react on a friction timescale t_f to changes between the relative particle velocity v and the local gas velocity u, $\begin{array}{ccc}{v̇}\mathrm{drag}\mathrm{=}\mathrm{-}\frac{\mathrm{1}}{{\mathit{t}}_{\mathrm{f}}}\mathrm{(}{v}\mathrm{-}{u}\mathrm{)}\mathrm{=}\mathrm{-}\frac{\mathit{\rho }{\mathit{c}}_{\mathrm{s}}}{{\mathit{\rho }}_{\mathrm{•}}\mathit{R}}\mathrm{(}{v}\mathrm{-}{u}\mathrm{)}\mathit{,}& & \end{array}$(5)where R and ρ_• are the radius and material density of the particle, while ρ is the local gas density. For particles in the vicinity of the midplane, with z < H, one can assume $\mathit{\rho H}\mathrm{\approx }\mathrm{\Sigma }\mathit{/}\sqrt{\mathrm{2}\mathit{\pi }}$, so that the particle size R in the MMSN can be recovered from its dimensionless friction time $\begin{array}{ccc}{\mathit{\tau }}_{\mathrm{f}}\mathrm{=}{\mathit{\Omega }}_{\mathrm{K}}{\mathit{t}}_{\mathrm{f}}& & \end{array}$(6)(also known as the Stokes number) as $\mathit{R}\mathrm{=}\mathrm{60}\hspace{0.17em}\mathrm{cm}\hspace{0.17em}{\mathit{\tau }}_{\mathrm{f}}{\left(\frac{{\mathit{\rho }}_{\mathrm{•}}}{\mathrm{2}{\hspace{0.17em}\mathrm{g}\hspace{0.17em}\mathrm{cm}}^{-3}}\right)}^{-1}{\left(\frac{\mathit{r}}{\mathrm{AU}}\right)}^{\mathrm{-}\mathrm{3}\mathit{/}\mathrm{2}}\mathrm{·}$(7)Figure 1 shows the relation between the orbital radius and the particle radius for different dimensionless friction times. Around 10 AU, a dimensionless friction time of τ_f = 0.1 corresponds to cm-sized particles, which we will refer to as pebbles. Close to the star, the gas density increases sufficiently for the particles to enter the Stokes drag regime, where ${\mathit{\tau }}_{\mathrm{f}}^{\mathrm{\left(}\mathrm{S}\mathrm{\right)}}\mathrm{=}\mathrm{(}\mathrm{4}\mathit{/}\mathrm{9}\mathrm{)}\mathrm{(}\mathit{R}\mathit{/}\mathit{\lambda }\mathrm{)}{\mathit{\tau }}_{\mathrm{f}}$ scales as  ∝ r^5/4. For a more complete description of different drag regimes, see e.g. Rafikov (2005) or Youdin (2010).

Fig. 1 The physical response of a particle reacting to the gas flow is set by the friction time tf. A particle of given size R has a dimensionless friction time τf = ΩKtf that depends on the orbital distance r. The red dash-dotted line marks the distance at which the particle size equals 9/4 of the mean free path λ of molecular hydrogen in the MMSN. Particles with R > (9/4)λ are located in the dashed region, and experience Stokes drag as opposed to Epstein drag. For the curves in the Stokes regime, we have ignored the transition into the non-linear Stokes regime, applicable for large particles close to the host star.

Fig. 2 Deviation Δ of the orbital velocity of a gas element with respect to an object orbiting with the full Keplerian frequency, normalised by the local sound speed, is plotted as function of orbital radius r in AU. The bold black line represents the traditional MMSN scaling, while the bold red line corresponds to the adapted MMSN as presented in Chiang & Youdin (2010). The shaded area connecting to thin curves indicates the effect of a strong pressure bump of strength δΔ =  −0.04. The adopted standard value of Δ = 0.05 in this paper, is accurate in a region around 5 AU, even without a strong pressure bump.

The gas component of the protoplanetary disc moves with a sub-Keplerian mean velocity, since the force due to the the solar gravity is reduced by the radially outwards pointing gas pressure force. The azimuthal velocity difference Δ = Δu_φ/c_s between the mean gas flow and a pure Keplerian orbit is given by $\begin{array}{ccc}\mathrm{\Delta }\mathrm{=}\mathit{\eta }\frac{{\mathit{v}}_{\mathrm{K}}}{{\mathit{c}}_{\mathrm{s}}}\mathrm{=}\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}\frac{{\mathit{c}}_{\mathrm{s}}}{{\mathit{v}}_{\mathrm{K}}}\frac{\mathit{\partial }\ln \mathrm{\left(}\mathit{P}\mathrm{\right)}}{\mathit{\partial }\ln \mathrm{\left(}\mathit{r}\mathrm{\right)}}\mathit{,}& & \end{array}$(8)where $\mathit{P}\mathrm{=}\mathit{\rho }{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}$ is the gas pressure and η is a measure of the gas pressure support (Nakagawa et al. 1986). In the MMSN, Δ has a weak radial dependency, Δ ≈ 0.05⁽r/AU^)1/4, as can be seen in Fig. 2. However, comparison of the MMSN model with observed protostellar accretion discs, (e.g. Bell et al. 1997) and studies of solar nebula metallicities (Lodders 2003) have prompted updated MMSN models, with a less steep pressure gradient, Δ = 0.036⁽r/AU^)2/7 (Chiang & Youdin 2010), as illustrated in Fig. 2.

The turbulent nature of an accreting protoplanetary disc can result in local pressure maxima (Johansen et al. 2009a; Fromang & Stone 2009). As can be seen from Eq. (8) these pressure bumps can locally reduce the headwind the pebbles experience. Reductions by δΔ ≈  − 0.02 are seen in shearing box simulations of the MRI (Johansen et al. 2009a; Fromang & Stone 2009) and global simulations (Lyra et al. 2008b). We have illustrated the effect of a strong pressure bump, with δΔ ≈  − 0.04, in Fig. 2.

Since particles face a headwind, they will drift radially and azimuthally as $\begin{array}{ccc}{\mathit{v}}_{\mathit{r}}& \mathrm{=}& \mathrm{-}\mathrm{2}\frac{{\mathit{\tau }}_{\mathrm{f}}}{{\mathit{\tau }}_{\mathrm{f}}^{\mathrm{2}}\mathrm{+}\mathrm{1}}\mathit{\eta }{\mathit{v}}_{\mathrm{K}}\mathit{,}\\ {\mathit{v}}_{\mathit{\phi }}& \mathrm{=}& \mathrm{-}\frac{\mathrm{1}}{{\mathit{\tau }}_{\mathrm{f}}^{\mathrm{2}}\mathrm{+}\mathrm{1}}\mathit{\eta }{\mathit{v}}_{\mathrm{K}}\mathit{,}\end{array}$as shown by Weidenschilling (1977) and Nakagawa et al. (1986). The total relative velocity between the particle and the core in pure Keplerian rotation is $\mathrm{\Delta }\mathit{v}\mathrm{=}\frac{\sqrt{\mathrm{4}{\mathit{\tau }}_{\mathrm{f}}^{\mathrm{2}}\mathrm{+}\mathrm{1}}}{{\mathit{\tau }}_{\mathrm{f}}^{\mathrm{2}}\mathrm{+}\mathrm{1}}\mathit{\eta }{\mathit{v}}_{\mathrm{K}}\mathit{,}$(11)which is well approximated by Δv ≈ ηv_k, since the particle sizes we consider, τ_f = (0.01,0.1,1), give us Δv/(ηv_K) = (1.0,1.0,1.1).

Fig. 3 Contour lines mark the core mass Mc (in units of Earth masses M ⊕ ), as function of the dimensionless mass unit ${\mathit{\mu }}_{\mathrm{c}}\mathrm{=}\mathrm{\left(}\mathit{G}{\mathit{M}}_{\mathrm{c}}\mathrm{\right)}\mathit{/}\mathrm{\left(}{\mathit{\Omega }}_{\mathrm{K}}^{\mathrm{2}}{\mathit{H}}^{\mathrm{3}}\mathrm{\right)}$ and orbital radius r in AU. The red contour line indicates the assumed minimal mass of a gas-giant core. The horizontal red dashed lines indicate the expected initial seed core mass from the streaming instability after planetesimal formation by gravitational collapse, μSI, and the transition mass μt (see Sect. 3.2).

Particles settle in the vertical direction (perpendicular to the orbital plane). The particle scale height H_p is a balance between midplane-directed gravity and turbulent diffusion parametrised by the coefficient δ_t (Youdin & Lithwick 2007), $\frac{{\mathit{H}}_{\mathrm{p}}}{\mathit{H}}\mathrm{\approx }\sqrt{\frac{{\mathit{\delta }}_{\mathrm{t}}}{{\mathit{\tau }}_{\mathrm{f}}}}\mathrm{\approx }\mathrm{0.01}\mathit{,}$(12)where we will make the approximation that this holds for the pebble size range we consider (τ_f = 0.01–1). The turbulence generated by the streaming instability (Youdin & Goodman 2005) self-regulates the particle midplane density to equal the gas density, independent of particle size. Since Z(H_p/H)^-1 ≈ ρ_p/ρ ≈ 1, a near-unity midplane density in a protoplanetary disc with metallicity Z = 0.01 sets the particle scale height to be H_p/H = 0.01. For turbulence generated through the MRI, a value of δ_t = 0.001 would give a ten times higher particle scale height, H_p/H ≈ 0.1, for particles of friction time τ_f = 0.1.

The aim of this paper is to investigate accretion onto cores of various masses, ranging from the expected initial masses of planetesimals to estimated final core masses of gas-giant planets. The core mass, or more precisely its gravitational parameter GM_c, is non-dimensionalised as ${\mathit{\mu }}_{\mathrm{c}}\mathrm{=}\frac{\mathit{G}{\mathit{M}}_{\mathrm{c}}}{{\mathit{\Omega }}_{\mathrm{K}}^{\mathrm{2}}{\mathit{H}}^{\mathrm{3}}}\mathrm{=}\frac{{\mathit{M}}_{\mathrm{c}}}{{\mathit{M}}_{\mathrm{\odot }}}{\left(\frac{\mathit{H}}{\mathit{r}}\right)}^{-3}\mathrm{\propto }{\mathit{r}}^{\mathrm{-}\mathrm{3}\mathit{/}\mathrm{4}}\mathit{,}$(13)with G the gravitational constant and M_⊙ the stellar mass. Figure 3 relates the dimensionless core mass μ_c to the orbital radius in the MMSN. Given the core mass, we can assign it an uncompressed radius of ${\mathit{R}}_{\mathrm{c}}\mathrm{=}\mathrm{890}{\left(\frac{\mathit{\rho }}{\mathrm{2}\hspace{0.17em}{\mathrm{g}\mathit{/}\mathrm{cm}}^{\mathrm{3}}}\right)}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{3}}{\left(\frac{{\mathit{M}}_{\mathrm{c}}}{{\mathrm{10}}^{-3}\hspace{0.17em}{\mathit{M}}_{\mathrm{\oplus }}}\right)}^{\mathrm{1}\mathit{/}\mathrm{3}}\mathrm{km}\mathit{.}$(14)The critical core mass for runaway accretion of a gaseous envelope is approximately 10 M _⊕  (Mizuno 1980), only weakly dependent on the orbital radius outside the terrestrial planet region (Rafikov 2006). At 5 AU this mass corresponds to μ_c ≈ 1, as can be seen in Fig. 3. However, as Hori & Ikoma (2011) point out, if the envelope can be significantly polluted by heavy elements from the accretion of icy bodies, the critical core mass will be reduced by up to two orders of magnitude.

Planetesimals are believed to have initial sizes in the 100−1000 km region (Johansen et al. 2007; Morbidelli et al. 2009). Johansen et al. (2012) find the characteristic clump mass by streaming instabilities to be μ_SI ≈ 5 × 10^-6 (see Fig. 3).

The dynamical equations of the particles, $\begin{array}{ccc}\frac{\mathrm{d}{v}}{\mathrm{d}\mathit{t}}\mathrm{=}\mathrm{-}\mathrm{2}{\Omega }\mathrm{K}\mathrm{\times }{v}\mathrm{+}\mathrm{3}{\mathit{\Omega }}_{\mathrm{K}}^{\mathrm{2}}\mathit{x}{e}\mathit{x}\mathrm{+}{g}\mathrm{c}\mathrm{-}\frac{\mathrm{1}}{{\mathit{t}}_{\mathrm{f}}}\mathrm{(}{v}\mathrm{-}{u}\mathrm{)}\mathit{,}& & \end{array}$(15)are solved with the Pencil Code¹ in the shearing box approximation (Goldreich & Tremaine 1980; Brandenburg et al. 1995). A Cartesian coordinate system is placed rotating at an arbitrary, but fixed orbital distance, with Keplerian frequency Ω_K. The x-axis points radially outwards, the azimuthal direction corresponds to the y-coordinate and the vertical z-direction is perpendicular to the midplane. The motion of the particles is described by Eq. (15), which includes the acceleration due to the core placed in the centre of the frame and the self-gravity of the particles solved for through the Poisson equation, $\mathrm{\nabla }\mathrm{·}{g}\mathrm{c}\mathrm{=}\mathrm{4}\mathit{\pi G}{\mathit{\rho }}_{\mathrm{p}}\mathit{.}$(16)Additionally, it includes the drag force term $\mathrm{-}\frac{\mathrm{1}}{{\mathit{t}}_{\mathrm{f}}}{}^{\mathrm{(}}{v}\mathrm{-}{u}{}^{\mathrm{)}}$, a term balancing the linearised gravity and the centrifugal force $\mathrm{3}{\mathit{\Omega }}_{\mathrm{K}}^{\mathrm{2}}\mathit{x}{e}{}_{\mathit{x}}$ and the Coriolis force  − 2Ω_K × v.

We perform simulations both with and without the backreaction term of the particles on the gas, with gas backreaction time (ρ_p/ρ)^-1t_f. The momentum equation for the fluid elements, $\begin{array}{ccc}\frac{\mathit{\partial }{u}}{\mathit{\partial t}}\mathrm{+}{u}\mathrm{·}\mathrm{\nabla }{u}& \mathrm{=}& \mathrm{-}\mathrm{2}{\Omega }\mathrm{K}\mathrm{\times }{u}\mathrm{+}\mathrm{3}{\mathit{\Omega }}_{\mathrm{K}}^{\mathrm{2}}\mathit{x}{e}\mathit{x}\mathrm{-}{\mathit{\Omega }}_{\mathrm{K}}^{\mathrm{2}}\mathit{z}{e}\mathit{z}\\ & & \end{array}$(17)includes the pressure gradient term,  − (1/ρ_p)∇P, and vertical gravity $\mathrm{-}{\mathit{\Omega }}_{\mathrm{K}}^{\mathrm{2}}\mathit{z}$. The continuity equation for the gaseous component of the protoplanetary disc is given by $\frac{\mathit{\partial \rho }}{\mathit{\partial t}}\mathrm{+}\mathrm{\nabla }\mathrm{·}\mathrm{\left(}\mathit{\rho }{u}\mathrm{\right)}\mathrm{=}\mathrm{0}\mathit{,}$(18)and we use artificial hyperdiffusivity for the gas to dissipate energy on the smallest scales.

We solve these equations non-dimensionalised by the Keplerian frequency Ω_K, the scale height H of the gas disc and ρ₀, the gas midplane density. This has the benefit that, when interpreting the normalised results, the orbital dependency is nearly fully recovered from these parameter’s MMSN orbital scalings (e.g. Figs. 1–3). The numerical results are however not completely scale-free. When solving for self-gravity with the Poisson equation, we must set $\mathit{\Gamma }\mathrm{=}\frac{\mathrm{4}\mathit{\pi G}{\mathit{\rho }}_{\mathrm{0}}}{{\mathit{\Omega }}_{\mathrm{K}}^{\mathrm{2}}}\mathit{,}$(19)the non-dimensionalised form of the gravity constant G, as an initial condition. It shows only a weak dependency on the orbital radius, $\mathit{\Gamma }\mathrm{\approx }\mathrm{0.04}{\left(\frac{\mathit{r}}{\mathrm{3} \mathrm{AU}}\right)}^{\mathrm{1}\mathit{/}\mathrm{4}}\mathit{,}$(20)in the MMSN. We therefore fix Γ = 0.04 for the remainder of this paper.

All simulations are performed in a three-dimensional shearing box, with a fixed particle scale height of H_p/H = 0.01. Run parameters for all simulations used in this paper can be inspected in Table 1. The core is fixed in the origin of the coordinate system, a valid approximation for the range of gas-decoupled seed masses we cover. When the escape velocity v_esc from the surface of the core is small compared the sound speed, variations in the gas density can be ignored, as can be seen from the hydrostatic equilibrium of an isothermal gas, $\frac{{\mathit{v}}_{\mathrm{esc}}^{\mathrm{2}}}{{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}}\mathrm{\approx }\frac{\mathit{\partial }\ln \mathit{\rho }}{\mathit{\partial }\ln \mathit{r}}\mathrm{·}$(21)We argue this approximation holds up to the largest cores we consider in Sect. 5. All simulations, with the exception of 1e-3_0.1_b (see Table 1), do include the gas drag on the particles, but lack the backreaction from the particles of the core on the gas. When omitting the backreaction term, the gas velocity equals to the sub-Keplerian velocity, u_y =  −ηv_K. However, simulation 1e-3_0.1_b includes the particle’s backreaction on the gas and follows the numerical scheme discussed in Youdin & Johansen (2007). When including the backreaction term, we also turn on the vertical gravity force for the particles, $\mathrm{-}{\mathit{\Omega }}_{\mathrm{K}}^{\mathrm{2}}\mathit{z}{e}{}_{\mathit{z}}$.

All runs have sheared periodic boundary conditions in the radial direction, but particles crossing azimuthal boundaries get removed from the simulation domain (with exception of 1e-3_0.1_b), in order to avoid accretion of particles already focused from their first passage past the core. We have also run simulations including collisions, with the scheme discussed in Johansen et al. (2012), and found no measurable difference in the accretion rates on the seed core masses. We therefore omit collisions from the simulations in this paper. Implications and limitations of the simulation set up are further discussed in Sect. 5.

3. Results

Inspection of the particle’s momentum equation, Eq. (15), reveals an important length scale. The Hill radius, ${\mathit{r}}_{\mathrm{H}}\mathrm{=}{\left(\frac{\mathit{G}{\mathit{M}}_{\mathrm{c}}}{\mathrm{3}{\mathit{\Omega }}^{\mathrm{2}}}\right)}^{\mathrm{1}\mathit{/}\mathrm{3}}\mathit{,}$(22)is set by the gravitational competition between the acceleration towards the core and the stellar tidal field in the radial direction. At a separation r_H from the core, the orbital time around the core approximately equals the orbital time around the star, $\mathrm{2}\mathit{\pi }{\mathit{\Omega }}_{\mathrm{K}}^{-1}$. The Hill sphere’s radius grows linearly with the orbital radius r_H ∝ r, placing more material in the gravitational region of influence of the core.

3.1. Drift accretion

When ignoring the stellar tidal field and the Coriolis force, the Bondi radius²${\mathit{r}}_{\mathrm{B}}\mathrm{=}\frac{\mathit{G}{\mathit{M}}_{\mathrm{c}}}{\mathrm{\Delta }{\mathit{v}}^{\mathrm{2}}}\mathit{,}$(23)marks the outer point at which particles approaching the core with relative velocity Δv get significantly gravitationally deflected ( ≳ 1 rad, e.g. Binney & Tremaine 1987). Ignoring the stellar tidal field is a valid approximation before the core mass grows to the point where the Bondi radius becomes comparable to the Hill radius (see Sect. 3.2), and we can associate a core mass, ${\mathit{M}}_{\mathrm{t}}\mathrm{=}\sqrt{\frac{\mathrm{1}}{\mathrm{3}}}\frac{\mathrm{\Delta }{\mathit{v}}^{\mathrm{3}}}{\mathit{G}{\mathit{\Omega }}_{\mathrm{K}}}\mathit{,}$(24)with this transition. Masses with M_c < M_t are in the drift regime and pebbles embedded in the gaseous disc approach the core with a mean velocity comparable to the gaseous headwind the core experiences, Δv ≈ ηv_K.

Depending on the balance between the gravitational attraction of the core and the drag force the pebbles experience, a particle can be pulled from the mean gas flow and accreted if it dissipates enough energy while being deflected. When the gas-free core-crossing time associated with the Bondi radius, ${\mathit{t}}_{\mathrm{B}}\mathrm{=}\frac{{\mathit{r}}_{\mathrm{B}}}{\mathrm{\Delta }\mathit{v}}\mathit{,}$(25)is similar to the friction time t_f, the drag force will cause all pebbles within the Bondi radius to spiral inwards. However the effective accretion radius, the drift radius r_d, shrinks with respect to the Bondi radius, when t_B ≈ t_f is not satisfied.

When t_B > t_f, the particle under consideration is strongly coupled to the gas. In this limit, only grazing particles deflected on timescales shorter than the friction time get pulled out of the flow. If we let g denote the gravitational attraction due to the core’s mass, the condition ${\mathit{t}}_{\mathrm{g}}\mathrm{=}\frac{\mathrm{\Delta }\mathit{v}}{\mathit{g}}\mathit{<}{\mathit{t}}_{\mathrm{f}}$(26)needs to be satisfied for accretion to occur. Since the deflection time t_g is given by Δvr²/(GM_c) = ⁽r/r_B⁾²t_B, the effective drift accretion radius is given by ${\mathit{r}}_{\mathrm{d}}\mathrm{=}{\left(\frac{{\mathit{t}}_{\mathrm{B}}}{{\mathit{t}}_{\mathrm{f}}}\right)}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}{\mathit{r}}_{\mathrm{B}}\mathit{,}$(27)in the strong coupling limit. This radius corresponds to the settling radius in Ormel & Klahr (2010) and is also equal to the radius found by Perets & Murray-Clay (2011) where drag forces shear apart bound binaries in the Epstein regime. We verified this power law by numerically integrating orbits of test particles in the 2-body problem including drag, $\begin{array}{ccc}\frac{\mathit{\partial }{\mathit{v}}_{\mathit{x}}\mathit{/}\mathrm{\Delta }\mathit{v}}{\mathit{\partial t}\mathit{/}{\mathit{t}}_{\mathrm{B}}}& \mathrm{=}& \mathrm{-}{\left(\frac{{\mathit{r}}_{\mathrm{B}}}{\mathit{r}}\right)}^{\mathrm{3}}\frac{\mathit{x}}{{\mathit{r}}_{\mathrm{B}}}\mathrm{-}\frac{{\mathit{t}}_{\mathrm{B}}}{{\mathit{t}}_{\mathrm{f}}}\frac{{\mathit{v}}_{\mathit{x}}}{\mathrm{\Delta }\mathit{v}}\mathit{,}\\ \frac{\mathit{\partial }{\mathit{v}}_{\mathit{y}}\mathit{/}\mathrm{\Delta }\mathit{v}}{\mathit{\partial t}\mathit{/}{\mathit{t}}_{\mathrm{B}}}& \mathrm{=}& \mathrm{-}{\left(\frac{{\mathit{r}}_{\mathrm{B}}}{\mathit{r}}\right)}^{\mathrm{3}}\frac{\mathit{y}}{{\mathit{r}}_{\mathrm{B}}}\mathrm{-}\frac{{\mathit{t}}_{\mathrm{B}}}{{\mathit{t}}_{\mathrm{f}}}\left(\frac{{\mathit{v}}_{\mathit{y}}}{\mathrm{\Delta }\mathit{v}}\mathrm{-}\mathrm{1}\right)\mathit{,}\end{array}$where we non-dimensionalised the particle equation of motion, ignoring disc dynamics. This is a valid approximation in the drift regime, where ${\mathit{t}}_{\mathrm{B}}\mathrm{\ll }{\mathit{\Omega }}_{\mathrm{K}}^{-1}$. Sample orbits can be investigated in the inset of Fig. 4, which shows the maximal particle-core separation leading to capture. The drift radius for strongly coupled particles falls of as  ∝ ${{}^{\mathrm{(}}\mathit{t}_{\mathrm{B}}\mathit{/}{\mathit{t}}_{\mathrm{f}}{\mathrm{)}}^{}}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}$ as predicted. Particles with t_f ≈ t_B, get efficiently accreted within a Bondi radius from the core.

Particles weakly coupled to the gas with respect to low-mass cores (t_B < t_f) are less aided by drag as they get deflected by the core. As seen in Fig. 4, a rapid fall-off occurs for particles with t_f ≈ 10²t_B. The orbits in the inset show these particles to be gravitationally scattered, similar to the case were no gas drag is present. Here, the physical radius of the core becomes relevant, since accretion now occurs through gravitational focusing of particles on the core’s surface, which we have not taken into account in Fig. 4.

Fig. 4 Accretion efficiency in the weak and strong coupling regime. When the Bondi time tB = GM/Δv3 is equal to the particle’s friction time tf, the drift accretion radius rd peaks and equals the Bondi radius rB. For a particle of fixed size, the ratio tB/tf on the horizontal axis increases as the core mass grows in time. When particles are strongly coupled to the gas (tB > tf), with respect to the gravitational attraction of the core, the drift radius decreases as ${\mathit{r}}_{\mathrm{d}}\mathrm{~}{{}^{\mathrm{(}}\mathit{t}_{\mathrm{B}}\mathit{/}{\mathit{t}}_{\mathrm{s}}{}^{\mathrm{)}}}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}$ (the analytical scaling of Eq. (27) is indicated with a full grey line). Near tB/tf ≈ 10-2 the drift radius rapidly decreases. The inset shows particle trajectories (grey curves) in this regime, which can be compared with those at tB/tf = 1 (black curves). Where the former are simply gravitationally deflected, in the latter case we see that particles inside the Bondi radius (marked by a red circle) are accreted by the central point source.

The accretion rate in the drift accretion regime is given by $\mathit{Ṁ}\mathrm{d}\mathrm{=}\mathit{\pi }{\mathit{\rho }}_{\mathrm{p}}{\mathit{r}}_{\mathrm{d}}^{\mathrm{2}}\mathrm{\Delta }\mathit{v,}$(30)when r_d is smaller than the particle scale height H_p. A representative simulation in this regime, performed with μ_c = 10^-5, is illustrated in Fig. 5. Pebbles drift with a sub-Keplerian velocity past the core and those entering the Bondi radius, here well inside the Hill radius, feed the growth of the embryo. Note that when r_d ≈ r_b, the core growth scales faster than exponential with mass, as $\mathit{Ṁ}{}_{\mathrm{d}}\mathrm{\propto }{\mathit{M}}_{\mathrm{c}}^{\mathrm{2}}$. Figure 6 shows the accretion rates calculated from simulations 1e-6_0.1, 2.5e-6_0.1 and 1e-5_0.1. Particles with friction time τ_f = 0.1 closely follow the maximal drift accretion efficiency, $\frac{\mathit{\mu ̇}\mathrm{d}}{{\mathit{\mu }}_{\mathrm{c}}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{4}}\frac{{\mathit{\rho }}_{\mathrm{p}}}{\mathit{\rho }}\frac{\mathit{\Gamma }{\mathit{\mu }}_{\mathrm{c}}}{{\mathrm{\Delta }}^{\mathrm{3}}}{\mathit{\Omega }}_{\mathrm{K}}\mathit{,}$(31)with r_d ≈ r_B. However the low-mass core in run 1e-6-0.1 comes close to the weak coupling limit and sees its accretion rate reduced.

Fig. 5 Accretion of pebbles with τf = 0.1 by the central core (μc = 10-5) in the drift regime. The colour coding shows the local particle surface density Σp, normalised by the average particle density ${}^{\mathrm{⟨}}\mathrm{\Sigma }_{\mathrm{p}}{\mathrm{⟩}}^{}$, in the simulated shearing box with Z = Hp/H = 0.01. Marked as a white dot, the central seed core can be seen in the first panel. Both the drift and Hill radii are plotted as white circles. The second panel illustrates the drift of the particles (Δ = 0.05) and creation of an accreting particle wake. In the third panel, particles within the drift radius rd ≈ rB are accreted. Particles further out may be carried out of the box by the sub-Keplerian gas or by the Keplerian shear. In the final panel the headwind has blown most pebbles past the core, with only a minority accreted.

We can envisage two effects reducing the accretion rate, if we were to continue to ignore the stellar tidal field even for higher mass cores. Firstly, when the core enters the strong coupling limit, growth slows down to exponential, $\mathit{Ṁ}{}_{\mathrm{d}}\mathrm{\propto }{\mathit{t}}_{\mathrm{B}}^{-1}{\mathit{r}}_{\mathrm{B}}^{\mathrm{2}}\mathrm{\propto }{\mathit{M}}_{\mathrm{c}}$. Secondly, when the accretion radius becomes comparable to the particle scale height, the appropriate expression for the accretion rate is given by $\mathit{Ṁ}\mathrm{d}\mathrm{=}\mathrm{2}{\mathit{r}}_{\mathrm{d}}{\mathrm{\Sigma }}_{\mathrm{p}}\mathrm{\Delta }\mathit{v,}$(32)where Σ_p is the particle column density. When r_d ≈ r_B, we get exponential growth Ṁ_d/M_c = 2Σ_p/Δv.

3.2. Hill accretion

When the core mass grows to the point where the Bondi radius ${\mathit{r}}_{\mathrm{B}}\mathrm{\propto }{\mathit{M}}_{\mathrm{c}}^{\mathrm{2}}$ is comparable to its Hill radius ${\mathit{r}}_{\mathrm{H}}\mathrm{\propto }{\mathit{M}}_{\mathrm{c}}^{\mathrm{1}\mathit{/}\mathrm{3}}$ (or identically v_H = Δv or ${\mathit{t}}_{\mathrm{B}}\mathit{/}{\mathit{t}}_{\mathrm{f}}\mathrm{=}{\mathit{\tau }}_{\mathrm{f}}^{-1}$), it will cross the transition mass, ${\mathit{M}}_{\mathrm{t}}\mathrm{=}\sqrt{\frac{\mathrm{1}}{\mathrm{3}}}\frac{\mathrm{\Delta }{\mathit{v}}^{\mathrm{3}}}{\mathit{G}{\mathit{\Omega }}_{\mathrm{K}}}\mathrm{\approx }\mathrm{3}\mathrm{\times }{\mathrm{10}}^{-3}{\left(\frac{\mathrm{\Delta }}{\mathrm{0.05}}\right)}^{\mathrm{3}}{\left(\frac{\mathit{r}}{\mathrm{5}\hspace{0.17em}\mathrm{AU}}\right)}^{\mathrm{3}\mathit{/}\mathrm{4}}{\mathit{M}}_{\mathrm{\oplus }}\mathit{,}$(33)defined earlier and see a change in pebble accretion mechanism. The dimensionless form of the transition mass, ${\mathit{\mu }}_{\mathrm{t}}\mathrm{=}\mathrm{7}\mathrm{\times }{\mathrm{10}}^{-5}{\left(\frac{\mathrm{\Delta }}{\mathrm{0.05}}\right)}^{\mathrm{3}}\mathit{,}$(34)

Fig. 6 Accretion rate as function of the initial core mass μ0. Theoretical curves for the branches corresponding to drift and Hill accretion are plotted as respectively dashed and full lines in grey. The grey full circle marks the transition mass. Black crosses represent the simulated results in a stratified shearing sheet, with Δ = 0.05. Triangles correspond to simulations with modified Δ, a grey upwards pointing triangle corresponds to Δ = 0.07 and a downwards pointing triangle corresponds to Δ = 0.03. In the Hill branch, the position of both triangles lie on top of the black crosses and are omitted for clarity. The grey square shows the result of simulation 1e-3_0.1_b, which includes the particle backreaction.

Fig. 7 Trajectories for particles with dimensionless friction time τf = 0.01,0.1,1 obtained from the 2D Hill equations including gas drag. Pebbles with τf = 0.1 and impact parameters below a Hill radius efficiently get accreted. Larger particles of τf = 1 are pulled in from wider separations, but cores lose particles on horseshoe orbits. Particles strongly coupled to the gas, with τf = 0.01, need close encounters well within the Hill sphere in order to fall onto the core.

Fig. 8 Normalised accretion rates, , for different particle sizes with friction time τf = 1 and τf = 0.01, as indicated by respectively asterisks and diamonds. The grey curves correspond to those shown in Fig. 6 and similarly the transition mass is indicated by a full grey circle. In the Hill branch, larger particles, τf = 1, get accreted as efficiently as particles with friction time τf = 0.1 (see Fig. 6), but in the drift branch they never get accreted at the full drift rate. On the other hand, small particles, τf = 0.01, get efficiently accreted in the drift branch, but less so in the Hill regime.

scales as the cube of the headwind parameter Δ (see Fig. 2). The Hill radius now sets the maximal impact parameter from which particles can be accreted. When M_c > M_t, pebbles at the edge of the Hill sphere approach the core with relative velocity ${\mathit{v}}_{\mathrm{H}}\mathrm{\equiv }{\mathit{\Omega }}_{\mathrm{K}}{\mathit{r}}_{\mathrm{H}}\mathit{.}$(35)The Keplerian shear v =  −(3/2)Ω_Kx dominates over the headwind in the Hill branch, since ${\mathit{v}}_{\mathrm{H}}\mathit{/}\mathrm{\Delta }\mathit{v}\mathrm{=}\sqrt{{\mathit{r}}_{\mathrm{B}}\mathit{/}{\mathit{r}}_{\mathrm{H}}}$. The inverse Keplerian frequency, ${\mathit{\Omega }}_{\mathrm{K}}^{-1}$, is the gravitational crossing timescale at the Hill radius, independent of core mass. For particles with friction times close to the orbital timescale (τ_f = 0.1–1), all particles entering the Hill sphere will be accrete, as illustrated in Fig. 7. Here we present particle orbits obtained from the Hill equations including drag $\begin{array}{ccc}\frac{\mathit{\partial }{\mathit{v}}_{\mathit{x}}\mathit{/}{\mathit{v}}_{\mathrm{H}}}{\mathit{\partial t}\mathit{/}{\mathit{\Omega }}_{\mathrm{K}}^{-1}}& \mathrm{=}& \mathrm{+}\mathrm{2}\frac{{\mathit{v}}_{\mathit{y}}}{{\mathit{v}}_{\mathrm{H}}}\mathrm{+}\mathrm{3}\frac{\mathit{x}}{{\mathit{r}}_{\mathrm{H}}}\mathrm{-}\mathrm{3}{\left(\frac{\mathit{r}}{{\mathit{r}}_{\mathrm{H}}}\right)}^{-3}\frac{\mathit{x}}{{\mathit{r}}_{\mathrm{H}}}\mathrm{-}\frac{\mathrm{1}}{{\mathit{\tau }}_{\mathrm{f}}}\frac{{\mathit{v}}_{\mathit{x}}}{{\mathit{v}}_{\mathrm{H}}}\\ \frac{\mathit{\partial }{\mathit{v}}_{\mathit{y}}\mathit{/}{\mathit{v}}_{\mathrm{H}}}{\mathit{\partial t}\mathit{/}{\mathit{\Omega }}_{\mathrm{K}}^{-1}}& \mathrm{=}& \end{array}$where $\mathit{r}\mathrm{=}\sqrt{{\mathit{x}}^{\mathrm{2}}\mathrm{+}{\mathit{y}}^{\mathrm{2}}}$ is the particle-core distance. Similar to the drift case (Sect. 3.1), when the gravitational deflection time (here independent of the core mass and  ~ Ω^-1) is similar to t_f, enough energy will be dissipated during the approach to mediate the accretion of the pebbles within the Hill sphere. The accretion rate is then given by $\mathit{Ṁ}\mathrm{H}\mathrm{=}\mathrm{2}{\mathit{r}}_{\mathrm{H}}{\mathrm{\Sigma }}_{\mathrm{p}}{\mathit{v}}_{\mathrm{H}}\mathrm{\propto }{\mathit{M}}_{\mathrm{c}}^{\mathrm{2}\mathit{/}\mathrm{3}}\mathit{,}$(38)since for the core masses under consideration, M_c > M_t, we accrete the total particle surface density, Σ_p (r_H > H_p). This growth mode is confirmed by our numerical simulations 1e-4_0.1, 1e-3_0.1, 1e-2_0.1, 1e-1_0.1 and 1e-3_1.0, 1e-2_1.0, 1e-1_1.0. They can be inspected in Figs. 6 and 8. An example of a simulation with a seed mass accreting at the Hill rate is illustrated in Fig. 9. Here, particles of τ_f = 1 entering the Hill sphere drive the growth of the core of mass μ_c = 10^-2. Accretion occurs through a particle disc, as was previously resolved in high-resolution 2D simulations by Johansen & Lacerda (2010).

Note that in the classical scenario of planetesimal accretion, one never captures objects from the full Hill sphere, but only by a fraction α^1/2r_H, with α ≈ r_c/r_H set by the physical radius of the core r_c. In the terminology of Rafikov (2011), slow accretion of planetesimals between the shear- and dispersion-dominated dynamical regime, from a part of the particle scale height H_p = v/Ω_K (e.g. Dodson-Robinson et al. 2009), goes as approximately $\mathit{Ṁ}\mathrm{\approx }\mathit{\pi \alpha }{\mathit{r}}_{\mathrm{H}}^{\mathrm{2}}{\mathit{\rho }}_{\mathrm{p}}\mathit{v}\mathrm{\approx }\mathit{\alpha }{\mathit{r}}_{\mathrm{H}}{\mathrm{\Sigma }}_{\mathrm{p}}{\mathit{v}}_{\mathrm{H}}\mathrm{\approx }\mathit{\alpha }\mathit{Ṁ}{}_{\mathrm{H}}$. At 5 AU, this gives a reduction in the accretion rate of α ≈ 10^-3(r/5   AU)^-1, for a standard solid density (Goldreich et al. 2004). Accretion of planetesimal fragments from a thin midplane, as discussed in Rafikov (2004) is more efficient, with the accretion rate being proportional to $\sqrt{\mathit{\alpha }}\mathit{Ṁ}{}_{\mathrm{H}}$. At 5 AU, this limits the growth by $\sqrt{\mathit{\alpha }}\mathrm{\approx }\mathrm{3}\mathrm{\times }{\mathrm{10}}^{-2}{\left(\mathit{r}\mathit{/}\mathrm{5}\hspace{0.17em}\mathrm{AU}\right)}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}$. Thus accretion of pebbles at a rate Ṁ_H from the entire Hill sphere is extremely efficient, compared to the classical gas-free case.

Fig. 9 When the core is massive enough, it can efficiently accrete particles with τf = 1 entering its Hill sphere (indicated by the white circle in the first panel). On this scale, the Keplerian shear dominates over the relative velocity difference between the gas and Keplerian velocity (Δ = 0.05). Accretion seems to occur through a particle disc, visible after a steady state has set in (panel 2 to 3). In the last panel, accretion and Keplerian shear have removed most particles in the box. The colour coding is similar to Fig. 5. The simulation was performed in a stratified shearing sheet box with Z = Hp/H = 0.01.

Only the smallest particles we consider, with τ_f = 0.01, have an accretion efficiency that is less than optimal in the Hill branch (Fig. 8). Similar to the case of the strongly coupled particles in the drift regime, accretion requires the gravitational deflection time to be shorter than the friction time, as previously expressed in Eq. (26). The relative velocity for particles approaching the Hill sphere is set by the Keplerian speed Δv ≈ Ωr. This allows us to rewrite the accretion criterion as $\mathit{\Omega r}\frac{{\mathit{r}}^{\mathrm{2}}}{\mathit{GM}}\mathit{<}{\mathit{t}}_{\mathrm{f}}\mathit{,}$(39)which gives us an effective accretion radius ${\mathit{r}}_{\mathrm{eff}}\lesssim {\mathit{\tau }}_{\mathrm{f}}^{\mathrm{1}\mathit{/}\mathrm{3}}{\mathit{r}}_{\mathrm{H}}\mathit{.}$(40)In this regime, $\mathit{Ṁ}{}_{\mathrm{H}\mathit{,}\mathrm{eff}}\mathrm{\propto }{\mathit{\tau }}_{\mathrm{f}}^{\mathrm{2}\mathit{/}\mathrm{3}}$, which compared to particles of τ_f = 0.1 would give a reduction of the accretion rate by  ≈ 0.2, as can be seen from comparing Figs. 6 and 8.

3.3. Influence of headwind reduction and particle feedback

In the above discussion, we have kept the relative velocity between core and the gas disc constant at Δ = 0.05. As previously mentioned, we have ignored the presence of pressure bumps, local extrema in the radial pressure force resulting in regions of a reduced headwind, as well as extreme orbital distances where Δ can change significantly. While the accretion rate Ṁ_H in the Hill regime is insensitive to Δ, the Bondi branch up to the transition mass M_t is not. Figure 6 illustrates the effect on changes in Δ for various core masses (Δ = 0.05 ± 0.02). For a core mass accreting approximately from the full Bondi branch, Eq. (31) indicates that the accretion rate will be modified by a factor (0.05/0.03)³ ≈ 5 for Δ = 0.03 compared to accretion with Δ = 0.05. This is in agreement with the measured accretion rates for μ = 2.5 × 10^-6 (2.5e-6_0.1_0.03). The increase is reduced for 1e-5_0.1_0.03 where r_B grows to R_H, and increased for 1e-6_0.1_0.03 where t_B/t_f grows sufficiently out of the weak coupling limit. Overall we see that even weak pressure bumps decreasing Δ by 0.02 lead to more rapid accretion. Vice versa, for Δ = 0.07 we expect a reduction of the accretion rate by a factor (0.05/0.07)³ ≈ 0.4 and we measure similar, but lower accretion rates (1e-5_0.1_0.03, 2.5e-6_0.1_0.03).

Fig. 10 Including the backreaction friction force on the gas, a core of μ = 10-3 accretes particles of τf = 0.1 (third and last panels), even in presence of turbulence caused by the streaming instability. The turbulence has first been given 20 orbits (~126 Ω-1) to saturate (first and second panels), before the core is inserted. From the full Hill sphere, rH/H ≈ 0.07, pebbles are attracted in a prograde motion to the core (third panel). As before, the colour bar in the first panel gives the particle surface density. The insets show the azimuthally averaged particle density and vertical extent of the particle layer. For clarity, in the insets the colour coding covers a twice as wide range in particle overdensity, compared to the surface density plots.

In simulation 1e-3_0.1_b, we have departed from a smooth gas velocity profile, by including friction on the gas and letting turbulence develop by the streaming instability. After approximately 20 orbits, we place the seed core mass in the centre of the simulated domain (Fig. 10). The measured accretion rate does not deviate measurably from the case not including particle backreaction, as can be seen in Fig. 6. This indicates that our results are robust for the Hill branch, even in a turbulent environment. More simulations should be carried out in the future to verify the validity in the drift regime. However, this requires very high resolution in order to resolve both the Bondi radius and the streaming instability wavelength simultaneously.

4. Implications for gas giant growth

Fig. 11 Core growth as function of time, plotted for various orbital distances (0.5,5 and 50 AU). The drift branch, marked by grey solid lines, assumes an initial core mass of M0 = 10-5M⊕ and Δ = 0.05. The drift growth continues until the transition mass Mt is reached (marked by a full grey dot). Accretion continues through the more efficient Hill branch, drawn in black. For clarity, we start the Hill growth from the transition mass at time t = 0 yr, instead of continuing from the time where drift accretion comes to a halt. The masses of Ceres and Pluto (located at respectively 2.7 and 39 AU) are marked on the vertical axis for reference. The grey dotted curves correspond to classical planetesimal accretion (PA), where the faster growth corresponds to 2D accretion of planetesimal fragments (Rafikov 2004) and the slower to 3D accretion of planetesimals (e.g. Dodson-Robinson et al. 2009). Note that drift accretion timescale at 50 AU takes more than 108 yr and its transition mass point is not plotted.

Having numerically confirmed the pebble accretion rates for low to high core masses, we can extrapolate our results and find the time necessary to grow a core to the critical mass needed to attract its gaseous envelope.

In the drift regime, the accretion rate implies a growth timescale of $\mathrm{\Delta }{\mathit{t}}_{\mathrm{d}}\mathrm{=}{\mathrm{\int }}_{{\mathit{M}}_{\mathrm{0}}}^{{\mathit{M}}_{\mathrm{c}}}\mathit{Ṁ}\begin{array}{}-1\\ \mathrm{d}\end{array}\mathrm{d}\mathit{m}\mathrm{\approx }\frac{\mathrm{\Delta }{\mathit{v}}^{\mathrm{3}}}{\mathit{\pi }{\mathit{\rho }}_{\mathrm{p}}{\mathit{G}}^{\mathrm{2}}}{\mathit{M}}_{\mathrm{0}}^{-1}$(41)to reach a core mass M_c from an initial seed mass M₀ ≪ M_t. This lower limit, since we assume optimal accretion from the full Bondi radius, sets the time until the accretion rate blows up hyperbolically as M_c ∝ (Δt_d − t)^-1 and thus does not depend on the final mass we wish to reach. However, from our numerical results we know that this growth is not sustained and turns off to the Hill branch. The characteristic time to reach the end point for drift accretion, the transition mass M_t, is given by $\mathrm{\Delta }{\mathit{t}}_{\mathrm{d}}\mathrm{\approx }\mathrm{8}\mathrm{\times }{\mathrm{10}}^{\mathrm{6}}{\left(\frac{\mathrm{\Delta }}{\mathrm{0.05}}\right)}^{\mathrm{3}}{\left(\frac{{\mathit{\rho }}_{\mathrm{p}}\mathit{/}{\mathit{\rho }}_{\mathrm{0}}}{\mathrm{0.01}}\right)}^{-1}{\left(\frac{{\mathit{M}}_{\mathrm{0}}}{{\mathrm{10}}^{-5}{\mathit{M}}_{\mathrm{\oplus }}}\right)}^{-1}{\left(\frac{\mathit{r}}{\mathrm{5}\hspace{0.17em}\mathrm{AU}}\right)}^{\mathrm{2}}\mathrm{yr}\mathit{,}$(42)which at 5 AU is comparable to the gas disc lifetime. If particles sediment to the midplane, the ratio of the particle to the gas density, ρ_p/ρ₀, would be of order unity. However, past 0.5 AU, a small seed core mass, M₀ = 10^-5M_⊕, accretes non-sedimented particles of size τ_f < 0.01 most efficiently. For reasonable values of the local headwind, Δ = 0.03–0.07, growth is too slow to form cores large enough to enter the Hill accretion regime, and without pressure bumps Δ is even larger at wide stellar separations.

Hill accretion on the other hand, has a growth timescale of $\mathrm{\Delta }{\mathit{t}}_{\mathrm{H}}\mathrm{=}{\mathrm{\int }}_{{\mathit{M}}_{\mathrm{t}}}^{{\mathit{M}}_{\mathrm{crit}}}\mathit{Ṁ}\begin{array}{}-1\\ \mathrm{H}\end{array}\mathrm{d}\mathit{m}\mathrm{\approx }\frac{{\mathrm{3}}^{\mathrm{5}\mathit{/}\mathrm{3}}{\mathit{\Omega }}_{\mathrm{K}}^{\mathrm{1}\mathit{/}\mathrm{3}}}{\mathrm{2}{\mathit{G}}^{\mathrm{2}\mathit{/}\mathrm{3}}{\mathrm{\Sigma }}_{\mathrm{p}}}{\mathit{M}}_{\mathrm{crit}}^{\mathrm{1}\mathit{/}\mathrm{3}}\mathit{,}$(43)which is only weakly dependent on the critical mass for gas envelope attraction M_crit and independent of the transition mass, when M_t ≪ M_crit. The core growth when accreting pebbles is fast in this regime, at 5 AU the critical mass is reached after $\mathrm{\Delta }{\mathit{t}}_{\mathrm{H}}\mathrm{\approx }\mathrm{4}\mathrm{\times }{\mathrm{10}}^{\mathrm{4}}{\left(\frac{{\mathit{M}}_{\mathrm{crit}}}{\mathrm{10}\hspace{0.17em}{\mathit{M}}_{\mathrm{\oplus }}}\right)}^{\mathrm{1}\mathit{/}\mathrm{3}}\left(\frac{\mathit{r}}{\mathrm{5}\hspace{0.17em}\mathrm{AU}}\right)\hspace{0.17em}\mathrm{yr}\mathit{.}$(44)Furthermore, the growth timescale Δt_H scales linearly with orbital distance r, as opposed to quadratic in the drift regime. This makes core formation possible in distant regions of the protoplanetary disc. Also, note that Hill accretion rate is maintained for a single particle size with friction times τ_f ~ 0.1–1, independent of the core mass, as opposed to the Bondi regime where one unrealistically needs to maintain t_B ≈ t_f to maintain the maximal accretion rate.

Figure 11 shows the core growth in both regimes and the dependency on the orbital distance. We conclude that fast core growth is possible through pebble accretion, provided that the initial seed mass for the core is above the local transition mass. A sufficiently large embryo can only be grown by drift accretion in pressure bumps with low headwind, Δ ≲ 0.05, or be the result of planetesimal formation by gravitational collapse after concentration by e.g. streaming instabilities (see also Fig. 3). It is interesting to note that both Ceres and Pluto have less than the critical mass needed for fast Hill accretion, which might explain why they failed to grow to gas or ice giants. Indeed, one can make the assumption that only those planetesimals that formed early enough in the high-mass tail of the initial planetesimal mass distribution could serve as the seed for gas-giant cores.

As an illustration of the rapid core growth by pebble accretion in the Hill regime, we compare it to the core growth time for planetesimal accretion in Fig. 12. As discussed in Sect. 3.2, the inability to accrete solids from the entire Hill sphere, as opposed to pebble accretion, leads to significant longer core formation times, in conflict with the observed dissipation time of protoplanetary discs.

Fig. 12 Time needed for core growth up to 10 M⊕ at various locations in the disc. The solid black line gives the formation time of the core for pebble accretion in the Hill regime, while grey lines give the time needed to form the critical 10-Earth-mass core by planetesimal accretion. The dashed grey line represents planetesimal fragment accretion from a thin midplane layer, as studied by Rafikov (2004). The red shaded area shows the approximate time interval in which the protoplanetary disc loses its gaseous component and encompasses for example the estimated age of gas giant LkCa 15b (Kraus & Ireland 2012). Core formation needs to occur before this time.

5. Discussion

We discuss here the assumptions and limitations of our results.

Midplane layer thickness. One component of the pebble accretion scenario is the presence of a thin particle disc (H_p = 0.01H). This low particle scale height is expected from turbulence driven by streaming instabilities, independent of particle size, as discussed in Sect. 2. A moderately higher particle scale height, as may be the case for turbulence caused by the magnetorotational instability, can result in a situation where r_H < H_p past the transition core mass. This would result in a temporarily reduced accretion rate, by a factor r_H/H_p = (H_p/H)^-1(r_H/H) = (1/3)^1/3(H_p/H)^-1μ^1/3, until the Hill radius grows beyond the particle scale height.

Particle size. The assumption of a single particle size in our simulations can be criticised, but as discussed in the introduction, observations of protoplanetary discs allow a large fraction of the solid mass to reside in the particle size range that we consider, τ_f = 0.01–1 (Wilner et al. 2005). A large abundance of particles larger than pebbles is not expected from coagulation models (Blum & Wurm 2008; Brauer et al. 2008; Windmark et al. 2012). However, as particles approach the core their icy component might sublimate; as friction would heat the particles, especially when a denser envelope starts forming around the core. It would be interesting to take this size-diminishing effect into account in a further investigation. On the other hand, particles might grow larger. In higher metallicity environments streaming instabilities become so effective in clumping solid material that one can fear particles to grow past the pebble size. However, we do not see this particle clumping in our simulations including the gas drag backreaction at the metallicity we consider (Z = 0.01). Strong clumping requires Z ≳ 0.02 (Johansen et al. 2009b; Bai & Stone 2010).

Gas structure. For the lower seed masses discussed in the paper, we previously argued (Sect. 2) that the gas density changes around the core are small. In the Hill regime, the ratio ${\mathit{v}}_{\mathrm{esc}}^{\mathrm{2}}\mathit{/}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathrm{\approx }\mathrm{2.3}\mathrm{\times }{\mathrm{10}}^{\mathrm{2}}{\mathit{\mu }}^{\mathrm{2}\mathit{/}\mathrm{3}}\mathrm{(}\mathit{r}\mathit{/}\mathrm{AU}\mathrm{)}$ (in the MMSN for standard solid density) can exceed unity for the highest core masses and the effects of an envelope should be taken into account. But, as also argued by Ormel & Klahr (2010), even if the direction of the flow moderately changes on scales within the Bondi radius due to stratification near the core, only particles with t_f ≪ t_B could be affected by it. Since these particles are too strongly coupled to the gas for accretion to take place in the first place (strong coupling limit), ignoring the core’s feedback on the gas is justified.

Keplerian orbits. In our analysis we assumed the core to be on a circular Keplerian orbit. The relative velocity between the core and the gas in Keplerian rotation could be significantly modified if competing cores would get excited by repeated close passages. However, as opposed to classical planetesimal growth, in our scenario gas damps the small particles and dynamical friction prevents the excitation of larger bodies, similar to the oligarchic growth regime. We do ignore gas-driven type-I migration of the core, important for core masses over 0.1 M _⊕  (Tanaka et al. 2002).

Random particle speed. In our simulations, particles approach the core in equilibrium with the gas flow. Particle interactions with the core last at most of the order ${\mathit{\Omega }}_{\mathrm{K}}^{-1}$, as in the Hill regime. Small particles (τ_f ≲ 1) are coupled to the gas on similar timescales. The passage of the core is quickly erased for the non-accreted particles, even when ignoring radial drift and turbulent diffusion. The core only catches up with the deflected particles after approximately ${\mathit{t}}_{\mathrm{pass}}\mathrm{=}\mathrm{2}\mathit{\pi r}\mathrm{\Delta }{\mathit{v}}^{-1}\mathrm{=}\mathrm{(}\mathrm{2}\mathit{\pi }\mathit{/}{\mathit{\Omega }}_{\mathrm{K}}\mathrm{)}\mathrm{(}\mathit{r}\mathit{/}\mathit{H}\mathrm{)}{\mathrm{\Delta }}^{-1}\mathrm{\approx }{\mathrm{10}}^{\mathrm{4}}{{}^{\mathrm{(}}\mathit{r}\mathit{/}{\mathrm{AU}}^{\mathrm{)}}}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}{\mathit{\Omega }}_{\mathrm{K}}^{-1}$. However, in the Hill accretion regime all particles that can be deflected are accreted, and it is the radial drift and diffusion of particles that fill up the feeding zone. Diffusion can be rapid, since the diffusion time associated with closing the Hill sphere ${\mathit{t}}_{\mathrm{H}}\mathit{\Omega }\mathrm{~}{\mathit{R}}_{\mathrm{H}}^{\mathrm{2}}\mathit{/}\mathrm{\left(}{\mathit{\delta }}_{\mathrm{t}}{\mathit{H}}^{\mathrm{2}}\mathrm{\right)}\mathrm{~}\mathrm{(}\mathrm{1}\mathit{/}\mathrm{3}{\mathrm{)}}^{\mathrm{2}\mathit{/}\mathrm{3}}{\mathit{\delta }}_{\mathrm{t}}^{-1}{\mathit{\mu }}^{\mathrm{2}\mathit{/}\mathrm{3}}$, is of order unity for a protoplanetary disc with δ_t = 0.01.

Particle drift. When the drag force responsible for radial drift is too small, particles could get trapped in mean motion resonances with the core. Weidenschilling & Davis (1985) studied large, τ_f ≥ 1, particles in the Stokes drag regime, and argued that particles smaller than these sizes feel large enough drag forces to escape resonant trapping around a Jupiter-mass planet at 5 AU. As shown by Tanaka & Ida (1997) inclusion of mutual planetesimal interactions breaks down the resonances, but dust gap formation still occurs for large planetesimals, where gas drag changes the semi-major axis of the the planetesimals after scattering with the protoplanet. The maximal particle size unaffected by particle trapping seems approximately inversely proportional to the planet’s mass, which is also seen in simulations performed by Paardekooper (2007). In fact, Weidenschilling & Davis (1985) argue that small pebbles are the only size that can be accreted by the core, since trapped larger planetesimals get dynamically excited and will be ground to fragments, which in their turn are capable of escaping the resonance. This picture is confirmed in simulations performed by Levison et al. (2010).

Dust gaps can open up before the core is massive enough to create a gap in the gas disc itself (Paardekooper & Mellema 2006). Muto & Inutsuka (2009) analytically show that the core has to be over a critical mass, ${\mathit{\mu }}_{\mathrm{c}}\mathit{>}\mathrm{\Delta }{\left(\frac{\mathit{H}}{\mathit{r}}\right)}^{-1}\mathrm{\approx }\mathrm{1}\mathit{,}$(45)for particles of τ_f ≤ 1 in order for a dust gap to emerge. Past r ≈ 1 AU, μ ≈ 1 is consistently above 10 M _⊕ , the critical core mass for gas and ice giants (see Fig. 3). Particles thus always drift radially fast enough to replenish the feeding zone of the core. Indeed, if the drift rate is set by R_d = 2πrΔvΣ_p, the requirement R_d ≥ Ṁ_H recovers the above criterion, Eq. (45). At the same time, as pointed out by Ormel & Kobayashi (2012), the particle drift can also be responsible of clearing up the entire reservoir of available pebbles in the disc.

Terrestrial planet formation. Growth at small orbital distances, r < 5 AU, is remarkably rapid in the pebble accretion model. Formation of rocky planets and possibly in situ formation of gas-giant planets in the terrestrial planet region seems problem-free from the perspective of the accretion rate. The growth timescales for both the drift and Hill accretion branch shrink to approximately 10⁵ yr at Earth-like separations from the host star. This could indicate that even terrestrial planet formation occurs rapidly during the gaseous disc phase. However, closer to the star the amount of material in an annulus of Hill-radius-width is small and the isolation mass by gap formation is lower. Also the optimally accreted particle size is large, around 10 cm, and ices are not available.

6. Conclusions

In this paper we have demonstrated that accretion of pebbles makes rapid formation of gas-giant cores possible. The growth timescale to reach the critical core mass for gas accretion is reduced by three orders of magnitude at 5 AU and four orders of magnitude at 50 AU, compared to the planetesimal accretion rate in between the shear- and dispersion-dominated dynamical regime. Compared to accretion of planetesimal fragments from a thin layer the formation time is shortened by approximately a factor 30 at 5 AU and a factor 100 at 50 AU. This is further support for the core accretion scenario, because cores can form by pebble accretion before gas dissipation after 1–10 million years, even at large orbital radii.

We can summarise the main numerical results as follows. Our simulations show gas drag to be a necessary ingredient for fast pebble accretion by the growing core. Omnipresent pebbles, particles with friction time around τ_f ≈ 0.1, are ideally suited for core growth. They are weakly enough bound to the gas to feel the gravitational pull from the core, but strongly enough to deposit their kinetic energy through drag forces, when passing the core. Low-mass cores, cores below the transition mass corresponding to a body of radius larger than approximately 1000 km, can accrete small particles drifting with the sub-Keplerian gas velocity past the core, but this process is slow, even in pressure bumps with reduced headwind and particle settling in a thin mid-plane layer. However growth in this regime could be important for a seed planetesimal formed just below the transition mass, where the accretion rate is high. Higher-mass cores can efficiently attract pebbles from the full Hill sphere, as was found by Johansen & Lacerda (2010) and Ormel & Klahr (2010). In this regime, the optimally accreted particle size is independent of the core mass, in contrast to the drift regime, where the particle size with the highest accretion rate increases linearly with radius of the growing core.

For the pebble accretion mechanism to be rapid, a significant fraction of the solid density needs to be in the form of pebbles close to the midplane, and some planetesimals need to form with sizes of 1000 km or larger. Theoretical models of planet formation show that these large seeds of approximately Ceres-size can form by self-gravity after clumping by the streaming instability (Johansen et al. 2007, 2012), so the pebble accretion scenario fits well with the formation of planetesimals by self-gravity.

The conditions for fast core growth are supported by observations. A large reservoir of pebbles is inferred in observations of many young protoplanetary discs (Wilner et al. 2005; Rodmann et al. 2006). Additionally, studies of the collisional evolution of the asteroid belt show that large asteroids must have formed early when gas was still present (Morbidelli et al. 2009). The early disappearance of mm-dust in protoplanetary discs on timescales shorter than 1 Myr (Lee et al. 2011) can be contributed to fast particle growth and rapid core formation. In fact, we see that gas giants form both rapidly (Kraus & Ireland 2012) and at large orbital radii (Marois et al. 2010). Our results predict that ice and gas giant planets, detectable with direct imaging surveys, will be abundant around young (~1 Myr) stars.

Further exploration of the pebble accretion mechanism by simulations with particle backreaction, around higher or lower metallicity discs in larger simulation domains, are needed to show the robustness of rapid core growth. Preferably, global simulations should be developed, including gap formation, radial drift, mean motion resonances and multiple cores, in order to get a full overview on the implications on fast core growth. When the core mass starts to approach the critical 10 Earth-masses for envelope attraction, we see a prograde particle disc emerge. Studying these discs, possibly the birthplaces of the regular satellites of gas giants, in conjunction with the emergence of a dense envelope around the core, will teach us more about the early growth and gas accretion thermodynamics of gas giants.

Pebble accretion provides us with a viable pathway to rapid formation of gas-giant and ice-giant cores. The necessary first step of the core accretion scenario can occur even at wide stellar separations, well within the lifetime of gaseous protoplanetary discs.

Acknowledgments

Computer simulations were performed at the Platon system of the Lunarc Centre for Scientific Computing at Lund University. Part of this work was supported by the Royal Physiographic Society. We thank Alessandro Morbidelli, Matthew Kenworthy and Chris Ormel for insightful discussions. The authors thank an anonymous referee for helpful comments.

References

All Tables

All Figures

Fig. 1 The physical response of a particle reacting to the gas flow is set by the friction time tf. A particle of given size R has a dimensionless friction time τf = ΩKtf that depends on the orbital distance r. The red dash-dotted line marks the distance at which the particle size equals 9/4 of the mean free path λ of molecular hydrogen in the MMSN. Particles with R > (9/4)λ are located in the dashed region, and experience Stokes drag as opposed to Epstein drag. For the curves in the Stokes regime, we have ignored the transition into the non-linear Stokes regime, applicable for large particles close to the host star.

In the text

Fig. 2 Deviation Δ of the orbital velocity of a gas element with respect to an object orbiting with the full Keplerian frequency, normalised by the local sound speed, is plotted as function of orbital radius r in AU. The bold black line represents the traditional MMSN scaling, while the bold red line corresponds to the adapted MMSN as presented in Chiang & Youdin (2010). The shaded area connecting to thin curves indicates the effect of a strong pressure bump of strength δΔ =  −0.04. The adopted standard value of Δ = 0.05 in this paper, is accurate in a region around 5 AU, even without a strong pressure bump.

In the text

Fig. 3 Contour lines mark the core mass Mc (in units of Earth masses M ⊕ ), as function of the dimensionless mass unit ${\mathit{\mu }}_{\mathrm{c}}\mathrm{=}\mathrm{\left(}\mathit{G}{\mathit{M}}_{\mathrm{c}}\mathrm{\right)}\mathit{/}\mathrm{\left(}{\mathit{\Omega }}_{\mathrm{K}}^{\mathrm{2}}{\mathit{H}}^{\mathrm{3}}\mathrm{\right)}$ and orbital radius r in AU. The red contour line indicates the assumed minimal mass of a gas-giant core. The horizontal red dashed lines indicate the expected initial seed core mass from the streaming instability after planetesimal formation by gravitational collapse, μSI, and the transition mass μt (see Sect. 3.2).

In the text

Fig. 4 Accretion efficiency in the weak and strong coupling regime. When the Bondi time tB = GM/Δv3 is equal to the particle’s friction time tf, the drift accretion radius rd peaks and equals the Bondi radius rB. For a particle of fixed size, the ratio tB/tf on the horizontal axis increases as the core mass grows in time. When particles are strongly coupled to the gas (tB > tf), with respect to the gravitational attraction of the core, the drift radius decreases as ${\mathit{r}}_{\mathrm{d}}\mathrm{~}{{}^{\mathrm{(}}\mathit{t}_{\mathrm{B}}\mathit{/}{\mathit{t}}_{\mathrm{s}}{}^{\mathrm{)}}}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}$ (the analytical scaling of Eq. (27) is indicated with a full grey line). Near tB/tf ≈ 10-2 the drift radius rapidly decreases. The inset shows particle trajectories (grey curves) in this regime, which can be compared with those at tB/tf = 1 (black curves). Where the former are simply gravitationally deflected, in the latter case we see that particles inside the Bondi radius (marked by a red circle) are accreted by the central point source.

In the text

Fig. 5 Accretion of pebbles with τf = 0.1 by the central core (μc = 10-5) in the drift regime. The colour coding shows the local particle surface density Σp, normalised by the average particle density ${}^{\mathrm{⟨}}\mathrm{\Sigma }_{\mathrm{p}}{\mathrm{⟩}}^{}$, in the simulated shearing box with Z = Hp/H = 0.01. Marked as a white dot, the central seed core can be seen in the first panel. Both the drift and Hill radii are plotted as white circles. The second panel illustrates the drift of the particles (Δ = 0.05) and creation of an accreting particle wake. In the third panel, particles within the drift radius rd ≈ rB are accreted. Particles further out may be carried out of the box by the sub-Keplerian gas or by the Keplerian shear. In the final panel the headwind has blown most pebbles past the core, with only a minority accreted.

In the text

Fig. 6 Accretion rate as function of the initial core mass μ0. Theoretical curves for the branches corresponding to drift and Hill accretion are plotted as respectively dashed and full lines in grey. The grey full circle marks the transition mass. Black crosses represent the simulated results in a stratified shearing sheet, with Δ = 0.05. Triangles correspond to simulations with modified Δ, a grey upwards pointing triangle corresponds to Δ = 0.07 and a downwards pointing triangle corresponds to Δ = 0.03. In the Hill branch, the position of both triangles lie on top of the black crosses and are omitted for clarity. The grey square shows the result of simulation 1e-3_0.1_b, which includes the particle backreaction.

In the text

Fig. 7 Trajectories for particles with dimensionless friction time τf = 0.01,0.1,1 obtained from the 2D Hill equations including gas drag. Pebbles with τf = 0.1 and impact parameters below a Hill radius efficiently get accreted. Larger particles of τf = 1 are pulled in from wider separations, but cores lose particles on horseshoe orbits. Particles strongly coupled to the gas, with τf = 0.01, need close encounters well within the Hill sphere in order to fall onto the core.

In the text

Fig. 8 Normalised accretion rates, , for different particle sizes with friction time τf = 1 and τf = 0.01, as indicated by respectively asterisks and diamonds. The grey curves correspond to those shown in Fig. 6 and similarly the transition mass is indicated by a full grey circle. In the Hill branch, larger particles, τf = 1, get accreted as efficiently as particles with friction time τf = 0.1 (see Fig. 6), but in the drift branch they never get accreted at the full drift rate. On the other hand, small particles, τf = 0.01, get efficiently accreted in the drift branch, but less so in the Hill regime.

In the text

Fig. 9 When the core is massive enough, it can efficiently accrete particles with τf = 1 entering its Hill sphere (indicated by the white circle in the first panel). On this scale, the Keplerian shear dominates over the relative velocity difference between the gas and Keplerian velocity (Δ = 0.05). Accretion seems to occur through a particle disc, visible after a steady state has set in (panel 2 to 3). In the last panel, accretion and Keplerian shear have removed most particles in the box. The colour coding is similar to Fig. 5. The simulation was performed in a stratified shearing sheet box with Z = Hp/H = 0.01.

In the text

Fig. 10 Including the backreaction friction force on the gas, a core of μ = 10-3 accretes particles of τf = 0.1 (third and last panels), even in presence of turbulence caused by the streaming instability. The turbulence has first been given 20 orbits (~126 Ω-1) to saturate (first and second panels), before the core is inserted. From the full Hill sphere, rH/H ≈ 0.07, pebbles are attracted in a prograde motion to the core (third panel). As before, the colour bar in the first panel gives the particle surface density. The insets show the azimuthally averaged particle density and vertical extent of the particle layer. For clarity, in the insets the colour coding covers a twice as wide range in particle overdensity, compared to the surface density plots.

In the text

Fig. 11 Core growth as function of time, plotted for various orbital distances (0.5,5 and 50 AU). The drift branch, marked by grey solid lines, assumes an initial core mass of M0 = 10-5M⊕ and Δ = 0.05. The drift growth continues until the transition mass Mt is reached (marked by a full grey dot). Accretion continues through the more efficient Hill branch, drawn in black. For clarity, we start the Hill growth from the transition mass at time t = 0 yr, instead of continuing from the time where drift accretion comes to a halt. The masses of Ceres and Pluto (located at respectively 2.7 and 39 AU) are marked on the vertical axis for reference. The grey dotted curves correspond to classical planetesimal accretion (PA), where the faster growth corresponds to 2D accretion of planetesimal fragments (Rafikov 2004) and the slower to 3D accretion of planetesimals (e.g. Dodson-Robinson et al. 2009). Note that drift accretion timescale at 50 AU takes more than 108 yr and its transition mass point is not plotted.

In the text

Fig. 12 Time needed for core growth up to 10 M⊕ at various locations in the disc. The solid black line gives the formation time of the core for pebble accretion in the Hill regime, while grey lines give the time needed to form the critical 10-Earth-mass core by planetesimal accretion. The dashed grey line represents planetesimal fragment accretion from a thin midplane layer, as studied by Rafikov (2004). The red shaded area shows the approximate time interval in which the protoplanetary disc loses its gaseous component and encompasses for example the estimated age of gas giant LkCa 15b (Kraus & Ireland 2012). Core formation needs to occur before this time.

In the text