© ESO 2018

1 Introduction

It is widely believed that turbulence plays an important role in the evolution of protoplanetary disks. For a long time, the magneto-rotational instability (MRI; Balbus & Hawley 1991) has been regarded as the leading candidate in driving the disk’s angular momentum transport. More recently, disk wind models have regained traction (Bai et al. 2016; Suzuki et al. 2016; Gressel 2017), where the turbulence in the midplane regions is limited to hydrodynamic instabilities such as the vertical shear instability (Nelson et al. 2013; Stoll & Kley 2014). Turbulence, in addition, is important in shaping the outcome of the early coagulation process. Already at low Mach numbers, turbulence dominates the relative velocity between particles (Völk et al. 1980; Ormel & Cuzzi 2007; Pan & Padoan 2010). It is also the only explanation for why we infer vertical structure (e.g., flared versus settled geometry), since turbulence allows small particles to be lifted from the disk midplane regions. Indeed, with the current state-of-the-art models small (micron-size) particles are produced in the midplane through collisions between pebbles and boulder sized particles before they diffuse upwards (Birnstiel et al. 2010, 2011; Krijt & Ciesla 2016).

Because of its subsonic nature, obtaining observational evidence of turbulence is hard. Disks like TW Hya, HD 163296, and DM Tau have been modeled by several groups (Hughes et al. 2011; Guilloteau et al. 2012; Flaherty et al. 2015, 2017) with turbulent Mach numbers inferred from the rather quiescent ~0.01 to the more vigorous ~0.1. However, it should be emphasized that constraining the turbulent rms velocity (σ) by these single-line profiles requires that the temperature profile be known to great precision. Generally, uncertainties affecting σ are limited by the absolute flux calibration and spectral resolving power (Teague et al. 2016). More indirect methods of obtaining σ employ the appearance of the dust disk in ALMA imagery. Applied to HL tau this indicates that the pebbles are settled into the midplane, resulting in a vertical turbulent diffusivity parameter α_z~10⁻⁴ (Pinte et al. 2016)¹. Flock et al. (2017) conclude this is consistent with a magnetized disk models that feature a “dead” midplane.

Both classical planetesimal-driven models for planet formation (Safronov 1969; Pollack et al. 1996) as well as the more recent pebble accretion model (Ormel & Klahr 2010; Lambrechts & Johansen 2012) are greatly affected by turbulence. Both models operate best under low-turbulence conditions. The runaway growth phase for the classical, planetesimal-driven accretion paradigm can only operate once the planetesimals start out with close to zerovelocity dispersions, but stochastic forcing by turbulence-triggered density fluctuations (Ida et al. 2008; Nelson & Gressel 2010; Gressel et al. 2011, 2012; Okuzumi & Ormel 2013) excites planetesimals to random velocities higher than their escape velocity. This implies that planetesimals have to be born large or that turbulence has to be weak (Ormel & Okuzumi 2013; Kobayashi et al. 2016). Similarly, the efficacy of pebble accretion to grow planets also depends on the turbulence. As pebbles will be stirred away from the midplane, it reduces the number of pebbles left to be accreted (Ormel & Klahr 2010; Guillot et al. 2014; Morbidelli et al. 2015). A second, less known, effect is that turbulent forcing may provide particles with an additional relative motion, which could also suppress accretion.

In this work, we consider simultaneously the effects of turbulent diffusivity and turbulent velocity. We do this by deriving a stochastic equation of motion (SEOM) for pebble-sized particles. Simply put, the SEOM is an extension of the Newtonian equation of motion, but with an additional stochastic component. For planets, stochastic forces have been invoked as a means to cross mean motion resonances (Rein & Papaloizou 2009; Paardekooper et al. 2013). As detailed by Rein & Papaloizou (2009), the stochastic motion is characterized by two parameters: the diffusivity D_P and the correlation time t_corr. The latter iscrudely the time over which the stochastic force changes its direction. For pebbles, we adopt a similar model, where now the stochastic motions are driven by aerodynamical coupling to the turbulent gas. However, in the few studies that have considered stochastic effects for pebble-sized particles, it is often assumed that turbulence does not feature a correlation time, i.e., t_corr is assumed less than any other timescale in the problem (Ciesla 2010; Zsom et al. 2011; Krijt & Ciesla 2016). This implies “white noise” behavior, i.e., that the particle is displaced in a random direction at every time. This approximation is known as the strong coupling limit (SCA).

A key goal of this paper is to test how the SCA fares in the light of the more accurate SEOM. We find that the SCA is generally applicable, as long as both the particle stopping time t_stop and the turbulent correlation time are sufficiently small. In addition, we will study the effect of a vertically varying turbulent gas diffusivity, D_zz(z), to investigate when it is viable to stir a fraction of pebble size particles to the disk surface.

Our main thrust will be to apply the SEOM and the SCA methods to calculate pebble accretion efficiencies in three-dimensional (3D) settings. In Liu & Ormel (2018, henceforth Paper I) we have defined ε as the probability that a pebble, drifting towards the star, will be accreted by a single planet(esimal)². A very small value of ε implies that a large number of pebbles are needed to grow the planet, while ε close to unity implies that pebble accretion is a very efficient accretion process. In Paper I we used planar 3-body calculation (star, planet, pebble) to calculate ε in two dimensions. We then investigated how this ε_2D changed as function of planet properties (mass and eccentricity), disk properties (position, radial drift velocity), and pebble properties (stopping time). In this work, we extend these calculation to the vertical dimension by including the planet’s inclination and disk turbulence. With the ε₃ D of this paper and the ε_2D of Paper I, we then obtain a general recipe for the pebble accretion efficiency (ε) of a single planet.

The plan of the paper is the following. In Sect. 2 we derive the SEOM. This section, as well as Appendix A, are rather technical and may be skipped by readers more interesting in the physical applications. In Sect. 3 we apply our newly developed SEOM to find vertical density distributions and show that our results are consistent with previous numerical and analytical studies. We apply the SEOM and SCA to pebble accretion in Sect. 4. We find ε for a variety of settings (planet mass and inclination, particle and disk properties) and present a framework to obtain ε under general circumstances (including the results found in Paper I). A comparison with previous studies is presented in Sect. 5. We summarize our findings in Sect. 6.

2 Model

2.1 Advection-diffusion equation

In this workwe model turbulence motion of particles and gas by an advection-diffusion equation: (1)

where ρ_P is the density of particles or gas species, v the systematic (drift) velocity, ρ_gas the gas density, and the particle diffusivity tensor whose elements are denoted D_ij. Importantly, the diffusion term acts on the gradient of the concentration (ρ_P ∕ρ_gas): it tends to erase concentration gradients and vanishes when the concentration is uniform.

In this work, we will restrict diffusion to operate only in the vertical (z) direction, considering only D_P,zz. Furthermore, we adopt the vertically isothermal solution for the gas density (2)

where Σ_gas is the gas surface density and H_gas the pressure scaleheight. Under these conditions, Eq. (1) can be manipulated: (3)

(Ciesla 2010)³. For small particles (including pebbles), the vertical velocity v_z equals the settling velocity, v_z = −zΩ²t_stop with t_stop the stopping time and Ω the Keplerian orbital frequency. The gas density no longer appears in Eq. (3), but the diffusivity appears at two places. On the RHS the diffusive term is responsible for spreading the particle concentration, resulting in a broader distribution of ρ_P (z). However, the additional advection term – a consequence of imposing Eq. (2) – counteracts this, enforcing the particle layer to remain stratified with a finite dispersion at all times (Ciesla 2010).

For a distribution of particles ρ_P(z)dz∕Σ gives the fraction of the particles within the interval [z, z + dz]. For a single particle, P(z) = ρ(z)∕Σ similarly denotes the probability of finding the particle within [z, z + dz] where P(z) is the probability density. We will use this identification below to obtain the correct, statistical properties of our single-particle (Lagrangian) stochastic model.

