Pebble Accretion in Turbulent Protoplanetary Disks

We use the Athena magnetohydrodynamics (MHD) code (Stone et al. 2008), which is a higher-order Godunov code with constrained transport to enforce the divergence-free constraint on the magnetic field, to perform local 3D MHD simulations under the shearing-sheet approximation (Goldreich & Lynden-Bell 1965). We use a Cartesian coordinate system in a corotating frame located at a fiducial radius with Keplerian frequency ${{\rm{\Omega }}}_{K}$. The dynamical equations are then written with $\hat{{\boldsymbol{x}}},\hat{{\boldsymbol{y}}},\hat{{\boldsymbol{z}}}$ denoting unit vectors pointing in the radial, azimuthal, and vertical directions respectively, where ${{\boldsymbol{\Omega }}}_{K}$ is along the $\hat{{\boldsymbol{z}}}$ direction. With gas density, gas velocity, and magnetic field denoted by ${\rho }_{g},{\boldsymbol{u}}$, and ${\boldsymbol{B}}$, the MHD equations can be written as follows in this non-inertial frame:

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\rho }_{g}}{\partial t}+{\rm{\nabla }}\cdot ({\rho }_{g}{\boldsymbol{u}})=0,\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\rho }_{g}{\boldsymbol{u}}}{\partial t}+{\rm{\nabla }}\cdot ({\rho }_{g}{{\boldsymbol{u}}}^{T}{\boldsymbol{u}}+{\mathsf{T}})\\ &&\quad =\,{\rho }_{g}[2{\boldsymbol{u}}\times {{\boldsymbol{\Omega }}}_{K}+3{{\rm{\Omega }}}_{K}^{2}x\hat{{\boldsymbol{x}}}+2{\rm{\Delta }}{v}_{K}{{\rm{\Omega }}}_{K}\hat{{\boldsymbol{x}}}-{\rm{\nabla }}{{\rm{\Phi }}}_{P}],\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\boldsymbol{B}}}{\partial t}={\rm{\nabla }}\times \left[{\boldsymbol{u}}\times {\boldsymbol{B}}+\displaystyle \frac{({\boldsymbol{J}}\times {\boldsymbol{B}})\times {\boldsymbol{B}}}{c\gamma {\rho }_{i}{\rho }_{g}}\right],\end{eqnarray} \tag{ 3 }$$

where ${\mathsf{T}}$ is the total stress tensor

$$\begin{eqnarray}&&{\mathsf{T}}=(P+{B}^{2}/8\pi )\ {\mathsf{I}}-\displaystyle \frac{{{\boldsymbol{B}}}^{T}{\boldsymbol{B}}}{4\pi },\end{eqnarray} \tag{ 4 }$$

${\mathsf{I}}$ is the identity tensor, P is the gas pressure, and ${{\rm{\Phi }}}_{P}$ is the gravitational potential of the planetary core (see Section 3.2). We assume an isothermal equation of state $P={\rho }_{g}{c}_{s}^{2}$, where c_s is the isothermal sound speed. The disk scale height is given by $H={c}_{s}/{{\rm{\Omega }}}_{K}$. In code units, we set ${c}_{s}={{\rm{\Omega }}}_{K}=H=1$. Disk vertical gravity in the gas is ignored and all of our simulations are vertically unstratified. This is because the length scale relevant to pebble accretion (the Hill radius of the planetary core) that we investigate in this paper is much smaller than H (see Section 3.2). The initial gas density ${\rho }_{0}$ is hence uniform, and we set ${\rho }_{0}=1$ in code units. Periodic boundary conditions are used in the azimuthal and vertical directions, while the radial boundary conditions are shearing periodic (Hawley et al. 1995).

The first and second terms on the right side of (2) are standard shearing-sheet source terms of Coriolis force and tidal gravity, leading to velocity shear along $\hat{x}$. In practice, these terms are modified to subtract background shear motion using the orbital advection algorithm of Stone & Gardiner (2010) that improves the accuracy.

The third term corresponds to a constant force pointing radially outward (which we have newly implemented), representing a global radial pressure gradient in the disk. As a result, the disk rotates more slowly than Keplerian by ${\rm{\Delta }}{v}_{K}$. This sub-Keplerian rotation is the source of the radial drift of dust particles, and sets the velocity scale of dust particles approaching the planetary core (which rotates at Keplerian velocity and is stationary in our simulation frame).⁷ It is convenient to normalize ${\rm{\Delta }}{v}_{K}$ by the sound speed c_s. We are mainly interested in the outer regions of PPDs, where pebble accretion is expected to play a dominant role in planetary core growth. In our simulations, we fix ${\rm{\Delta }}{v}_{K}=0.1{c}_{s}$, which is generally applicable in the outer PPDs.⁸

The last term in (3) represents AD, with γ denoting the coefficient of momentum exchange in ion–neutral collisions, and ${\rho }_{i}$ being the ion density. In disks, AD can be most conveniently parameterized by the Elsässer number, defined as

$$\begin{eqnarray}&&{Am}\equiv \displaystyle \frac{\gamma {\rho }_{i}}{{{\rm{\Omega }}}_{K}}.\end{eqnarray} \tag{ 6 }$$

The value of Am marks the importance of AD: the ideal MHD regime corresponds to ${Am}\to \infty$, where the magnetic field is frozen into the gas, while AD significantly affects gas dynamics when ${Am}\lesssim 10$ (Bai & Stone 2011). In PPDs, AD is the dominant non-ideal MHD effect in the outer region ($R\gtrsim 30$ au), where ${Am}\approx 1$ is found to be widely applicable (Bai 2011).

We impose a net vertical magnetic field B₀ in our simulations. The strength of this net vertical field is measured by the plasma β:

$$\begin{eqnarray}&&{\beta }_{0}\equiv \displaystyle \frac{{\rho }_{0}{c}_{s}^{2}}{{B}_{0}^{2}/8\pi },\end{eqnarray} \tag{ 7 }$$

which is the ratio of gas pressure to the magnetic pressure of the net vertical field. Realistic simulations of PPD gas dynamics in the outer disk suggested that net vertical magnetic field is needed for the MRI turbulence to be sustained in the presence of strong AD (Bai & Stone 2011; Simon et al. 2013). It has also been found that ${\beta }_{0}\sim {10}^{4}$ is needed to achieve the desired accretion rates consistent with observations (Simon et al. 2013; Bai 2015). In all MHD simulations described in this paper, we set ${\beta }_{0}=1.28\times {10}^{4}$.

We perform simulations with three different turbulence levels. First, we perform pure hydrodynamic simulations of pebble accretion for reference. These runs are labeled "hyd" in Table 1. We then perform non-ideal MHD simulations with Am = 1, which represent a realistic level of turbulence in the midplane regions for the outer PPDs, labeled by "AD" in Table 1. Finally, we perform simulations in ideal MHD, corresponding to an exaggerated level of turbulence, labeled by "idl" in the Table. These simulations are described in detail in Section 3.2.

Table 1.  Characteristic Parameters of the Simulations in this Work

Simulation	Resolution	Box Size (H)	${\text{}}{M}_{c}/{\text{}}{M}_{T}$	${r}_{H}(H)$	Hp(H)	${\beta }_{0}$	Am
hyd3	192/H	1 × 4 × 1	3 × 10−3	0.1	0.01	$\infty$	⋯
AD3	192/H	2 × 4 × 1	3 × 10−3	0.1	S.C.	12,800	1
idl3	192/H	3 × 6 × 1	3 × 10−3	0.1	S.C.	12,800	$\infty$
hyd2	96/H	2 × 4 × 1	3 × 10−2	0.23	0.01	$\infty$	⋯
AD2	96/H	2 × 4 × 1	3 × 10−2	0.23	S.C.	12,800	1
idl2	96/H	3 × 4 × 1	3 × 10−2	0.23	S.C.	12,800	$\infty$

Note. All simulations have a particle number of 24,576 of each type per H³, and assume ${\rm{\Delta }}{v}_{K}=0.1{c}_{s}$. "S.C." stands for self-consistent.

Download table as:  ASCIITypeset image

Dust particles are included in our simulations using the particle module described in Bai & Stone (2010). The dust particles are treated as test particles that passively respond to the gas flow without exerting a backreaction on the gas, which allows us to cleanly separate the effect of the MRI turbulence from further complications due to particle feedback. We also note that Lambrechts & Johansen (2012) found that pebble accretion rate are essentially unaffected when the particles' backreaction is included.

In the shearing-sheet framework, the equation of motion for particle i reads

$$\begin{eqnarray}&&\displaystyle \frac{d{{\boldsymbol{v}}}_{i}}{{dt}}=2{{\boldsymbol{v}}}_{i}\times {{\boldsymbol{\Omega }}}_{K}+3{{\rm{\Omega }}}_{K}^{2}{x}_{i}\hat{{\boldsymbol{x}}}-{{\rm{\Omega }}}_{K}^{2}{z}_{i}\hat{{\boldsymbol{z}}}-\displaystyle \frac{{{\boldsymbol{v}}}_{i}-{\boldsymbol{u}}}{{t}_{{\rm{s}}}}+{{\boldsymbol{a}}}_{P},\end{eqnarray} \tag{ 8 }$$

where v_i is particle velocity, t_s is the particle stopping time characterizing the drag force, and ${{\boldsymbol{a}}}_{P}=-{\rm{\nabla }}{{\rm{\Phi }}}_{P}$ is the acceleration due to the planetary core (see Section 3.2). Note that unlike the gas, particles are not affected by the radial pressure gradient in the gas. The resulting difference between gas and particle velocities causes the particles to feel gas drag, which leads to their radial drift. In addition, we have included vertical gravity for the particles, allowing them to settle toward the disk midplane (see Section 3.1).

A particle's stopping time t_s depends on properties of both the gas and the dust. In the low-density outer region of the disk, the relevant drag law is the Epstein law and ${t}_{s}={\rho }_{s}a/{\rho }_{g}{c}_{s}$ (Epstein 1924), where ${\rho }_{s}$ is the density of the dust material and a is the particle radius, which is smaller than the mean free path of surrounding gas. It is most convenient to use a dimensionless stopping time ${\tau }_{s}\equiv {{\rm{\Omega }}}_{K}{t}_{s}$. Particles with ${\tau }_{s}\ll 1$ are strongly coupled to the gas, and particles with ${\tau }_{s}\gg 1$ are loosely coupled to it. Note that in the standard MMSN model, centimeter-sized particles have a dimensionless stopping time of ${\tau }_{s}\sim 0.1$ at 10 au in the midplane, and ${\tau }_{s}\sim 1$ at ∼40 au (see Chiang & Youdin 2010 for a review).

In our pebble accretion simulation (see Section 3.2), particles of seven different sizes are considered with stopping times ranging from ${\tau }_{s}={10}^{-2}$ to ${\tau }_{s}=10$. In our simulations, we fix ${\tau }_{s}$ to be constant for individual particles. This is because gas density is largely constant in the local disk midplane regions of interest, as in our unstratified simulations.

Because of vertical gravity, dust particles settle toward the disk midplane, balanced by turbulent diffusion. Therefore, before proceeding to pebble accretion simulations, we first conduct MRI simulations without a planetary core to determine the vertical profile of dust particles under different levels of turbulence in a steady state. This will be used to initialize pebble accretion simulations to be described in the next subsection.

Two simulations are conducted, in either ideal MHD or non-ideal MHD with Am = 1. Both runs use the same box size of $H\times 4H\times H$ in $x,y$, and z dimensions, resolved by $192\times 384\times 192$ cells. The use of such a high resolution first guarantees that the most unstable MRI wavelength ${\lambda }_{m}$ is resolved. In ideal MHD, we have ${\lambda }_{m}\approx 9.18{\beta }_{0}^{-1/2}H\approx 0.081H$ for ${\beta }_{0}={\rm{12,800}}$. This corresponds to ∼16 cells per scale height, which is reasonably well resolved. For non-ideal MHD with Am = 1, ${\lambda }_{m}$ is roughly doubled (Wardle 1999; Bai & Stone 2011), and we expect the MRI to be very well resolved in this run. In addition, the high resolution we use here is necessary to properly resolve the Hill sphere of the planetary core in our pebble accretion simulations.