2.2 Stochastic equation of motion (SEOM)

The stochastic equation of motion is given by the following set of stochastic differential equations (SDEs): (4a) (4b) (4c)

where x is the position andv is the velocity of a particle. The particle is subject to gravitational forces F_g and gas drag forces. The latter have been expressed in terms of the stopping time, Δv∕t_stop, where Δ v is the relative gas-particle velocity. In Eq. (4c), W_t denotes a Wiener process (Brownian motion) and the corresponding differential is where is the normal distribution with zero mean and unity variance.

Apart from v, three velocity terms appear in Eq. (4b):

• 1.

  A laminar gas velocity v_gas. In our case, the gas velocity operates in the azimuthal direction

  (5)

  where and η represents the disk radial pressure gradient:(6)

  is assumed constant (Nakagawa et al. 1986).

• 2.

  A turbulent velocity v_turb. In Eq. (4b) this has been written in terms of an rms value () and a non-dimensional stochastic variable ζ_t. Here t_corr and D_zz are, respectively, the correlation time and diffusivity of the turbulent gas. In Eq. (4b) the turbulent forcing acts only in the vertical dimension.

• 3.

  A correction term v_hs

  (7)

  where . This is needed to enforce that Eq. (4) satisfies the hydrostatic balance condition, which assumption has entered the advection-diffusion equation (Eq. (3)). It also accounts for spatial gradients in D_zz. We derive it below.

Finally, Eq. (4c) is a stochastic differential equation (SDE) about a quantity ζ_t. This can be thought of as the normalized strength of the turbulent velocity. Specifically, Eq. (4c) describes an Ornstein–Uhlenbeck process (Uhlenbeck & Ornstein 1930) with zero mean (⟨ζ_t ⟩ = 0), unity variance () and correlation time t_corr. On long timescales ζ_t will be normally distributed; events separated by Δt ≫ t_corr are uncorrelated. However, times separated by Δt ≪ t_corr will feature a similar value of ζ_t and hence a similar turbulent gas velocity.

2.3 Turbulent correlation time

In Eq. (4) we are at liberty to choose t_corr, which can be identified with the correlation time (or lifetime) of the turbulent eddies. A smaller t_corr (while keeping D_zz fixed) implies a more vigorous turbulent forcing (larger σ_z), while a long t_corr implies that the turbulence is characterized by weaker but larger and longer-lived eddies. It is customary to adopt the Shakura & Sunyaev (1973) α-parameterization for the turbulent viscosity: (8)

We will adopt a similar parameterization for the gas diffusivity, i.e., where α_z reflects the diffusivity of the gas, not angular momentum transport. Using , the turbulent rms velocity becomes (9)

Following Dubrulle et al. (1995), Cuzzi et al. (2001) and Johansen et al. (2006), we usually adopt t_corr = Ω⁻¹ and hence . In Sect. 3.2 we also consider models where t_corr is longer.

The following qualifications will be adopted towards the turbulence strength:

• laminar for α_z = 0;

• weakly turbulent for α_z < 10⁻⁴;

• moderately turbulent for 10⁻⁴ < α_z < 10⁻²;

• strongly turbulent for α_z > 10⁻².

2.4 Strong coupling approximation (SCA)

The strong coupling applies when t_stop is small. In that limit, it can be shown that Eq. (4) simplifies into a single SDE for the position: (10)

(see Appendix A for the derivation). In Eq. (10), we have allowed D_zz to depend on position, which would give rise to an additional advective term (the fourth term in the square brackets). Equation (10) is analogous to Smoluchowski (1916) equation for the stochastic motion of a massless particle subject to a fluctuating force.

We arenow in a position to obtain the hydrostatic correction term v_hs. SDEs of the form (11)

can equivalently be cast in terms of an equation for the evolution of the probability density P(x, t) – i.e., a Fokker–Planck equation (12)

(van Kampen 1992)⁴. Applied to Eq. (10), the Fokker–Planck equation for the probability density P(z, t) reads (13)

where we consider the limit t_stop ≪ t_corr. The RHS of Eq. (13) can be expanded as . Identifying the probability density P(z) with the density ρ_P of Eq. (3), F_zt_stop with v_z, and using that D_P,zz = D_zz for strongly coupled particles (Völk et al. 1980), we obtain the hydrostatic correction term, Eq. (7). With this correction term, the SCA for the particle position in 1D reads (14)

as was already derived by Ciesla (2010) and also used in Zsom et al. (2011) and Charnoz et al. (2011). Comparing the second and third terms on the RHS, we obtain that the turbulent gradient effect becomes important when D_zz changes on scales less than .

We reflect on our findings. Equation (10) with v_hs equal to Eq. (7) describesthe stochastic motion of a particle experiencing drag and turbulent forces, with the turbulence characterized by a correlation time t_corr and a (possibly spatially dependent) gas diffusivity D_zz. Under the assumption of small t_stop and small t_corr, we obtain Eq. (14), consistent with Ciesla (2010). However, these equations provide no model for the particle velocity; and they will fail when the SCA conditions no longer materialize (long t_stop or long t_corr) – i.e., when the particle’s inertia matters. In these cases, Eq. (4) provides a more general description of stochastic motion of pebble-sized particles.

3 Vertical diffusion

We test our algorithms – the stochastic equation of motion (SEOM; Eq. (4)) and the strong coupling approximation (SCA; Eq. (14)) – for tracer particles (t_stop = 0) in Sect. 3.1 and massive particles (Sect. 3.2) for a variety of stopping times and α_z.

3.1 Tracer particle

Tracer particles should have a vertical distribution identical to the gas, Eq. (2). Because t_stop = 0 for tracer particles, the SEOM, as described in Eq. (4), contains a singularity and is not applicable. Therefore, we adopt the SCA method. In Eq. (10) we take t_stop = 0, η = 0, . The choices for H_gas and Ω are arbitrary.

In Fig. 1a we show the distribution of the vertical position of a single particle for α_z = 10⁻⁴, 10⁻³ and 10⁻². These have been obtained by storing the vertical positions after every 1 Ω⁻¹ for a total time of t_max = 10⁵ Ω⁻¹. In addition, we plot the expected distribution according to Eq. (2). The distributions are normalized such that they integrate to unity. Clearly, the distributions among the α_z differ, with α_z = 10⁻² best matching the expected normal distribution α_z = 10⁻⁴ the worst. The origin of these differences is the different number of independent samples that are obtained among the α_z. Since our sampling time is only Δ t = 1 Ω⁻¹, much smaller than the diffusion time t_diff = 1∕α_zΩ, sequential samples (in time) will be strongly correlated; only t_max∕t_diff = 10⁵α_z samples will be independent. Hence, the higher α, the better the correspondence to a Gaussian distribution.

Fig. 1 Top panel: normalized distributions of the vertical position z for the integration of the tracer case, using the strong coupling approximation (SCA) with tstop = 0. The vertical height z is recorded after every 1 Ω−1 for t = 105 Ω−1. Bars give the simulated distribution while the analytic – normal – distribution is shown by the black dashed line. Bottom panel: temporal evolution of the vertical position for αz = 10−2 (black) and αz = 10−2 (blue).

3.2 Massive particle

Next we consider the vertical distribution of massive particles (t_stop > 0), obtained bythe SEOM and the SCA methods. This means that in Eq. (4) we put F_g = −Ω²ze_z. We further take t_corr = 1 Ω⁻¹ unless mentioned otherwise. Integrations ran for 10⁶ Ω⁻¹ and the sampling period was Δt = 10 Ω⁻¹. The results are shown in Fig. 2. We fix α_z = 10⁻² but vary τ_s = t_stopΩ = 10⁻² (left panel), 10⁻¹ (center panel), and 10⁰ (right panel).