We first run simulations without particles to time $t=120{{\rm{\Omega }}}_{K}^{-1}$ when the MRI is expected to fully saturate. We then inject nine types of particles, each characterized by a constant dimensionless stopping time ${\tau }_{s}$ spanning from 0.01 to 100 with two particle types per decade in ${\tau }_{s}$. We inject 28,800 particles per type, whose spatial distribution follows a Gaussian profile $\propto \exp (-{z}^{2}/2{H}_{p0}^{2})$, where H_p0 is an initial guess of particle scale height ranging from 0.003 for the most loosely coupled particles to 0.3 for the most strongly coupled. We further wait for another $\gtrsim 20$ orbits for particles to interact with turbulence before starting to measure particle vertical profiles.

The profiles can be well fitted by a Gaussian, characterized by particle scale height H_p, which is expected from the balance between vertical settling and turbulent diffusion. The value of H_p is related to the strength of turbulence, and particle stopping time τ. Theoretically, we expect (Youdin & Lithwick 2007)

$$\begin{eqnarray}&&{H}_{p}\approx \sqrt{\displaystyle \frac{{D}_{g,z}}{{\rm{\Omega }}{\tau }_{s}}}=\sqrt{\displaystyle \frac{{\alpha }_{z}}{{\tau }_{s}}}H,\end{eqnarray} \tag{ 9 }$$

where ${D}_{g,z}\equiv {\alpha }_{z}{c}_{s}H$ is the turbulent diffusion coefficient in the gas.

The results for the normalized particle scale heights ${\text{}}{H}_{p}/{\text{}}H$ in our simulations as a function of ${\tau }_{s}$ are shown in Figure 1. The dependence of H_p on ${\tau }_{s}$ in both simulations nicely follows a power law with index close to (but not exactly) $-1/2$, in reasonable agreement with theoretical expectations, as well as with previous simulation results (Carballido et al. 2011; Zhu et al. 2015). The particle scale heights obtained here (based on the dashed lines shown in Figure 1) will be used for initializing the pebble accretion simulations. By fitting the results using (9), we find a vertical diffusion coefficient ${\alpha }_{z}=7.8\times {10}^{-4}$ for the AD simulation and ${\alpha }_{z}=4.4\times {10}^{-3}$ for the ideal MHD simulation.

Zoom In Zoom Out Reset image size

For reference, we also measure the α parameter for disk angular momentum transport (Shakura & Sunyaev 1973), which is obtained by evaluating the time- and volume-averaged sum of the Maxwell stress and Reynolds stress, normalized by thermal pressure:

$$\begin{eqnarray}&&\alpha =\left\langle \displaystyle \frac{-{B}_{x}{B}_{y}+\rho {v}_{x}{v}_{y}}{\rho {c}_{s}^{2}}\right\rangle .\end{eqnarray} \tag{ 10 }$$

We find $\alpha =6.8\times {10}^{-4}$ in the AD simulation and $\alpha =0.045$ in the ideal MHD simulation. We note that the difference in ${\alpha }_{z}$ in between the ideal MHD and AD runs is much smaller than their difference in α. This is related to the longer turbulent correlation time in the AD case, as studied in Zhu et al. (2015). More specifically, we find the mean square of the vertical turbulent velocity to be $\langle {v}_{z}^{2}/{c}_{s}^{2}\rangle =3.0\times {10}^{-4}$ and $\langle {v}_{z}^{2}/{c}_{s}^{2}\rangle \,=0.0121$ in the AD and ideal simulations, respectively. This is consistent with the results of Zhu et al. (2015) where ${D}_{g,z}\sim (3\mbox{--}4){{\rm{\Omega }}}^{-1}\langle {v}_{z}^{2}\rangle$ in the AD case and ${D}_{g,z}\lesssim 0.5{{\rm{\Omega }}}^{-1}\langle {v}_{z}^{2}\rangle$ in the ideal MHD case.

We conduct pebble accretion simulations with three different levels of turbulence as mentioned before, and two different planetary core masses M_c, with a total of six runs. Their simulation parameters are listed in Table 1 and described in more detail below.

We choose the planetary core mass to be such that it is sufficiently massive for pebble accretion to proceed in the most efficient "Hill regime" (Lambrechts & Johansen 2012) while not too massive to open a gap in the disk. The core mass is best normalized by the thermal mass (Lin & Papaloizou 1993), defined as

$$\begin{eqnarray}&&{M}_{T}=\displaystyle \frac{{c}_{s}^{3}}{G{{\rm{\Omega }}}_{K}}\approx 160\,{M}_{\oplus }{\left(\displaystyle \frac{R}{30\mathrm{au}}\right)}^{3/4},\end{eqnarray} \tag{ 11 }$$

where the number in the approximate equality is obtained based on the temperature profile of the MMSN disk model around a protostar of solar mass. With this definition, the Hill radius of the planetary core is given by

$$\begin{eqnarray}&&{r}_{H}={\left(\displaystyle \frac{{M}_{c}}{3{M}_{T}}\right)}^{1/3}.\end{eqnarray} \tag{ 12 }$$

Lambrechts & Johansen (2012) defined the "transition mass" for pebble accretion as

$$\begin{eqnarray}&&{M}_{t}=\displaystyle \frac{{\rm{\Delta }}{v}_{K}^{3}}{G{{\rm{\Omega }}}_{K}}\approx 0.160\,{M}_{\oplus }{\left(\displaystyle \frac{R}{30\mathrm{au}}\right)}^{3/4}{\left(\displaystyle \frac{{\rm{\Delta }}{v}_{K}}{0.1{c}_{s}}\right)}^{3},\end{eqnarray} \tag{ 13 }$$

where they showed that beyond this core mass, particles with stopping time ${\tau }_{s}\sim 0.1\mbox{--}1$ passing the core within the Hill radius can be accreted by the core (which defines the "Hill regime"). For our choice of ${\rm{\Delta }}{v}_{K}=0.1{c}_{s}$, we have ${M}_{t}={10}^{-3}\,{M}_{T}$. In our simulations, we choose the core masses to be ${M}_{c}=3\times {10}^{-3}\,{M}_{T}$ and $3\times {10}^{-2}\,{M}_{T}$, and we attach the numbers 3 and 2 to the names of the corresponding runs. These core masses guarantee that pebble accretion occurs in the Hill regime. We avoid choosing an even larger core mass, which would notably affect the bulk flow structure and potentially open gaps in the vicinity of the core (Dong et al. 2011; Ormel 2013), leading to additional effects that make it more difficult to isolate the role of turbulence.

The Hill radii for the two choices of core mass are ${r}_{H}=0.1H$ and ${r}_{H}=0.23H$, respectively. In the former case, we choose the grid resolution to be 192 cells per H, so that r_H is resolved by about 20 cells. In the latter case, we reduce the resolution to 96 cells per H and the Hill radius is resolved by about an equal number of cells.

In the simulations, the planetary core is placed at the center of the simulation box, with its gravitational potential given by

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{p}=-{{GM}}_{c}\displaystyle \frac{{r}^{2}+3{R}_{s}^{2}/2}{{({r}^{2}+{R}_{s}^{2})}^{3/2}},\end{eqnarray} \tag{ 14 }$$

where r is the distance to the core, and R_s is the softening length. The smoothing function is accurate to fourth order at $r\gg {R}_{s}$ while it avoids divergence at $r\lt {R}_{s}$ (e.g., Dong et al. 2011). The softening lengths R_s for the gas are chosen to be about $0.02H$ for most of the simulations, except for runs idl3 and AD3, which use $0.01H$. We use half the gas softening length for particles,⁹ which is only slightly larger than the grid scale and is much smaller than r_H (by a factor ∼10). This ensures that the process of pebble accretion is largely unaffected by the smoothed gravitational potential. For runs idl3 and AD3, using a greater softening length would make some particles in our smallest size bin (${\tau }_{s}=0.01$) be stripped from the artificially shallow potential well after they are accreted.¹⁰

We also note that the Bondi radius of the core ${R}_{B}={{GM}}_{c}/{c}_{s}^{2}=({M}_{c}/{M}_{T})H\ll H$ is roughly at grid resolution for ${M}_{c}=3\times {10}^{-3}\,{M}_{T}$, and it is resolved by about three cells for ${M}_{c}=3\times {10}^{-2}\,{M}_{T}$. While this is not sufficient to properly resolve the gas flow in the vicinity of the planetary core as well as its atmosphere (e.g., Ormel et al. 2015), our choice of particle smoothing length R_s guarantees that particles of the sizes considered here (${\tau }_{s}\gtrsim 0.01$) are strongly bound to the core once accreted.

At the beginning of our pebble accretion simulations, we first gradually introduce the gravity from the planetary core over a period of ∼10 orbits ($60{{\rm{\Omega }}}_{K}^{-1}$) without adding particles. This allows the gas to adapt to the presence of the core. For all MHD runs, we further run the simulations to time $t=180{{\rm{\Omega }}}_{K}^{-1}$ to allow the MRI turbulence to fully develop. We then start to inject particles. We use seven particle types spanning dimensionless stopping time from ${\tau }_{s}=0.01$ to ${\tau }_{s}=10$, and each particle type contains 24,576 particles per H³. For MHD simulations, particles are injected uniformly in the horizontal domain, and their vertical distribution follows the Gaussian profile $\propto \exp (-{z}^{2}/2{H}_{p}^{2})$ with particle scale heights H_p for each particle type derived from previous diffusion simulations. For hydrodynamic simulations, all the particles are initialized to be in the midplane. All particles are initialized with velocities following the solution in the laminar disk given by Equations (15) and (16) below. As in Lambrechts & Johansen (2012), we use open/outflow boundary conditions for particles (not gas) in the radial and azimuthal directions so that any particles that leave the box will no longer re-enter from the other side. Correspondingly, all particles transiently pass the core in the simulations, and either get captured or get lost.

In the MHD runs, because of the stochastic nature of the MRI turbulence, we repeat the particle injection process multiple times while we continue running the MRI simulations, which allows us to improve the statistics by sampling different realizations of the MRI turbulence. Note that the timescale for particle passage through the core is of the order of $H/{\rm{\Delta }}{v}_{K}\sim 10{{\rm{\Omega }}}_{K}^{-1}$. In practice, we repeat particle injection for six cycles for the most MRI simulations, with each cycle spanning ${\rm{\Delta }}{t}_{\mathrm{cycle}}=30{{\rm{\Omega }}}_{K}^{-1}$, except for run idl3, where we cover 10 cycles with each cycle spanning ${\rm{\Delta }}{t}_{\mathrm{cycle}}=24{{\rm{\Omega }}}_{K}^{-1}$.

In our simulations, we adopt box sizes from $H\times 4H\times H$ to $3H\times 6H\times H$ (see Table 1). They are chosen in the main to enable a larger fraction of particles in the simulation box to be accreted, and in particular, we use the relatively extended domain in the radial dimension in order to access particles carried by large turbulent eddies on scales of $\sim H$.

Figure 2 shows snapshots of projected particle densities from initial to later ($T=8.5{{\rm{\Omega }}}_{K}^{-1}$) times in the central part of ideal MHD and AD simulations for particle sizes ${\tau }_{s}=0.1$ and 1. The snapshots correspond to the first cycle of particle injection. In the absence of turbulence, particles drift relative to Keplerian orbits in both radial and azimuthal directions because of sub-Keplerian motion of the gaseous disk. For particles with dimensionless stopping time ${\tau }_{s}$, their drift velocity relative to the local Keplerian velocity is given by

$$\begin{eqnarray}&&{v}_{r}=-2\displaystyle \frac{{\tau }_{s}}{{\tau }_{s}^{2}+1}{\rm{\Delta }}{v}_{K},\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}&&{v}_{\phi }=-\displaystyle \frac{1}{{\tau }_{s}^{2}+1}{\rm{\Delta }}{v}_{K}\end{eqnarray} \tag{ 16 }$$

(Weidenschilling 1977a; Nakagawa et al. 1986). In our simulation frame, their velocity also includes Keplerian shear:

$$\begin{eqnarray}&&{v}_{{sh}}(x)=-\displaystyle \frac{3}{2}{{\rm{\Omega }}}_{K}x\hat{y}.\end{eqnarray} \tag{ 17 }$$

In addition, particles are subject to the gravitational attraction of the core, as well as random kicks from the MRI turbulence.

Zoom In Zoom Out Reset image size

In simulations with AD, turbulence is relatively weak. Unless particles are strongly coupled, they largely follow trajectories in the laminar flow to zeroth order. This is best seen in the snapshots for ${\tau }_{s}=1$ particles shown in Figure 2(b), where bulk particle motion is dominated by radial drift and shear. Most particles are captured by the core as they pass by, leading to density enhancement at the center, while others exit the box. The overall process is relatively brief, and takes much less than $10{{\rm{\Omega }}}_{K}^{-1}$ because we run out of particles in the simulation box. For ${\tau }_{s}=0.1$ particles, while the overall process is similar, evolution (and hence accretion rate) is slower because the radial drift speed is lower for more strongly coupled particles. We also see that the distribution of particles is more blurry, again because these particles are more strongly coupled to the gas and hence are more strongly affected by turbulence.