In all panels we compare the results for the SEOM (gray curves) with the SCA of Eq. (10) (blue). The thin black curve corresponds to a normal distribution with pebble aspect ratio: (15)

(Youdin & Lithwick 2007)⁵, where h_gas = H_gas∕r and (16)

In the limit of Ωt_corr ≪ 1 or , ξ ≈ 1 and the pebble aspect ratio reduces to (Dubrulle et al. 1995).

For small stopping times (τ_s = 10⁻²; left panel) the SCA and the SEOM overlap. From a numerical perspective the SCA is preferable as it is computationally much less intensive than the SEOM-method. For τ_s = 10⁻¹, the distributions slightly differ, as can best be seen from the lower panels. For τ_s = 1 particles the differences between the two methods amount to several tens of percents at z = 0, while towards the tails of the distribution the relative difference is larger even. The SEOM, however, is in perfect agreement with the Youdin & Lithwick (2007) theory on diffusive transport. The reason is that, like Youdin & Lithwick (2007), the SEOM accounts for the vertical oscillation (epicyclic motion) of particles, whereas the SCA does not. The SCA does not account for the pebble’s inertia and also does not involve a correlation time. From Eq. (15), we deduce that the SCA becomes invalid for .

In the above, we assumed that t_corr ≈Ω⁻¹, which is applicable in the ideal limit of MRI-turbulence (Sano et al. 2004; Johansen et al. 2006; Carballido et al. 2011). However, for non-ideal effects as ambipolar diffusion, the correlation time is expected to be longer (Bai & Stone 2011; Zhu et al. 2015). From Eq. (15), it is clear that pebbles will be more strongly stratified for large values of the turbulent correlation time t_corr. In Fig. 2b, a case with a 10 times longer correlation times (but still the same α_z = 10⁻² vertical diffusivity, implying larger, longer-lived but less vigorous eddies) is presented. Its vertical distribution is much narrower than the canonical t_corr = 1 Ω⁻¹ turbulence. A case with t_corr = 10² Ω⁻¹ is also shown in Fig. 2c. Clearly, a degeneracy between turbulent correlation time t_corr, diffusivity (α_z), and stopping time (t_stop) is present.

Fig. 2 Vertical distribution for particles of different stopping times: τs = tstopΩ = 10−2 (panel a), 10−1 (panel b), and 1 (panel c). Histograms plot the numerically obtained distributions with the stochastic equation of motion (SEOM; gray and red) and strong coupling approximation (SCA; blue) methods. Long tcorr runs are shown with red histograms (the tcorr = 102 Ω−1, τs = 0.1 run is displayed in panel c). Thin curves gives the normal distribution with the scaleheight of Eq. (15) (Youdin & Lithwick 2007). Note the different scaling among the panels.

3.3 Vertical gradient in the diffusivity

As an application of a more convoluted model, we consider a vertically dependent diffusion. In terms of α_z we adopt (17)

where α_mid is the diffusivity in the midplane and α_surface is the diffusivity in the upper regions. Such layered accretion (Gammie 1996) when the turbulence is confined to the upper regions, although our parameterization in Eq. (17) is completely arbitrary. We choose α_surface = 0.1 and α_mid = 10⁻³. This profile is plotted in Fig. 3 by the red curve. We further choose τ_s = 10⁻², such that τ_s ∕α_z > 1 in the midplane (indicating settling) and τ_s ∕α_z < 1 in the upper regions (indicating coupling to the gas).

In Fig. 3 the gray and black histograms show the normalized distribution of z obtained with the SCA and SEOM methods, respectively. The methods give consistent result. Clearly, the high |z| regions are sparsely sampled as the probability to find a particle at these heights is low. However, the fact that pebbles can be stirred to these heights at all may be surprising given the low α_mid. This is illustrated with the dark blue curve in Fig. 3, which plots P(z, h_mid) where . In fact, pebbles at |z| ≳ 2H_gas follow the gas distribution (light blue curve). Altogether P(z) can be approximated as the sum of two distributions. The majority of the pebbles (≈90%) follow the midplane distributions (given by h_mid), but about 10% of the pebbles follow the distribution given by the gas scaleheight.

Fig. 3 Density distribution P(z) obtainedfrom a vertically varying diffusivity. The diffusivity profile in terms of αz is given by the red curve. The particle stopping time is τs = 10−2 and is taken independent of z. The distributions obtained from integrating the stochastic equation of motions (Eq. (4); black histogram) and the strong coupling approximation (Eq. (10); gray histograms) are consistent. The dark blue curve gives the normal distribution, for the particle scaleheight evaluated in the midplane (i.e., Eq. (15) with α = 10−3). The light blue curves gives the gas distribution, scaled by a factor 0.1.

3.4 Local replenishment of small grains?

The ability of turbulence to stir ~mm-sized pebbles from the midplane to many gas scaleheights may offer an explanation for the persistent presence of small, (sub)micron size particles in the disk surface, as deduced from near-IR observations (e.g., Juhász et al. 2010). From a theoretical perspective, the presence of small particles is problematic as they should quickly coagulate among themselves and then settle to the disk midplane (Nakagawa et al. 1986; Tanaka et al. 2005; Dullemond & Dominik 2005). This implies that the grains are replenished. The most common theory postulates that the replenishment occurs in the disk midplane regions. Here grains are produced by high-velocity collisions among pebble-sized particles, which are subsequently transported (by diffusion) to the disk surface (Birnstiel et al. 2010). However, this is a rather indirect route to replenish small grains. First, it is doubtful if pebbles in the midplane will fragment; the gas may not be sufficiently turbulent. Second, it takes grains a time ~ 1∕α_zΩ to diffuse, which is rather long, again when α_z is small; these small grains may simply collide before reaching the surface (Krijt & Ciesla 2016).

Alternatively, in a disk with a suitable diffusivity profile (α_z (z)), it is possible to diffuse a small number of pebbles to the disk surface. There, due to the much stronger turbulent velocity field as compared to the midplane, collisions will undoubtedly be catastrophic. To cement these ideas, a coupled transport-collision/fragmentation model need to be considered.

4 3D pebble accretion

Following Paper I we calculate the accretion efficiency (ϵ) by conducting a series of N-body integrations to follow the trajectory of pebbles as they drift from orbits exterior to the planet to orbits interior to it. The pebble accretion efficiency ϵ is then found simply by counting the fraction of particles that settle to the planet. While Paper I investigated the role of the planet’s eccentricity, we fix e_p = 0 here and instead investigate the role of turbulence and the planet inclination (Sect. 4.5). To this effect, we let the pebble experience a stochastic motion in the vertical direction, as outlined in Eq. (4). The initial vertical position z_ini is given by the steady-state distribution characterized by the pebble scaleheight h_P. The initial radial position is set by (18)

where is the impact parameter for pebble accretion in the Hill regime and Δr_syn(τ_s, η, N) the distance a pebble drifts after N synodical orbits⁶. We choose N = 3. When the pebble radius has drifted to a distance r_p − r_Hill, we stop the calculation. It is then counted as a miss.

The fact that particles are only kicked in the z-direction, allows us to restrict the computations to a narrow ring. In contrast, when we would have considered the general case (turbulence operating inall dimensions), a much larger computational domain would be required because of the possibility of multiple encounters (similar to the eccentric case in Paper I). This complication is the key reason why we consider only turbulence in the vertical dimension.

Similar to Paper I, we express lengths in terms of the disk radius r_p and times in units of Ω⁻¹.

The keyparameters are as follows:

• q_p = M_p∕M_⋆ the planet-to-stellar mass ratio;

• h_gas = H_gas∕r the disk aspect ratio at the location of the planet;

• η, a measure of the radial drift velocity of the pebbles (Eq. (6));

• τ_s = t_stopΩ, the dimensionless stopping time;

• α_z, a proxy for the gas vertical diffusivity D_zz;

• t_corr, the turbulent correlation time, which together with α_z determines the magnitude of the turbulent rms velocities, Eq. (9).