Snapshots from ideal MHD simulations show dramatic differences. For both ${\tau }_{s}=0.1$ and ${\tau }_{s}=1$ particles shown in Figures 2(c) and (d), the motion of all particles is strongly perturbed by turbulence, with their velocity very effectively randomized. A large fraction of particles linger around the core for a much longer time. While accretion still occurs, it is much less obvious what fraction of particles is accreted.

Additionally, we see that in the AD simulation with ${\tau }_{s}=1$ (Figure 2(b)), the particle scale height is well below r_H, and hence pebble accretion proceeds essentially in a 2D manner: the entire particle layer is potentially within the reach of the core, and most of the physics can be understood by restricting particle orbits to the disk midplane. In other cases, however, with particle scale height ${H}_{p}\gtrsim {r}_{H}$, pebble accretion is inherently 3D. One might think that only particles located within $z=\pm {r}_{H}$ can potentially be accreted. However, because particle scale height is maintained by turbulent diffusion, individual particles follow complex trajectories and can travel through very different vertical heights during their passage of the core. As a result, it is not obvious what fraction of particles can be captured, and whether one can use theories developed in the 2D case (Ormel & Klahr 2010; Lambrechts & Johansen 2012) to predict the rate of pebble accretion. We will address this question through more detailed analysis.

To measure the instantaneous particle accretion rate, we first discuss the criterion for a particle to be counted as being accreted. We adopt an energy criterion based on the total energy of individual particles relative to the planetary core, defined as

$$\begin{eqnarray}{E}_{b} & = & {E}_{p}+{E}_{k}=\,-\,{{GM}}_{c}\displaystyle \frac{{r}^{2}+3{R}_{s}^{2}/2}{{({r}^{2}+{R}_{s}^{2})}^{3/2}}\\ & & +\displaystyle \frac{1}{2}({v}_{x}^{{\prime} 2}+{v}_{y}^{{\prime} 2}+{v}_{z}^{{\prime} 2}),\end{eqnarray} \tag{ 18 }$$

where ${v}_{x}^{{\prime} },{v}_{y}^{{\prime} },{v}_{z}^{{\prime} }$ are particle velocities with Keplerian shear motion subtracted. For accreted particles, the total energy is negative and gradually decreases as the particle spirals into the core due to the energy dissipation. For the escaped particles, the total energy will remain positive until the end of the simulation. There is a lower limit to the total energy, which is set by the softening length, given by setting r = 0 and E_k = 0:

$$\begin{eqnarray}&&{E}_{p0}=-{\mu }_{c}{\left.\displaystyle \frac{{r}^{2}+3{R}_{s}^{2}/2}{{({r}^{2}+{R}_{s}^{2})}^{3/2}}\right|}_{r=0}=-\displaystyle \frac{3{\mu }_{c}}{2{R}_{s}},\end{eqnarray} \tag{ 19 }$$

corresponding to particles falling into the center of the gravitational potential well. Because we do not resolve the Bondi radius of the planet, the depth of the potential well E_p0 is limited by the softening length R_s, associated with grid resolution.

We count a particle as being accreted when its total energy falls below an energy threshold, which is taken to be 3/4 of E_p0:

$$\begin{eqnarray}&&{E}_{b0}=\displaystyle \frac{3}{4}{E}_{p0}=-\displaystyle \frac{9{{GM}}_{c}}{8{R}_{s}}.\end{eqnarray} \tag{ 20 }$$

This corresponds to the gravitational potential energy at $r\approx 2{R}_{s}/3$, at which the softening of the gravitational potential becomes significant. With our choice of the softening length, we confirm that in almost all runs, the vast majority of particles of all sizes that meet this energy threshold end up being accreted into the core in all our simulations. The only exceptions involve ${\tau }_{s}=0.01$ particles in the runs idl3 and AD3, where some fraction ($\lesssim 20 \%$) of these particles that meet our criterion (20) are later kicked out of the core by the MRI turbulence. This is largely owing to the softening of the gravitational potential used in our simulations. In reality, we still expect many of these particles to fall into the core, although the accretion rates measured based on criterion (20) likely represent an overestimate.

We then determine the steady accretion stage to calculate the accretion rate. In Figure 3, we show the measured time sequence of instantaneous particle accretion rates for ${\tau }_{s}=0.1$ and ${\tau }_{s}=1$ particles (averaged over all cycles). Because particles are initially released without random velocities, it takes a short period of time ($\sim {\tau }_{s}^{-1}$) before they catch up with turbulent motion in the gas, as well as for them to respond to the gravity of the core. In addition, we expect a steady state to be achieved after particles initially located right outside the Hill sphere reach the core, which roughly takes a time of around ${r}_{H}/{\rm{\Delta }}{v}_{K}\sim {{\rm{\Omega }}}^{-1}$. Overall, it should take up to a few ${{\rm{\Omega }}}_{K}^{-1}$ (longer for more loosely coupled particles) for the particle accretion rate to reach a steady state. To better illustrate how a steady state is achieved, we have removed particles whose initial total energies are below the threshold energy (20). From Figure 3 we see that approximate steady-state accretion is indeed achieved for essentially all runs; it is exhibited as a plateau following a rapid rise in accretion rate. The steady-state period is relatively brief (a few ${{\rm{\Omega }}}_{K}^{-1}$); it is followed by a rapid drop in particle accretion rate as we run out of particles available within the simulation box.

Zoom In Zoom Out Reset image size

We use the accretion rate curves of hydrodynamic simulations as the standard of reference, and mark the beginning and ending times of these plateaus for each core mass and particle size as the time interval for steady accretion. The same interval is used for the corresponding AD and ideal MHD simulations, and we can see from Figures 3(a) and (b) that, in general, the time intervals chosen from pure hydrodynamic simulations match reasonably well for AD and ideal MHD simulations. We measure the steady-state particle accretion rate by averaging the accretion rates over the marked time intervals.

Although we include particles with stopping time as large as ${\tau }_{s}=10$ in our simulations, it takes much longer for these particles to respond to turbulence, during which time a large fraction of them have left our simulation box. Therefore, we consider our measured accretion rates for these particles to be less reliable. Moreover, since grains typically continue to grow up to about centimeter size (Birnstiel et al. 2011), where they have ${\tau }_{s}\lt 10$ in the bulk of an MMSN disk, we expect our study to cover most of the relevant parameter space for PPDs.

4.3.1. Absolute Accretion Rate Coefficient

A natural scale to normalize particle accretion rate in the Hill regime is (Lambrechts & Johansen 2012)

$$\begin{eqnarray}&&{\dot{M}}_{2D}=2{{\rm{\Sigma }}}_{p}{r}_{H}{v}_{H}=3{{\rm{\Sigma }}}_{p}{{\rm{\Omega }}}_{K}{r}_{H}^{2},\end{eqnarray} \tag{ 21 }$$

where ${{\rm{\Sigma }}}_{p}$ is the particle surface density, and ${v}_{H}\equiv (3/2){{\rm{\Omega }}}_{K}{r}_{H}$ is the Hill velocity, characterizing the velocity with which particles approach the core at the Hill radius. Effectively, it is the characteristic accretion rate in the 2D case, with accretion radius $\sim {r}_{H}$. We define the absolute accretion rate coefficient ${k}_{\mathrm{abs}}$ as

$$\begin{eqnarray}&&{k}_{\mathrm{abs}}=\displaystyle \frac{{\dot{M}}_{\mathrm{sim}}}{{\dot{M}}_{2D}}=\displaystyle \frac{{\dot{M}}_{\mathrm{sim}}}{3{{\rm{\Sigma }}}_{p}{{\rm{\Omega }}}_{K}{r}_{H}^{2}},\end{eqnarray} \tag{ 22 }$$

where ${\dot{M}}_{\mathrm{sim}}$ is the accretion rate calculated in our simulations. This coefficient is independent of particle size, and provides a direct, dimensionless measure of absolute particle accretion rates as a function of ${\tau }_{s}$.

4.3.2. 2D/3D-modified Accretion Rate Coefficient

When the thickness of the particle layer reaches or exceeds r_H, as mentioned earlier, the accretion process is inherently 3D. The effect is noted and distinguished in Morbidelli et al. (2015). Here, we provide another way to normalize the measured particle accretion rate that takes into account the finite thickness of the particle layer, which transitions smoothly from the 2D regime to the 3D one as H_p varies.

Our procedure is a direct generalization of the 2D normalization factor, by dividing the particle layer into 2D sheets at individual heights, applying (21), and then integrating over height:

$$\begin{eqnarray}{\dot{M}}_{\mathrm{mod}} & = & {\displaystyle \int }_{-{r}_{H}}^{{r}_{H}}2{\rho }_{p}(z){R}_{H}(z){v}_{\mathrm{sh}}(z){dz}\\ & = & {\displaystyle \int }_{-{r}_{H}}^{{r}_{H}}3{\rho }_{p}(z){{\rm{\Omega }}}_{K}({r}_{H}^{2}-{z}^{2}){dz},\end{eqnarray} \tag{ 23 }$$

where ${R}_{H}(z)\equiv \sqrt{{r}_{H}^{2}-{z}^{2}}$, ${v}_{{sh}}(z)\equiv (3/2){{\rm{\Omega }}}_{K}{R}_{H}(z)$, ${\rho }_{p}(z)\,={\rho }_{p0}{e}^{-{z}^{2}/(2{H}_{p}^{2})}$ is the local particle density, and ${\rho }_{p0}\,={{\rm{\Sigma }}}_{p}/(\sqrt{2\pi }{H}_{p})$ is the midplane particle density. Note that Morbidelli et al. (2015) used the mass flux ${\dot{M}}_{F}$ of pebbles due to radial drift as the normalization factor, which is dependent on particle size. Our choice of normalization factor is independent of ${\tau }_{s}$, which is more convenient for comparing the pebble accretion rate among different particle sizes.

Similar to Section 4.3.1, we now define the 2D/3D-modified accretion rate coefficient, ${k}_{\mathrm{mod}}$, as

$$\begin{eqnarray}&&{k}_{\mathrm{mod}}=\displaystyle \frac{{\dot{M}}_{\mathrm{sim}}}{{\dot{M}}_{\mathrm{mod}}}.\end{eqnarray} \tag{ 24 }$$

This coefficient better characterizes the intrinsic efficiency of pebble accretion when the particle layer is puffed up, and allows us to directly compare simulations with different levels of turbulence.

Theoretically, the particle accretion rate is largely determined by the accretion radius r_a, the maximum impact radius below which particles can be accreted by the core. Lambrechts & Johansen (2012) qualitatively derived scaling relations between r_a and particle stopping time ${\tau }_{s}$. In the so-called "drift regime" with core mass below transition mass M_t (13), the accretion radius is given by (see also Baruteau et al. 2016)

$$\begin{eqnarray}&&{r}_{a}\sim {\tau }_{s}^{1/2}{\left(\displaystyle \frac{{M}_{c}}{{M}_{t}}\right)}^{1/6}{r}_{H}\propto {\tau }_{s}^{1/2},\end{eqnarray} \tag{ 25 }$$

for ${\tau }_{s}\lesssim 1$. In the "Hill regime" that we consider here, these authors find

$$\begin{eqnarray}&&{r}_{a}\sim {\tau }_{s}^{1/3}{r}_{H},\end{eqnarray} \tag{ 26 }$$

again for ${\tau }_{s}\lesssim 1$. We provide scalings derived from these approximate asymptotic relations for comparison with our results.¹¹

The rate of pebble accretion also depends on the velocity ${v}_{\mathrm{rel}}$ at which particles approach the core. For ${\tau }_{s}\lesssim 1$, the particles are well coupled to the gas, so that ${v}_{\mathrm{rel}}$ is roughly the relative velocity between the core and the gas:

$$\begin{eqnarray}&&{v}_{\mathrm{rel}}=\max [{\rm{\Delta }}{v}_{K},{v}_{\mathrm{sh}}({r}_{a})].\end{eqnarray} \tag{ 27 }$$

In the drift and Hill regimes, ${v}_{\mathrm{rel}}$ is dominated by drift (${\rm{\Delta }}{v}_{K}$) or shear (${v}_{\mathrm{sh}}$), respectively. Since ${r}_{a}\lesssim {r}_{H}$ and ${v}_{\mathrm{sh}}\propto {r}_{a}$, even with ${M}_{c}\gt {M}_{t}$, drift can dominate for sufficiently strongly coupled particles.