Analogous to Paper I, we consider two integration methods:

• The SEOM, which integrates the particle velocity (Eq. (4)).

• The hybrid method, which uses the SCA, but switches to the SEOM when the particle is in the vicinity of the planet. Here we take the criterion to switch to the SEOM to be , where Δx and Δy are the distances between the planet and pebble in the x and y Cartesian coordinates. Note the absence of the vertical position in this criterion. The reason is that for some parameter combinations (small planet mass and high α_z) the vertical step size can easily become larger than the Hill radius in the SCA method.

As explained in Paper I, the SCA assumes that gas and particles are well coupled. It does not capture effects that take place on timescales Δ t less than t_stop + t_corr. On small Δt the particle moves ballistically, while the SCA model keep exhibiting random walk behavior at all (time)scales. Within the Hill sphere, where the numerical timestep will become small, the particle trajectories are therefore incorrect. The strong fluctuations in velocity space also complicates the numerical integration.

Different from Paper I, we do not account for the ballistic regime. Ballistic encounters are encounters in which the pebble is not captured by the settling mechanism, but where accretion occurs by virtue of the pebble hitting the surface of the target. Computationally, we can easily distinguish between ballistic and settling encounters by assigning an arbitrary small physical size to the planet (while keeping its mass). All accretion then occurs through the settling mechanism.

4.1 Standard model and analytical fits

For the standard model we take q_p = 3 × 10⁻⁷ (a 0.1 M_⊕ mass planet for a solar-mass star), η = 10⁻³, h_gas = 0.03, t_corr = Ω⁻¹ and vary α_z and τ_s. The choices for η and h_gas approximately correspond to a disk location of 1 au; both values will generally be higher in the outer disk. In Fig. 4 symbols give the mean value of ε obtained from our numerical integrations, ε = N_set∕N_tot, where out ofN_tot integrations N_set pebbles settled to the planet. Error bars correspond to the Poisson error, . The pebble accretion efficiency can be converted into the pebble growth mass, M_P,grw, defined as (21)

(values labeled on the right y-axis). This is the amount of pebbles needed to e-fold the planet’s mass. Finally, solid curves gives our fit to the data, which will be discussed in the subsequent sections and summarized in Sect. 4.7. The fitting expression is appropriate only for τ_s ≲ 1.

Clearly, where ε is higher, fewer integrations are needed to obtain a minimum signal-to-noise. In our integrations N_tot is not fixed, but different for each run in order to obtain a certain signal-to-noise ratio, which is determined by the number of settling encounters. Hence, most of the computational effort is spent in the runs where ε is small. Results of the SEOM method are shown by the open circles, while the crosses (slightly offset) denote the results from the hybrid method. The two methods give consistent results (see also Paper I). For small τ_s the SEOM method becomes computationally inefficient as it takes the pebble a long time to drift to the interior disk and the integration becomes very stiff. The hybrid method removes this latter problem and is the method of choice for small τ_s. Remarkably, the hybrid method gives acceptable results up to τ_s = 0.5.

Results for the non-turbulent (2D) limit (α_z = 0; black and gray symbols) were already discussed in Paper I. The efficiency decreases towards increasing particle stopping times – τ_s = 1 particles are accreted at the lowest efficiency – because the faster drift by the higher τ_s particles outweighs the larger linear cross section. Consequently, the pebble growth mass M_P,grw is large; many pebbles are needed in order to grow the planet. The steepening of the slope that can be noticed at small τ_s is the result of a transition from the shear regime (velocities dominated by the Keplerian shear) at high τ_s to the headwind regime (velocities dominated by the gas sub-Keplerian motion, ηv_K) at small τ_s.

When α_z > 0 (colored symbols) pebbles are stirred to higher regions, reducing the local density of pebbles in the midplane. This reduces ε with respect to the 2D case. As can be seen in Fig. 4 pebbles are most affected when they are small (τ_s ≪ 1) and when the turbulence is strong (high α_z), which is of course natural. Hence, in the 3D case there is a preferred pebble aerodynamic size where ε(τ_s) peaks, which occurs approximately at the point when the pebble accretion impact parameter equals the pebble scaleheight, i.e., at the transition of the 2D and 3D regimes. For heavier particles ε decreases because of more rapid radial drift, whereas for small particles ε decreases because of a reduced local density. However, when α_z ≳ τ_s the pebble scaleheight has reached that of the gas, h_P ≈ h_gas, and no further reduction is possible. As a result, the curves eventually converge when τ_s becomes very small, as is seen in Fig. 4 in the bottom-left corner.

In Ormel (2017), as well as Paper I, we derived that the pebble accretion efficiency in the 3D limit reads (22)

where A₃ is a numerical constant, h_P the pebble aspect ratio, and f_set – the settling fraction – a modulation factor which becomes less than unity when the settling criteria (slow encounters) is no longer fulfilled. The characteristic velocity beyond which settling encounters disappear is (23)

(Ormel & Klahr 2010; Paper I). Qualitatively, when the approach velocity⁷ Δ v ≪ v_* settling is fully operational (f_set = 1), while for Δ v ≫ v_* settling (and therefore pebble accretion) are no longer viable (f_set = 0). Quantitatively, the settling modulation function is fitted empirically by an exponential function (Ormel & Klahr 2010; Visser & Ormel 2016). In Paper I we adopted (24)

with a_set = 0.5. The subscript “I” indicates that this expression is valid for a laminar disk (Paper I). We will refine it below (Sect. 4.6) accounting for a turbulence velocity field.

The reduction of ε by f_set enters quadratically in 3D because the cross section is two-dimensional. By virtue of the rather high planet mass, f_set nevertheless evaluates to unity for most runs in Fig. 4. In Fig. 4 the dotted lines give Eq. (22), where we took f_set = 1 and evaluated h_peb according to Eq. (15). With A₃ = 0.39, this matches the numerical results well for small τ_s. At high τ_s, Eq. (22) clearly overestimates the pebble accretion efficiency; the settling efficiency is then given by its planar limit.

Only forα_z = 0.1 tend the efficiencies to lie below the expression given by Eq. (22). The reason is that now the pebble velocity becomes dominated by turbulent motions since , which suppresses accretion through settling as encounters become too fast for settling, f_set < 1, because of a high turbulent velocity. We present a model to include for turbulence effect in f_set in Sect. 4.6.

4.2 Pebble accretion in the outer disk

For planets on circular orbits, the pebble accretion efficiency is fully determined by the five dimensionless parameters q_p, h_gas, η, τ_s, and α_z. In the outer disk, the aspect ratio h_gas and (as a consequence) η are usually higher, resulting in much lower efficiencies. This is illustrated in Fig. 5, which shows ε(τ_s, α_z) for the inner disk (left panel) and outer disk (right panel) for a q_p = 3 × 10⁻⁸ planet (0.01 M_⊕ for a solar-type star). For the outer disk run, we have increased h_gas and η by factors of ≈2 and ≈5, respectively. For a standard passively irradiated disk model, this would correspond to an increase by a factor 30 in r_p; e.g., we contrast the situation at 1 au (left) with 30 au (right).

Efficiencies in the outer disk are always lower – in particular, the 3D efficiencies. In addition, the transition to the 2D limit occurs at a longer stopping time. Both effects are caused by the larger gas scaleheight. A further consequence of a larger h_gas is that turbulence velocities become higher compared to the critical threshold v_* (both velocities are lower in the outer disk, but whereas σ_z ∝ c_s ∝ h_gasv_K, v_* ∝ v_K). Consequently, for the 30 au run, settling already fails for α = 10⁻², which manifests itself by the flattening and decrease of the curves.

It has been suggested that pebble accretion is an effective mechanism to grow planets in the outer disk (Ormel & Klahr 2010; Lambrechts & Johansen 2012; Bitsch et al. 2015; Johansen & Lambrechts 2017). Compared to planetesimal accretion the key advantage is that capture radii are large due to the large Hill radii, whereas planetesimals suffer from scattering – a negative feedback to the growth of planets (Kobayashi et al. 2010). Nevertheless, as Fig. 5 illustrates, the efficiency of pebble accretion in the outer disk is lower than in the inner disk. For example, growing planets in strongly turbulent (α_z > 10⁻²) disks may require hundreds, if not thousands, of Earth masses in pebbles – numbers that seem rather large in the light of recent ALMA observations (Ansdell et al. 2017; Miotello et al. 2017). Still, invoking (a combination of) low turbulence (as suggested by HL tau; Pinte et al. 2016), pressure bumps (Pinilla et al. 2012), or massive disks (Bitsch et al. 2018a), there is enough leeway to grow planets through pebble accretion in the outer disks. But from an efficiency perspective, it is more conducive to grow them in the inner disk.

Fig. 4 Pebble accretion efficiency vs dimensionless stopping time for several values of the vertical turbulence strength, parameterized by αz (colors). Open symbols give ϵ obtained from directly integrating the stochastic equation of motion, while crosses give the results from the hybrid algorithm. Crosses are offset by 20% to the right for clarity. Error bars correspond to the Poisson counting error on the number of hits. For αz = 0 the triangles give ε resulting from the local calculations and the black line gives our fit to the 2D limit (ϵ2D), which we obtained in Paper I. The dotted lines gives ϵ3D accounting just for the density correction Eq. (22). The colored solid lines give ε accounting for all 3D/2D-effects (Eq. (40)) as explained in the main text. The right y-axis converts ε into the pebble growth mass – the total amount of pebbles needed to e-fold the mass of the planet – assuming a solar-mass star.

Fig. 5 Efficiency of pebble accretion of a 0.01 M⊕ planet (qp = 3 × 10−8) at 1 au (left panel) and at 30 au (right panel). At 1 au hgas = 0.03 and η = 10−3, while at 30 au we take hgas = 0.07 and η = 5 × 10−3. Right axis gives the planet growth mass. At 30 au turbulence more significantly affects the pebble accretion efficiency.

4.3 Efficiency at H₂O iceline for different stellar mass

There have been a number of studies that argue for the H₂O iceline as a preferential site for the formation of the (first) generation of planetesimals and planetary embryos (Cuzzi & Zahnle 2004; Ros & Johansen 2013; Banzatti et al. 2015; Ida & Guillot 2016; Schoonenberg & Ormel 2017; Dra̧żkowska & Alibert 2017). In many of these works, the abundance of ices increases just outside the iceline by condensation of H₂O vapor that diffused back over the iceline. In addition the surface density increases through a “traffic jam” effect when evaporating ice boulders liberate much smaller grains. For these reasons it is worthwhile to consider the efficiency of pebble accretion at the snowline. However, the snowline locations vary with stellar mass. In disks of lower mass stars the snowline will be much closer in, where the disk aspect ratio is likely to be smaller.

In Fig. 6 the pebble accretion efficiency is plotted at the H₂O iceline, contrasting a solar type star (bottom) with a late M-star (top). The planet mass is 0.1 M_⊕ in both cases. Two effects conspire to render efficiencies much higher for icelines around low-mass stars. First, the aspect ratio and η are lower because the iceline lies further in. Second, a planet(esimal) of the same mass around a lower mass star will have a higher mass ratio (q_p). The pebble accretion cross section is then larger because of the reduced Keplerian shear and headwind velocities. Figure 6b illustrates the effect for a 0.1 M_⊙ M-star for (the same) 0.1 M_⊕ mass planet with h_gas = 0.03 and η = 10⁻³; i.e., q_p is higher by a factor of 10, h_gas lower by a factor of 1.6, and η lower by a factor 3 compared to the solar-type star. Clearly, M-star pebble accretion efficiencies are both higher and are less sensitive to variations inτ_s and α_z. Therefore, (late type) M-stars can efficiently convert their pebble-sized building blocks into planetary systems. These finding confirm our earlier analytical estimates on a pebble formation origin of the TRAPPIST-1 system (Ormel et al. 2017).

Fig. 6 Pebble accretion efficiency for a 0.1 M⊕ planet at thelocation of the H2O iceline. Top panel: iceline around a 0.1 M⋆ mass star (hgas = 0.03, η = 10−3; qp = 3 × 10−6). Bottom panel: iceline of a solar mass star (hgas = 0.05, η = 3 × 10−3, and qp = 3 × 10−7).

4.4 Dichotomy of the solar system

It was argued by Morbidelli et al. (2015) that pebble accretion is more efficient beyond the H₂O iceline because icy pebbles, being less prone to collisional fragmentation, are larger (higher τ_s). Specifically, Morbidelli et al. (2015) adopted τ_s = 10^−1.5 for icy pebbles outside the iceline and τ_s = 10^−2.5 for silicate pebbles (just) interior to the iceline and adopted α_z = 10⁻³. The curve corresponding to these parameters is the purple curve in the bottom panel of Fig. 6. The positive slope indeed indicates that larger (icy) pebbles have a higher ε than the smaller (silicate) pebbles, a situation that generally applies to the 3D limit (Eq. (22)) where the pebble aspect ratio h_P decreases with higher τ_s. Morbidelli et al. (2015) concluded that embryos just outside the H₂O iceline will outcompete inner embryos⁸, arguing that the great dichotomy of the solar system – small Mars next to big Jupiter – is well explained under a pebble accretion scenario.

There is a corollary to this hypothesis, not (explicitly) addressed by Morbidelli et al. (2015). As Fig. 6 shows, the efficiencies of the τ_s = 10^−1.5 pebbles (ε ≈ 0.4%) are very small: >99% of the pebbles drift past the planet to enrich the inner solar system. This means that the formation of Jupiter’s core is associated with hundreds of Earth masses in pebbles drifting into the terrestrial planet region. Evidently, in the scenario outlined by Morbidelli et al. (2015), most pebbles could not have been accreted by the bodies in the inner solar system but must have ended up in the young Sun. This puts a constraint on the structure of the inner disk, i.e., it had to be transparent to pebble drift. No dense planetesimals belts or long-lived pressure maxima (which could trigger formation of super-Earths; Chatterjee & Tan 2014) could have existed in the inner solar system during Jupiter’s formation.

4.5 Inclined planets

When planet(esimal)s move on orbits inclined with respect to the disk, the pebble accretion efficiency can be significantly reduced (Johansen et al. 2015; Levison et al. 2015). Planets of high inclination i_p, such that i_p > h_P, only interact with pebbles overa fraction ~i_p∕h_P of their orbits. Inclined planets, therefore, have a similar effect on the settling efficiencies as turbulence: the higher the inclination, the less it interacts with the pebbles. In addition, inclined planets encounter pebbles at an additional velocity (~ i_pv_K), which suppresses accretion once it exceeds v^*.

These effects are illustrated in Fig. 7, where ε is plotted as function of inclination i_p (curves) for the laminar disk (α_z = 0; all pebbles in the midplane), α_z = 10⁻⁴, and α_z = 10⁻². In general, the higher the inclination, the lower the accretion efficiencies. For α_z = 0 inclination effects already become visible for i_p~10⁻³, whereas for α_z = 10⁻² the curves only diverge for i_p > 10⁻² (inclinations are given in radians). Accretion of large τ_s particles in particular are suppressed because of the increase in the approach velocity and the decrease in f_set. The planet moves too fast through the pebble plane.

Since the planet’s inclination has a similar effect on ε as turbulent stirring of pebbles – both reduce the amount of interaction between the two-components – an effective scaleheight may be defined as h_eff ≃ i_p + h_P. A more precise estimate is obtained by averaging the pebble density over the phase of the planet’s orbit (25)