The corresponding accretion rate in the 2D case is approximately given by

$$\begin{eqnarray}&&\dot{M}\approx 2{{\rm{\Sigma }}}_{p}{r}_{a}{v}_{\mathrm{rel}}.\end{eqnarray} \tag{ 28 }$$

Equations (25) and (26) with appropriate choices for v_rel imply that $\dot{M}\propto {\tau }_{s}^{1/2}$ in the drift regime, and $\dot{M}\propto {\tau }_{s}^{2/3}$ in the Hill regime. Note that the 2D accretion rate directly applies for a laminar disk. Generalization to 3D can be made using the 2D/3D-modified normalization factor described in Section 4.3.2. If the efficiency of pebble accretion is similar in 2D and in 3D, then we would expect the same scaling relation ($\dot{M}$ on ${\tau }_{s}$) to hold for ${k}_{\mathrm{mod}}$ when normalizing the measured accretion rates by (24).

However, there is no proof that such a generalization is valid, though it has been implicitly assumed to hold in the pebble accretion prescriptions in global models of planet formation (e.g., Morbidelli et al. 2015). The main caveat is that our normalization prescription is obtained essentially by assuming that particle motion is restricted at constant height, while in reality, particles follow complex trajectories as a result of turbulence. One of the main goals of this paper is to test whether pebble accretion efficiency remains the same in the turbulent 3D case as in the laminar 2D case.

In addition, given the qualitative nature of the derivation, the scaling relations are not necessarily very accurate, particularly near the transition between the drift and Hill regimes. More rigorous calculations can be found in Ormel & Klahr (2010), who calculated the accretion rates via direct integration of individual particle orbits. Therefore, a more meaningful comparison is to directly compare ${k}_{\mathrm{mod}}$ among the hydro, AD, and ideal MHD runs.

In Figure 4, we show the steady-state accretion rate coefficients measured from our simulations (averaged over all particle injection cycles). In our simulations with turbulence, individual cycles show significant temporal fluctuations in the pebble accretion rate. For the ideal MHD simulations in particular, no obvious steady state is achieved in individual cycles (this is also because of the limitation on the size of our simulation boxes, so that the time of continuous pebble flux is relatively short). With several cycles in each run, we expect that the standard deviation among the cycles approximately reflects the statistical uncertainties in our measurement of the pebble accretion rate. These are the uncertainties that we quote in Figure 4. We may underestimate the uncertainties for the pebble accretion rate because of the uncertainty within each of the cycles, e.g., the fluctuation of accretion rate within the steady-state accretion stage in each cycle.

Zoom In Zoom Out Reset image size

The absolute accretion coefficients shown in Figures 4(a) and (c) clearly indicate that there is a modest reduction in accretion rate in turbulence runs toward strongly coupled particles. Stronger turbulence leads to greater reduction. However, when looking at modified accretion rate coefficients shown in Figures 4(b) and (d), the reduction is mostly eliminated. This fact suggests that turbulence does not strongly affect the intrinsic efficiency of pebble accretion. The reduction in absolute accretion rate with stronger turbulence is mainly caused by larger particle scale height. Other effects cancel out, as we discuss in Section 5.2. This is the most important conclusion from this work.

Overall, the efficiency of pebble accretion peaks at $0.3\lesssim {\tau }_{s}\lesssim 3$, consistent with theoretical expectations. We also note from Figures 4(b) and (d) that the dependence of ${k}_{\mathrm{mod}}$ on ${\tau }_{s}$ at $\tau \lt 1$ in pure hydrodynamic runs follows more closely a power law of ${\tau }_{s}^{1/2}$ instead of ${\tau }_{s}^{2/3}$. The main reason for the apparent discrepancy is likely because small particles with $\tau \lesssim 0.3$ (0.03) in the ${M}_{c}=3{M}_{t}$ (30M_t) case are already in the drift regime of pebble accretion. Therefore, our simulations only probe a relatively narrow range in ${\tau }_{s}$ where the scaling relations are potentially applicable, and the difference between ${\tau }_{s}^{1/2}$ and ${\tau }_{s}^{2/3}$ is not very easily distinguishable. In addition, as discussed in Section 4.4, the analytical scaling relations are unlikely to be accurate near the transition from the drift regime to the Hill one.

With relatively small core mass ${M}_{c}=3\times {10}^{-3}{M}_{T}=3{M}_{t}$ (near the lower mass end of the Hill regime), we see that ${k}_{\mathrm{mod}}$ is slightly reduced toward smaller ${\tau }_{s}$ in runs with turbulence compared with the pure hydrodynamic run. Therefore, turbulence still leads to a small-to-modest reduction in the pebble accretion efficiency toward more strongly coupled particles (by a factor of at most 2–3 for ${\tau }_{s}\gtrsim 0.01$). With larger core mass (i.e., ${M}_{c}\sim 30{M}_{t}$), on the other hand, we see that ${k}_{\mathrm{mod}}$ is insensitive to the strength of the MRI turbulence at least down to ${\tau }_{s}=0.01$. In other words, stronger gravity from more massive cores helps overwhelm the effect of turbulence, as one would naturally expect.

It is well known that pebble accretion is the most efficient for marginally coupled particles with $0.1\lesssim \tau \lesssim 1$ in the Hill regime. We find that in some cases MRI turbulence even slightly enhances the accretion rate of particles in this size range. The reason will be discussed in the following subsections with more detailed analysis.

The fact that the efficiency of pebble accretion is largely unaffected by turbulence is somewhat surprising, especially given that particles follow complex trajectories in the 3D case. In this subsection, we analyze the trajectories of accreted and non-accreted particles and address how turbulence affect the process of pebble accretion.

5.2.1. Particle Trajectories

In Figure 5, we show some typical trajectories of accreted particles with ${\tau }_{s}=0.1$ and 1. Particles in each panel are randomly chosen to present the trajectories, so that the numbers of different "types" of lines shown in the figure roughly reflect the actual probability of each "type" of trajectory being populated. The Hill sphere is marked by a black circle in each panel. Projected to the x–y plane (lower panels in each subfigure), particle trajectories are color-coded based on their initial positions $({x}_{0},{y}_{0})$ marked by the black dots, with warm (cold) colors corresponding to ${y}_{0}\lt 0$ (${y}_{0}\gt 0$). We note that while initial particle distributions are spatially uniform, in practice, the pebble accretion process itself would inhibit some particle trajectories coming from the ${y}_{0}\lt 0$ region (which is one limitation of our local simulations). Nevertheless, we note from Figure 5 that these particles represent a minority among accreted particles, and we confirm that their contribution to our measured pebble accretion rates is well below $\sim 50 \%$. Projected to the x–z plane (upper panels in each subfigure), each trajectory is colored from cold (blue) to warm (red) according to time.

Zoom In Zoom Out Reset image size

We first see that for pure hydrodynamic runs, because particle orbits are sub-Keplerian and undergo radial drift, most accreted particles approach the core from the upper right side in each panel (${x}_{0}\gt 0,{y}_{0}\gt 0$) corresponding to the upwind direction, and follow a parabolic shape before reaching the core due to shear. As core mass increases, on the other hand, particles from the lower side can be captured as well, which is a result of a larger Hill radius (${\rm{\Omega }}{r}_{H}\gt {\rm{\Delta }}{v}_{K}$), and can be deduced by directly integrating individual particle orbits (Ormel & Klahr 2010).

Turbulence strongly affects particle trajectories in both the horizontal and the vertical dimensions. In AD simulations with weak turbulence, while most trajectories are similar to those in the hydrodynamic run, we already see that the trajectories are more spatially spread, and more particles initially located in ${y}_{0}\lt 0$ regions end up being accreted, implying that the zone feeding the planetary core becomes broader in the presence of turbulence. In ideal MHD simulations, particles show much more significant random motion, and their trajectories occupy a much larger volume of our simulation box (which requires us to choose a relatively large and broad box). Many particles exhibit swirling trajectories around $x\sim 0$ before reaching the core as they are scattered by turbulence. Overall, we see that stronger turbulence can bring particles from a broader range of locations to the core, which can potentially enhance the rate of pebble accretion.

On the other hand, not all particles that enter the Hill sphere can be accreted, especially in the presence of strong turbulence. In Figure 6, we show some typical trajectories of particles that have entered the Hill sphere but eventually escaped capture by the core. We again choose particle stopping times ${\tau }_{s}=0.1$ and 1, whose accretion radii r_a are expected to be close to r_H. We show the trajectories only from runs AD2 and idl3 because they represent the two extreme situations with the most significant contrast. In run AD2, with the combination of a high-mass core and weak turbulence, we see that only a very small fraction of particles that enter the core escape, and these particles all barely touch the Hill sphere at the edge before exiting. By contrast, in run idl3, with the combination of a low-mass core and strong turbulence, we see that a large population of particles enter the Hill sphere from all directions. Their trajectories fill the entire Hill sphere yet they simply pass the core without being captured. Overall, we see that stronger turbulence reduces the probability for particles that enter the Hill sphere to be accreted.

Zoom In Zoom Out Reset image size

In addition, we also see from Figures 5 and 6 that when ${H}_{p}\gtrsim {r}_{H}$, particle trajectories experience significant vertical motion. This is the case for both accreted and non-accreted particles that enter the Hill sphere. The trajectories of ${\tau }_{s}=0.1$ and 1 particles shown here span the entire vertical thickness of the particle layer. Therefore, even when accretion proceeds in a 3D manner, the core has potential access to a large fraction of particles in the entire particle column, instead of just those initially located within $| z| \lesssim {r}_{H}$ (as assumed in evaluating the normalization factor in Section 4.3.2).

In sum, we have shown that turbulence affects particle trajectories in ways that can both enhance and reduce the rate of pebble accretion. The outcome of the pebble accretion efficiency is determined by the level to which these effects cancel. Enhancement or reduction of pebble accretion efficiency observed in Figure 4(b) can be understood as imperfect cancellation, but overall, these effects cancel to a large degree and this leads to a more or less unchanged efficiency for pebble accretion in the presence of turbulence compared with the laminar case.

5.2.2. Map of Accretion Probability

We can reinforce the conclusion reached in the previous subsection by looking at the map of accretion probability, shown in Figure 7 for particle sizes ${\tau }_{s}=0.1$ and 1. The probability map is obtained by identifying the initial positions of particles that are accreted to the core, binning them to the grid, and normalizing them by the initial particle density. In practice, they are projected to the x–y plane and the x–z plane.

Zoom In Zoom Out Reset image size

For each core mass and particle size, the accretion probability in hydrodynamic simulations is either 0 or 1, with clear boundaries that indicate the shape of the feeding zones. Higher core mass leads to a wider feeding zone due to a larger Hill radius, and also enables the accretion of particles initially located on the downwind side of the flow, as mentioned earlier. The ${\tau }_{s}=1$ particles also have a slightly wider feeding zone than ${\tau }_{s}=0.1$ particles, which is consistent with theoretical expectations (see Section 4.4). Besides, the feeding zones of ${\tau }_{s}=0.1$ particles are more aligned with the y-axis due to the particles' smaller radial drift velocities.

As the turbulence gets stronger, the map of accretion probability becomes smoother, and the feeding zone becomes broader. This is accompanied by the reduction in accretion probability within the feeding zone, especially for the ideal MHD case. Because our measurements show that k_mod in the ideal MHD simulations is comparable with that in the hydrodynamic case, the two aforementioned effects approximately cancel each other. More specifically, with our fiducial core mass ($\mu =3\times {10}^{-3}$), we find that for ${\tau }_{s}=0.1$ particles, the fractions of particles entering the Hill sphere that are eventually accreted are ${f}_{\mathrm{hyd}3}=0.77$, ${f}_{\mathrm{AD}3}=0.61$, and ${f}_{\mathrm{idl}3}=0.11$ in our hydrodynamic, AD and ideal MHD simulations, respectively; for ${\tau }_{s}=1$ particles, the corresponding fractions are ${f}_{\mathrm{hyd}3}=0.91$, ${f}_{\mathrm{AD}3}=0.94$, and ${f}_{\mathrm{idl}3}=0.66$. In the vertical direction, we see again that the core accretes particles initially located over most of the vertical column.

Overall, our analysis of accretion probability supports our conclusion that the efficiency of pebble accretion is not strongly affected by turbulence because positive and negative effects largely cancel each other.