where is the normal distribution with standard deviation h_P and t the phase (mean anomaly). In the limit of i_p ≪ h_P, and h_eff is equal to the pebble scaleheight. More generally, the formal solution to Eq. (25) reads (26)

where I₀(x) is the modified Bessel function of the first kind. For practical purposes, Eq. (26) may be approximated as (27)

In the limit i_p ≫ h_P, , implying that our first guess (h_eff = i_p) is off by 25%. The reason is that the inclined planet spends most of its time near its end points, whereas it moves quickly through the midplane regions.

Fig. 7 Planet inclination dependence on pebble accretion efficiency. We plot ε for our standard parameters (qp = 3 × 10−7, η = 10−3 and hgas = 0.03) as function of τs (x-axis), αz (panels) and planet inclination ip (colors). Curves give our fit (Sect. 4.7). The plotted planet inclination values are ip = 0 (black, top), 10−3, 2 × 10−3, 5 × 10−3, 10−2, 2 × 10−2, 5 × 10−2 and 0.1 radians. Many of the low ip points and curves overlap.

4.6 Role of turbulence – refinement of f_set

We refine the expression for f_set (Eq. (24)) accounting for a, possibly anisotropic, turbulent velocity field. Let Δ v be the non-turbulent component of the approach velocity with its component, Δ v_i, pointing in the, ith direction. Concerning the planar direction, we already obtained Δv in Paper I, which we now relabel Δv_y: (28)

where (29)

is the approach velocity in the circular limit. It combines the headwind (approach velocity dominated by ηv_K) and the shear (approach velocity determined by the Keplerian shear velocity) regimes. In addition, (30)

is the eccentric velocity. In the above formulae, a_cir, a_e and a_sh are all fit constants (see Table 1).

Similarly, for the vertical approach velocity, we have (31)

form the planet’s inclination. The order of unity prefactor a_i is again obtained numerically.

We assume a tri-axial Gaussian velocity distribution of width σ_P, centered on Δv. Explicitly, for component i the approach velocity v_i is normally distributed: (32)

where σ_P,i is the turbulent rms velocity in direction i. Following Youdin & Lithwick (2007), we take (33)

for the pebble rms velocity where ξ is defined in Eq. (16). From Eq. (24), we can write for the accretion probability (34)

The new, distribution-averaged f_set is then obtained, by integration over the velocity distribution: (35)

Formally, the integration gives a_turb = 2a_set, but we relax this constant in order to obtain the best fit to the simulated data.

Equation (35) features the following limits:

• 1.

  σ_P,x, σ_P,y, σ_P,z ≪ v_*. Turbulence is unimportant and the laminar form of f_set is retrieved, Eq. (24).

• 2.

  Isotropic turbulence, σ_P,x = σ_P,y = σ_P,z = σ_P. This simplifies Eq. (35) to

  (36)

• 3.

  σ_P,x = σ_P,y = 0 with turbulence only operating in the vertical direction, which is (by construction) the case considered in this paper and may be applicable to the vertical shear instability (Stoll et al. 2017). Hence, we have used here

  (37)

• 4.

  σ_P,x ≃ σ_P,y ≃ σ_P,z ≃ σ_P ≫ v_* + Δv, the turbulence-dominant limit. In this case

  (38)

  Turbulence reduces the accretion, but not exponentially. In a turbulence-dominated velocity field there is always a fraction of particles with velocities low enough to accrete by settling.

4.7 Summary of the pebble accretion efficiency fit

We provide an executive summary on how to generally obtain ε. Quantities involving the recipe and corresponding numerical constants are given in Table 1.

The settling efficiencies in the 2D and 3D limits read (39a) (39b)

Apart from the f_set term, the 3D expression is independent of the pebble-planet relative velocity; a characteristic feature of pebble accretion in the 3D limit (Ormel 2017; Paper I). The effective scaleheight, h_eff, Eq. (27), accounts for the reduced interaction between pebble and planetesimal by either turbulent stirring of pebbles or planetesimal inclination. The non-turbulent component of the relative motion Δ v is given by Eqs. (28) and (31), for the azimuthal and vertical component respectively.

The other velocity-dependent term is the settling fraction f_set, which becomes important when either laminar or turbulent velocities exceed the critical settling velocity . In the laminar case (σ = 0) f_set is given by Eq. (24), while in the turbulent case (σ≠0) it is given by Eq. (35). When v_* significantly exceeds σ + Δv the settling fraction evaluates to unity.

Finally, we find that the 3D (this Paper) and 2D (Paper I) can be combined as (40)

which ensures a smooth transition. All curves shown in Figs. 4–9 follow this recipe.

When f_set ≪ 1 (e.g., Δv ≫ v_* or σ ≫ v_*) the ballistic regime (see Paper I) takes over. In Paper I we found that the combined efficiency is well-fitted by (41)

where ε_bal is the efficiency in the ballistic limit (see Paper I). However, in 3D growth by ballistic encounters is generally slow.

4.8 Neglected effects

We list several neglected effects, which may change ε:

• 1.

  Aerodynamical deflection. For planetesimals, we have not accounted for the change in v_gas in the vicinity of the planetesimals and correspondingly ignored aerodynamic deflection (Visser & Ormel 2016). This reduction, however, only becomes important when particles are very small, t_stop < R∕v_hw, and only when the encounters operate in the ballistic regime. Under these conditions, turbulence may in fact overcome the aerodynamic barrier (Homann et al. 2016).

• 2.

  Pre-planetary atmospheres. We have not accounted for the effects of the early primordial atmosphere that forms around massive bodies in gaseous disks. As a rule of thumb, this affects the density and flow pattern out to a Bondi radius,. Planets as massive as the thermal mass (h³M_⋆) will have such extensive envelopes that our constant-density assumption will break down. A complex flow structure may prevent small particles from accreting to the planet (Ormel 2013), perhaps after their evaporation (Chambers 2017; Brouwers et al. 2018).

• 3.

  Pressure maxima and resonances. Such massive planets also affect the pressure profile of the disk, resulting in a pressure maximum at a distance ~H_gas from the planet. For this pebble isolation mass, accretion will terminate (Lambrechts et al. 2014; Bitsch et al. 2018b). Similarly, τ_s > 1 pebbles can be stopped at resonant location (Weidenschilling & Davis 1985; Picogna et al. 2018).

• 4.

  Gas radial flow. In viscous disks, gas flows radial at a velocity ~− ν∕r and adds to the drift velocity of pebbles. This effect hence starts to dominate radial drift motions for α_ν > τ_s. It is straightforward to adjust expressions for ε (see, e.g., Ida et al. 2016). Similarly, the drift velocity depends on vertical position, v_r(z), because τ_s increases towards the disk surface. We can correct ε accordingly, i.e., by using a vertically averaged v_r (Takeuchi & Lin 2002; Kanagawa et al. 2017).

• 5.

  By design, the simulations of this work only considered turbulence in the vertical direction – a simplification that allowed us to carry out a numerical parameter study in a controlled way. By modeling vertical turbulence we have accounted for the reduction of the midplane pebble density. On the other hand, the reduction of ε by fast encounters is underestimated in the numerical simulations, particularly when the turbulence is strongly anisotropic (σ_x, σ_y ≫ σ_z). Nevertheless, in Sect. 4.6 we formulated a general recipe for f_set in a anisotropic velocity field.

5 Discussion

Table 2 compiles expressions for ε that have been used in, or derived from, recent studies. For simplicity, we ignore the f_set term and list the numerical prefactor belonging to the leading term for the 3D and the 2D limits. In the 2D limit the shear-dominated velocity regime has been adopted (valid for large planets), while in the 3D limit the expression is independent of Δ v. We emphasize that this work gives the correct numerical prefactor, as it is calibrated against numerical simulations. On the other hand, the existing literature expressions often employed scaling arguments, where typically the 3D rate is estimated to be a fraction b_set ∕H_P of the 2D rate where b_set is the pebble accretion impact parameter. It is therefore quite remarkable that the existing literature prefactors lie so close to our calculated values. Nevertheless, in the 2D-limit the literature expression turn out to be too high, up to 50%. This may still be significant, since a factor of two difference means that planet formation by pebble accretion takes twice as long and requires twice the number of pebbles.

Like our study, Chambers (2014) considered the cases where planet eccentricity (inclinations) dominate and accounts for the velocity effect of turbulence. However, he only considered the turbulent rms velocity for the relative motion between planet and pebble, which suppresses pebble accretion exponentially (f_set ≫ 1) once σ ≫ v_*. This is illustrated in Fig. 8 by the gray curve, where the settling fraction is plotted as function of planet mass. Parameters are chosen such that turbulent rms velocities are similar to the laminar headwind velocity, ηv_K. When f_set is calculated by adding in quadrature the laminar and turbulent rms velocities (as in Eq. (24)), it ensures exponential behavior at low q_p. However, accounting for a velocity distribution changes the functional behavior of f_set and, for the adopted parameters, we obtain the counter intuitive result that turbulence increases f_set (at low q_p; solid curve). Accounting for the distribution, there are always a few particles slow enough for the settling mechanism to operate. This finding is important especially in the outer disk, where even at f_set ~10⁻³ settling interactions dominate over ballistic interactions.

Recently, Xu et al. (2017) measured pebble accretion rates from laminar and MRI-turbulent hydrodynamical simulations. They expressed their result in terms of a dimensionless quantity k_abs, which is the ratio of the mass accretion rate (Ṁ) normalized to the Hill accretion rate . In these units, our expressions for the 2D and 3D efficiencies (Eq. (39)) read⁹: (42a) (42b)

where we did insert Δv in Eq. (42a), but have for clarity omitted the f_set modulation factor.

As a note in passing, Eq. (42a) depends, apart from τ_s, only on the quantity q_p∕η³. The reason is that the pebble equation of motion, expressed in Hill units, only contains this parameter¹⁰. Likewise, the 3D-rates also depend on a single, but different, parameter (assuming h_eff can be quantified in terms of the stopping time). In the general case, then, two parameters (or three if τ_s is included) are necessary to calculate the Hill accretion rate. In their local shearing box simulations, Xu et al. (2017) choose to fix the thermal mass and the quantity η∕h_gas. For the Hill accretion rate (k_abs), the problem is then fully specified. However, to calculate the efficiencies (ε) the degeneracy that existed among q_p, η, and h_gas is broken – q_p, η, and h_gas each need to be specified. This is because the accretion rate is a local quantity, whereas in order to calculate ε the radial pebble flux must be known, i.e., the disk circumference should be specified.

In Fig. 9 symbols correspond to the runs conducted by Xu et al. (2017, their Fig. 4a) for a thermal mass of and a headwind parameter of η = 0.1h_gas. Blue symbols correspond to the hydrodynamic runs (non-turbulent), green symbols to the ambipolar diffusion runs, and red symbols to the ideal MRI runs. As expected, the accretion rate decreases in the turbulent case and more so in the ideal MRI run than in the AD-restive run. Solid curves give the corresponding k_abs obtained from Eq. (42), f_set, and the 2D and 3D averaging formula (Eq. (40))¹¹. Fits are plotted until τ_s = 1 beyond which they loose their validity.

The results of our analytical model agree well with Xu et al. (2017). For the turbulent runs, we find that the reduction of the accretion rate is mostly due to turbulent diffusion (turbulence lofting pebbles away from the midplane) rather than the turbulent velocity effect. This explains why the ideal MRI run features lower accretion rates than the ambipolar diffusion runs. Towards τ_s = 0.1–1 we see that the ambipolar diffusion run and the hydrodynamic run converge, as the accretion becomes 2D and turbulent diffusion does not enter k_abs. However, we also find that the ideal MRI run does not entirely converge on the hydrodynamic run. This can be attributed to the turbulent velocity effect; f_set ≲ 1 because σ ≳ v_* for the τ_s ≃ 0.1–1 particles in the MRI run. In other words, settling is marginally failing. We anticipate a much stronger reduction for either smaller planets or more vigorous turbulence.

Fig. 8 Settling fraction fset as function of planet mass for αz = 10−2, hgas = 0.05, η = 5 × 10−3, τs = 0.1 and M⋆ = 1 M⊙. Turbulence is assumed to be isotropic and the correlation time tcorr = Ω−1. Curves give fset when it is calculated based on the laminar-only velocity (dashed), turbulent rms velocity (gray), and for a distribution of velocities (solid).

Fig. 9 Hill-normalized accretion rates for a planet of thermal mass and η∕hgas = 0.1. Data points give kabs obtained from hydrodynamical laminar (2D), ideal MRI (ID), and resistive (ambipolar diffusion) MRI (AD) simulations reported by Xu et al. (2017). Solid lines give our fits, where we adopt σ from the same simulations. Data kindly provided by Ziyan Xu.

6 Summary

In this work we have developed a general framework for the stochastic equation of motion (SEOM) of pebble-sizedparticles. The SEOM is described by the particle’s stopping time, by the gas diffusivity D_gas and correlation time t_corr, and by gravitational forces. Using the SEOM we have investigated the vertical transport of particles in disks.

• 1.

  From the SEOM we obtain the strong coupling approximation (SCA) by taking the limit of t_stop → 0, as was used in previous studies (Ciesla 2010; Zsom et al. 2011). For small t_stop and t_corr the SEOM becomes consistent with the SCA.

• 2.

  Exploring the effect of a vertically dependent diffusivity, α_z(z), such that the midplane diffusivity α_mid ≪ τ_s and the surface diffusivity α_surface ≫ τ_s, we find a two component distribution for the pebbles. Although most of the pebbles are concentrated in the midplane (following the scaleheight given by α_mid), a small fraction of pebbles are distributed according to the gas scaleheight. Then, high-velocity collisions among these pebbles could explain the prolonged presence of small, micron-sized grains in the disk surface.

As its main application, we have integrated pebble trajectories to find the pebble accretion efficiency ε.

• 1.

  Compared to the laminar disk, turbulence reduces ε in two ways: it diminishes the local density of pebbles at the midplane through diffusion and it increases the rms velocity of particles. The latter becomes important for low planet masses, where fulfilling the settling condition becomes more difficult.

• 2.

  Because of the disk geometry, pebble accretion is more efficient in the inner disk. From an efficiency perspective, accretion around low-mass stars is also favored because pebble capture radii are larger around low-mass stars.

• 3.

  Together with the 2D expressions already derived in Paper I, we have formulated a general prescription for ε as a function of pebble properties (aerodynamical size), planet properties (mass, inclination, eccentricity), and disk properties (pressure profile, gas density, turbulence properties). This prescription is summarized in Sect. 4.7.

Finally, we remark that ε has been defined with respect to a single planet. A small ε therefore does not necessarily imply that pebble accretion is globally inefficient. Indeed, in the case of a planetesimal belt (or multiple planets), the total filtering efficiency may well reach unity, whereas the individual ε ≪1 (Guillot et al. 2014). This raises the question of how pebble accretion proceeds when multiple seeds are present.

While it is likely to result in a more efficient pebble sweep-up, a multi-seed scenario could also suppress planet growth because of the mutual dynamical excitation among the protoplanets (Levison et al. 2015). In a following work, we will investigate the efficacy of pebble accretion when pebbles interact with a narrow planetesimal belt (Liu et al., in prep.). In such cases, the combined planetesimal and pebble coagulation can best be studied with N-body techniques, where the pebble accretion rate on the N-bodies is obtained using the prescription for ε summarized in Sect. 4.7.

Acknowledgements

The authors are supported by the Netherlands Organization for Scientific Research (NWO; VIDI project 639.042.422). We thank Xuening Bai and Ziyan Xu for sharing and discussing their simulations, and Ramon Brasser and Sebastiaan Krijt for proofreading the manuscript. The comments from the referee improved the clarity of the manuscript.

Appendix A Proof of Equation (10)

In deriving Eq. (10), we follow the proof outlined by Hottovy et al. (2015). The first step is to write Eq. (4) as (A.1a) (A.1b)

In this equation, the stochastic variable ζ_t has been combined with the velocity into one vector v = (v, ζ_t). Similarly, x contains an additional parameter, say υ, but this is entirely dummy. Comparing with Eq. (4), we therefore have (A.2a) (A.2b) (A.2c)

It must be emphasized that we only treat here the z coordinate. When the full 3D equation of motion is considered, v will be vector of length six and γ a matrix of 36 elements, containing all entries of the diffusion tensor D_ij.

In the limit of t_stop → 0, Hottovy et al. (2015) proves that x can be described by the SDE (A.3)

where γ⁻¹ is the inverse of γ and S is the noise-induced drift term – a vector whose ith component is defined (A.4)

with J the solution of the Lyapunov equation (A.5)

Since we consider here only the z coordinate, Eq. (A.3) reduces too (A.6)

The equations can now be solved. Inverting γ we find (A.7)

Furthermore, solving the Lyapunov equation, we find (A.8)

with which the noise-induced drift term becomes (A.9)

and we retrieve Eq. (10).

A.1 Stratonovich interpretation

A feature peculiar to stochastic integrals is that equations of the form (A.10)

are ill-defined when B is a function of position. In contrast to ODEs, it does matter whether the integrand is evaluated at t (i.e., B(x) = B(x[t]) – Ito’s choice), at t + Δt, or whether we let (A.11)

(Stratonovich’ choice). These definitions will produce different results when B is not constant (van Kampen 1992). Specifically, with Stratonovich’, rather than Ito’s, interpretation for the stochastic integral, a term (A.12)

should be added to the RHS of the Fokker–Planck Eq. (12).

In deriving Eq. (10), as well as the conversion from Eq. (10) to Eq. (13), we have followed Ito’s interpretation. On the other hand, under Stratonovich interpretation, the RHS of Eq. (10) and the RHS of Eq. (13) gain an additional term : (A.13)

Hence, the SDE for the strong coupling approximation depends on the interpretation rule (a fact not highlighted by Ciesla 2010 or Zsom et al. 2011). The underlying reason between the Ito and Stratonovich interpretations reflects the nature of the stochastic forcing (van Kampen 1992). Ito’s interpretation would hold, for example, when the stochastic forcing amounts to a series of infinitely short “pulses” with every pulse completely independent. This interpretation is often used in finance. However, it may be argued that for physical problems, where the correlation time is never really infinitely small, Stratonovich amounts to the more correct model (van Kampen 1992).

More practically, these differences are expressed in our choice of the numerical integration scheme. When Ito’s interpretation is adopted, the stochastic equation must be integrated with a corresponding numerical scheme, of which the Euler method is the simplest example. Similarly, the Stratonovich equation should be integrated with an appropriate numerical scheme, where the midpoint scheme (Heun’s method) is the simplest example. In our code, where we use an Runge–Kutta method, which is a generalization of a midpoint scheme, we therefore adopt Eq. (A.13), instead of Eq. (10).

References

All Tables

All Figures

Fig. 1 Top panel: normalized distributions of the vertical position z for the integration of the tracer case, using the strong coupling approximation (SCA) with tstop = 0. The vertical height z is recorded after every 1 Ω−1 for t = 105 Ω−1. Bars give the simulated distribution while the analytic – normal – distribution is shown by the black dashed line. Bottom panel: temporal evolution of the vertical position for αz = 10−2 (black) and αz = 10−2 (blue).

In the text

Fig. 2 Vertical distribution for particles of different stopping times: τs = tstopΩ = 10−2 (panel a), 10−1 (panel b), and 1 (panel c). Histograms plot the numerically obtained distributions with the stochastic equation of motion (SEOM; gray and red) and strong coupling approximation (SCA; blue) methods. Long tcorr runs are shown with red histograms (the tcorr = 102 Ω−1, τs = 0.1 run is displayed in panel c). Thin curves gives the normal distribution with the scaleheight of Eq. (15) (Youdin & Lithwick 2007). Note the different scaling among the panels.

In the text

Fig. 3 Density distribution P(z) obtainedfrom a vertically varying diffusivity. The diffusivity profile in terms of αz is given by the red curve. The particle stopping time is τs = 10−2 and is taken independent of z. The distributions obtained from integrating the stochastic equation of motions (Eq. (4); black histogram) and the strong coupling approximation (Eq. (10); gray histograms) are consistent. The dark blue curve gives the normal distribution, for the particle scaleheight evaluated in the midplane (i.e., Eq. (15) with α = 10−3). The light blue curves gives the gas distribution, scaled by a factor 0.1.

In the text

Fig. 4 Pebble accretion efficiency vs dimensionless stopping time for several values of the vertical turbulence strength, parameterized by αz (colors). Open symbols give ϵ obtained from directly integrating the stochastic equation of motion, while crosses give the results from the hybrid algorithm. Crosses are offset by 20% to the right for clarity. Error bars correspond to the Poisson counting error on the number of hits. For αz = 0 the triangles give ε resulting from the local calculations and the black line gives our fit to the 2D limit (ϵ2D), which we obtained in Paper I. The dotted lines gives ϵ3D accounting just for the density correction Eq. (22). The colored solid lines give ε accounting for all 3D/2D-effects (Eq. (40)) as explained in the main text. The right y-axis converts ε into the pebble growth mass – the total amount of pebbles needed to e-fold the mass of the planet – assuming a solar-mass star.

In the text

	Fig. 5 Efficiency of pebble accretion of a 0.01 M⊕ planet (qp = 3 × 10−8) at 1 au (left panel) and at 30 au (right panel). At 1 au hgas = 0.03 and η = 10−3, while at 30 au we take hgas = 0.07 and η = 5 × 10−3. Right axis gives the planet growth mass. At 30 au turbulence more significantly affects the pebble accretion efficiency.
In the text

	Fig. 6 Pebble accretion efficiency for a 0.1 M⊕ planet at thelocation of the H2O iceline. Top panel: iceline around a 0.1 M⋆ mass star (hgas = 0.03, η = 10−3; qp = 3 × 10−6). Bottom panel: iceline of a solar mass star (hgas = 0.05, η = 3 × 10−3, and qp = 3 × 10−7).
In the text

Fig. 7 Planet inclination dependence on pebble accretion efficiency. We plot ε for our standard parameters (qp = 3 × 10−7, η = 10−3 and hgas = 0.03) as function of τs (x-axis), αz (panels) and planet inclination ip (colors). Curves give our fit (Sect. 4.7). The plotted planet inclination values are ip = 0 (black, top), 10−3, 2 × 10−3, 5 × 10−3, 10−2, 2 × 10−2, 5 × 10−2 and 0.1 radians. Many of the low ip points and curves overlap.

In the text

	Fig. 8 Settling fraction fset as function of planet mass for αz = 10−2, hgas = 0.05, η = 5 × 10−3, τs = 0.1 and M⋆ = 1 M⊙. Turbulence is assumed to be isotropic and the correlation time tcorr = Ω−1. Curves give fset when it is calculated based on the laminar-only velocity (dashed), turbulent rms velocity (gray), and for a distribution of velocities (solid).
In the text

	Fig. 9 Hill-normalized accretion rates for a planet of thermal mass and η∕hgas = 0.1. Data points give kabs obtained from hydrodynamical laminar (2D), ideal MRI (ID), and resistive (ambipolar diffusion) MRI (AD) simulations reported by Xu et al. (2017). Solid lines give our fits, where we adopt σ from the same simulations. Data kindly provided by Ziyan Xu.
In the text