Gas-assisted Growth of Protoplanets in a Turbulent Medium

We begin by discussing generally how we model growth of planets in the presence of nebular gas. The details of how specific quantities are calculated are deferred to subsequent sections.

Our calculation proceeds in an order-of-magnitude manner—i.e., the approximations made and the neglected effects mean that quoted values should be correct to within an order of magnitude. In what follows, we consider two types of bodies. The large bodies, or protoplanetary cores, are assumed to be massive enough that they are unaffected by gas drag and thus move at the local Keplerian orbital velocity. This constrains our cores to have radii ≳10 km, in which case their velocity will deviate from Keplerian by 10⁻³ of the gas velocity, at most, with a weak dependence on stellocentric distance for fiducial disk parameters. In practice, we rarely consider cores of such small size. We note also that, while these cores are insensitive to aerodynamical gas drag, large bodies with masses in the range ${10}^{21}\,{\rm{g}}\lt M\lt {10}^{25}\,{\rm{g}}$ can still be affected by gas dynamical friction (Grishin & Perets 2016). We do not include these effects here. The second type of particles considered are "small" bodies, which can be substantially affected by gas drag. The growth timescale in pebble accretion is strongly dependent on the size of the small body under consideration, unlike canonical core accretion where the size of the planetesimals enters only through its effect on the small body velocity dispersion. Thus, all calculations are performed as a function of small body radius, r_s. Quantities of interest can later be averaged over size by assuming a size distribution for the small bodies. Note that quoted values of the growth timescale t_grow implicitly assume that all of the surface density is contained in particles of the given value of r_s. For a size distribution where most of the mass is in the largest sizes of particles present (e.g., a Dohnyani size distribution, Dohnanyi 1969), this value of t_grow is approximately equal to the growth timescale for a distribution of small body sizes where the maximum size present is the given value of r_s. In this paper, we will not explicitly perform integrations over small body size; see our Paper II for examples of this process, as well as a discussion of the effects of altering the size distribution.

In gas-assisted growth, the interaction between the small body and the nebular gas substantially modifies the accretion process. As the small body approaches the core, gas drag will dissipate the kinetic energy of the small body relative to the large body. This loss of energy can cause small bodies on non-collisional trajectories to become bound to the core and eventually be accreted. This process is similar to the ${L}^{2}s$ mechanism identified by Goldreich et al. (2002) for formation of Kuiper Belt binaries, with gas drag as the source of dissipation in the place of dynamical friction. Gas drag can also stop small bodies from accreting—if particles couple strongly to the gas as they flow around the core, then they will be unable to accrete.

Which of these processes occur(s) depends on the relative size of two different length scales: the stability radius, R_stab, and the Bondi radius, R_b. The stability radius is the smallest radius at which stable orbits by the small body about the large body are possible: outside of R_stab, interactions between the small body and either the nebular gas or the central star will shear the small body away from the large body's gravity. Within R_stab, the small body can safely inspiral onto the core. The details of how R_stab is calculated are discussed in Section 4.1. The Bondi radius, on the other hand, is approximately the radius at which the escape velocity from the core is equal to the speed of sound in the gas:

$$\begin{eqnarray}&&{R}_{b}=\displaystyle \frac{{GM}}{{c}_{s}^{2}},\end{eqnarray} \tag{ 1 }$$

where M is the mass of the core and ${c}_{s}=\sqrt{{kT}/\mu }$ is the isothermal sound speed of the gas. Here, k is Boltzmann's constant and μ is the mean molecular weight of the nebular gas. We consider the Bondi radius because it tells us roughly the length scale interior to which the the core can have a stable atmosphere, which has substantial effects on the flow pattern. For the lowest-mass cores we consider, the Bondi radius may be less than the physical radius of the core, R. This occurs roughly at a core mass of:

$$\begin{eqnarray}&&{M}_{a}\equiv \displaystyle \frac{{c}_{s}^{3}}{G}{\left(\displaystyle \frac{3}{4\pi G{\rho }_{p}}\right)}^{1/2}\\ &&\,\approx \,2\times {10}^{-4}{M}_{\oplus }{\left(\displaystyle \frac{a}{30\mathrm{au}}\right)}^{-9/14}{\left(\displaystyle \frac{{\rho }_{p}}{2{\rm{g}}{\mathrm{cm}}^{-3}}\right)}^{-1/2},\end{eqnarray} \tag{ 2 }$$

where ρ_p is the density of the protoplanet (e.g., Rafikov 2006), and we have used our fiducial disk parameters (see Section 5.1) for the expression in the second line. If ${R}_{b}\lt R$, then the effects discussed below are unchanged, with R taking the place of R_b. In what follows, we will discuss accretion for R_b < R_H, where R_H is the core's Hill radius (see Section 4.1.1). We discuss accretion in the regime R_b > R_H in Section 7.

Given these considerations, we center our model around two main ideas about accretion in the presence of gas, which are summarized in Figure 1:

• 1.  
  If the radius for stable orbits exceeds the Bondi radius, i.e., ${R}_{\mathrm{stab}}\gt {R}_{b}$, then the flow pattern of gas is not substantially altered in the region where particles can stably orbit the core. In this case, any small bodies that deplete their kinetic energy relative to the core within R_stab will inspiral onto the core and be accreted. On the other hand, any particles that are unable to dissipate their kinetic energy in this regime will pass out of R_stab and will not be accreted.
• 2.  
  If instead the Bondi radius exceeds the stable orbit radius, i.e., R_b > R_stab, then particles that are able to deplete their kinetic energy relative to the core will not accrete. This is due to the fact that, in this regime, well-coupled particles will tend to flow around the core's atmosphere, which extends up to R_b. If particles are instead unable to dissipate their kinetic energy, such that they penetrate into the atmospheric radius, the increase in density as the particle enters the growing planet's atmosphere is taken to be so substantial that the particle will now be able to dissipate its kinetic energy and will accrete onto the core.

Zoom In Zoom Out Reset image size

The first point is supported not only by order of magnitude considerations, but also by detailed numerical simulations of the growth of protoplanets in the presence of gas (e.g., Ormel & Klahr 2010; Lambrechts & Johansen 2012). Analytic calculations also show that even small bodies several orders of magnitude larger than the sizes with which we will be concerned inspiral on times shorter than the disk lifetime. Thus, we are justified in neglecting this part of the accretion process (Perets & Murray-Clay 2011). We also neglect the possibility that the core's envelope is periodically replenished by the protoplanetary disk; see Ormel et al. (2015) for a discussion of this possibility, and Alibert (2017) for an application of this replenishment to pebble accretion.

Our second point invokes the classical solution of flow around an obstacle: for R_b < R_H, the flow of the nebular gas is subsonic, meaning the core's atmosphere is approximately incompressible. See Ormel (2013) Figure 5(A) for an example of this flow pattern in the context of a planet embedded in a protoplanetary disk. We note here that more detailed simulations of the flow show structure that may produce circulation into R_b (Ormel 2013; Fung et al. 2015), which we assume we can neglect for the level of accuracy we desire in this model. In addition, we do not explicitly calculate the work done on particles interior to R_b—our assumption that particles with ${KE}\gt W$ will always be accreted for ${R}_{\mathrm{stab}}\lt {R}_{b}$ will eventually be violated for large enough particle sizes. See Inaba & Ikoma (2003) for an in-depth discussion of accretion in this regime.

While the criteria above are used throughout parameter space, in order to better illustrate how we model gas-assisted growth, we also provide a simplified "sketch" of how accretion in our model operates—one that accurately describes pebble accretion over a large amount of parameter space.⁵ Figure 2 shows a "cartoon" of gas-assisted growth for a fixed core mass and semimajor axis in the disk, with small body radius increasing from left to right. Each panel shows the core's Bondi radius, R_b, as well as the two radii that are used to determine R_stab—the Hill radius R_H (see Section 4.1.1), which is the distance past which the stellar gravity will pull particles off the core, and the radius ${R}_{\mathrm{WS}}^{{\prime} }$ (see Section 4.1.2), beyond which gas drag pulls particles off the core. Because smaller particles have a higher ratio of surface area to volume and therefore experience larger gas drag accelerations, ${R}_{\mathrm{WS}}^{{\prime} }$ is smaller for smaller particles. In the far left panel, the particles have low mass, meaning gas drag can easily pull them off of the core, and ${R}_{\mathrm{WS}}^{{\prime} }$ lies inside the core's atmosphere. Because these particles have low mass, they easily dissipate their kinetic energy during the encounter with the core, meaning they are in the regime in the lower left-hand panel of Figure 1 and do not accrete. As particle size grows, ${R}_{\mathrm{WS}}^{{\prime} }$ increases as well, until it exceeds R_b, the scale of the core's atmosphere; i.e., there now exists a region exterior to the core's atmosphere where particles can stably orbit the core. Because these particles are still of relatively low mass, they are able to dissipate their kinetic energy through gas drag. We therefore fall into the upper left-hand panel of Figure 1, which signals the onset of pebble accretion (middle panel of Figure 2). Finally, as particle size continues to increase, we eventually reach a point where the particles are so massive that they no longer dissipate their kinetic energy. As particle size increases, ${R}_{\mathrm{WS}}^{{\prime} }$ will continue to increase, meaning we clearly still have R_stab > R_b. These particles are therefore in the upper right-hand panel of Figure 14 and will not accrete (right-hand panel of Figure 2).

Zoom In Zoom Out Reset image size

Given this formalism, for the purposes of discussing whether a given small body will accrete, we simply need to compare the magnitude of the kinetic energy of the particle relative to the core and the work done on the particle during its encounter with the core. The kinetic energy of the particle relative to the core before the encounter is

$$\begin{eqnarray}&&{KE}=\displaystyle \frac{1}{2}{{mv}}_{\infty }^{2},\end{eqnarray} \tag{ 3 }$$

while we take the work done by gas drag to be simply

$$\begin{eqnarray}&&W=2{R}_{\mathrm{acc}}{F}_{D}({v}_{\mathrm{enc}}),\end{eqnarray} \tag{ 4 }$$

where ${R}_{\mathrm{acc}}=\max ({R}_{\mathrm{stab}},{R}_{b})$, and ${F}_{D}({v}_{\mathrm{enc}})$ is the drag force on a small body moving at a velocity v_enc, which is the relative velocity between the small body and the large body during the encounter. Section 3.6 gives our calculation of v_enc and a discussion of why R_acc is used for determining the work done by gas drag.

Our consideration of the ranges of particle sizes that can be accreted is an important aspect to our modeling that is often omitted in other works. See Section 5.5 for a detailed discussion.

Given the uncertainties in the size distribution of small bodies present in protoplanetary disks, understanding the extent of particle radii for which gas-assisted growth is possible is an important facet of studying the role of pebble accretion in planet formation. Furthermore, an important effect of nebular turbulence is to change the range of small body sizes that can be accreted (see Section 6.1), which makes a detailed consideration of the small body sizes where pebble accretion can operate an important aspect of our model.

In summary, a particle will be able to accrete if one of the following criteria are met:

• 1.  
  ${R}_{\mathrm{stab}}\gt {R}_{b}$ and $2{F}_{D}({v}_{\mathrm{enc}}){R}_{\mathrm{acc}}\gt \tfrac{1}{2}{{mv}}_{\infty }^{2}$
• 2.  
  ${R}_{\mathrm{stab}}\lt {R}_{b}$ and $2{F}_{D}({v}_{\mathrm{enc}}){R}_{\mathrm{acc}}\lt \tfrac{1}{2}{{mv}}_{\infty }^{2}.$

Small bodies that do not satisfy either of these criteria will not be able to accrete, i.e., we set ${t}_{\mathrm{grow}}=\infty$ for these particles. We emphasize that setting ${t}_{\mathrm{grow}}=\infty$ refers only to the timescale for growth by pebble accretion: it is possible these particles could still be accreted via less-efficient processes, such as capture by gravitational focusing unassisted by gas drag or by collisions unaffected by gravity (see Equation (128)). In other words, having ${t}_{\mathrm{grow}}=\infty$ does not necessarily imply that growth actually halts. See Section 5.5 for further discussion. While our model could, in principle, be extended to include these effects, in this work we are concerned primarily with gas-assisted growth; therefore, we do not explicitly include these other timescales.

We will now elaborate upon how the growth timescales for cores are computed. In order to calculate the growth timescale for the large bodies for a given core mass M, we use the usual expression (see, e.g., Goldreich et al. 2004, hereafter GLS)

$$\begin{eqnarray}&&{t}_{\mathrm{grow}}\equiv {\left(\displaystyle \frac{1}{M}\displaystyle \frac{{dM}}{{dt}}\right)}^{-1}.\end{eqnarray} \tag{ 5 }$$

The rate at which small bodies encounter the core is given by $n{\sigma }_{\mathrm{acc}}{v}_{\infty }$, where n is the volumetric number density of small bodies of a given size, σ_acc is the cross section for accretion of the small body by the large body, and v_∞ is the velocity at which small bodies encounter the large body. Note that ${v}_{\infty }$ is not necessarily the velocity of the small body during its encounter with the core, because this encounter can change the relative velocity; ${v}_{\infty }$ is the relative velocity between the two bodies at large separations. The number density of solids is simply $n={\rho }_{s}/m={f}_{s}{\rm{\Sigma }}/(2\,{{mH}}_{p})$ where ρ_s is the mass density of the small bodies, Σ is the surface density of the gaseous component of the disk, f_s is the solid-to-gas mass ratio, m is the mass of the small body, and H_p is the scale height of solids in the disk. Because each accretion of a small body increases the mass of the large body by m, the growth timescale is given by

$$\begin{eqnarray}&&{t}_{\mathrm{grow}}=\displaystyle \frac{{{M\; H}}_{p}}{2{f}_{s}{\rm{\Sigma }}{v}_{\infty }{R}_{\mathrm{acc}}{H}_{\mathrm{acc}}},\end{eqnarray} \tag{ 6 }$$

where we have decomposed ${\sigma }_{\mathrm{acc}}$ into the product of length scales in the plane of the disk and perpendicular to it: ${\sigma }_{\mathrm{acc}}=(2{R}_{\mathrm{acc}})(2{H}_{\mathrm{acc}})$. Thus, the aim of our calculation is to determine the quantities R_acc, H_acc, H_p, and v_∞. Once these quantities are known, we can immediately determine the growth timescale. The role each of these quantities plays in the growth timescale is illustrated graphically in Figure 3.

Zoom In Zoom Out Reset image size

In the next few sections, we will expound upon how each of the above quantities are calculated in the context of pebble accretion. Our algorithm is summarized in Appendix A.

Gas drag is generally treated by breaking the drag force, F_D, into a number of different regimes. We summarize our choices below; for a more in-depth discussion, see Batchelor (2000). First, we distinguish between the "diffuse regime," which applies for particles with ${r}_{s}\lt 9\lambda /4$, and the "fluid regime," r_s > 9λ/4. Here, r_s is the radius of the small body and λ is the mean free path of the gas particles. In the diffuse, non-supersonic,⁶ regime, the drag force on the particle is given by the Epstein drag law

$$\begin{eqnarray}&&{F}_{D}=\displaystyle \frac{4}{3}\pi {\rho }_{g}{v}_{\mathrm{th}}{v}_{\mathrm{rel}}{r}_{s}^{2},\end{eqnarray} \tag{ 7 }$$

where ${\rho }_{g}={\rm{\Sigma }}/(2{H}_{g})$ is the density of the gas, ${v}_{\mathrm{th}}=\sqrt{8/\pi }\,{c}_{s}$ is the thermal velocity of the gas, and v_rel is the relative velocity between the gas and the object. For a particle in the fluid regime, we must consider an additional parameter—the Reynolds number of the particle, given by ${Re}\,=2{r}_{s}{v}_{\mathrm{rel}}/(0.5\,{v}_{\mathrm{th}}\lambda )$. For ${Re}\,\lt \,1$, the particle is in the Stokes drag regime, and the drag force is given by

$$\begin{eqnarray}&&{F}_{D}=3\pi {\rho }_{g}{v}_{\mathrm{th}}{v}_{\mathrm{rel}}\lambda {r}_{s},\end{eqnarray} \tag{ 8 }$$

where λ is the mean free path of the gas particles. For Re > 1, the particle is in the Ram pressure regime, and

$$\begin{eqnarray}&&{F}_{D}=\displaystyle \frac{1}{2}{\rho }_{g}\pi {r}_{s}^{2}{v}_{\mathrm{rel}}^{2}.\end{eqnarray} \tag{ 9 }$$

To mitigate discontinuities in the drag force, we use a smoothed drag force law in the fluid regime given by Cheng (2009):

$$\begin{eqnarray}&&{F}_{D}=\displaystyle \frac{1}{2}{C}_{D}({Re})\pi {r}_{s}^{2}{\rho }_{g}{v}_{\mathrm{rel}}^{2},\end{eqnarray} \tag{ 10 }$$

where

$$\begin{eqnarray}{C}_{D}({Re}) & = & \displaystyle \frac{24}{{Re}}{(1+0.27{Re})}^{0.43}\\ & & +0.47[1-\exp (-0.04{{Re}}^{0.38})].\end{eqnarray} \tag{ 11 }$$

As ${Re}\,\to \,0$, we have ${C}_{D}\to 24/{Re}$, so F_D reduces to the Stokes drag law given in (8). As ${Re}\,\to \,\infty$, ${C}_{D}\to 0.47$, in which case F_D becomes the Ram drag force given in (9) (with a slightly different prefactor).

We can define a timescale from the drag force, known as the "stopping time" of the particle

$$\begin{eqnarray}&&{t}_{s}\equiv \displaystyle \frac{{{mv}}_{\mathrm{rel}}}{{F}_{D}}.\end{eqnarray} \tag{ 12 }$$

The Epstein and Stokes drag laws are both linear in velocity—in these linear drag regimes, it is straightforward to calculate, from Equations (7) and (8), for spherical particles of uniform density ρ_s,

$$\begin{eqnarray}{t}_{s}=\left\{\begin{array}{ll}\displaystyle \frac{{\rho }_{s}}{{\rho }_{g}}\displaystyle \frac{{r}_{s}}{{v}_{\mathrm{th}}}, & \mathrm{Epstein}\\ \displaystyle \frac{4}{9}\displaystyle \frac{{\rho }_{s}}{{\rho }_{g}}\displaystyle \frac{{r}_{s}^{2}}{{v}_{\mathrm{th}}\lambda }, & \mathrm{Stokes}\end{array}\right.\end{eqnarray} \tag{ 13 }$$

In these regimes, the stopping time of the particle is a function only of the properties of the particle and the gas; in particular, it is independent of the particle's velocity. Hence, t_s is often used as a parameterization of the particle's size, r_s, in terms of how well the small body is coupled to the gas flow (e.g., Chiang & Youdin 2010). If the drag law is instead quadratic in velocity, as in Equation (10), we numerically solve for the stopping time using Equations (18), (23), and (25). In practice, we solve these equations iteratively to calculate a self-consistent solution.

Due to the internal pressure of the gas, the gaseous component of the protoplanetary disk will move at a sub-Keplerian orbital velocity. Weidenschilling (1977) gives the difference in velocity Δv as Δv ≈ ηv_k, where v_k is local Keplerian orbital velocity, v_k = aΩ, and η is a measure of the local gas pressure support, with approximate value $\eta \,\approx {c}_{s}^{2}/(2{v}_{k}^{2})$. Due to this sub-Keplerian rotation, small bodies experience a "headwind" from the gas, which produces a drag force on the small bodies. If we use a polar coordinate system such that $\hat{{\boldsymbol{r}}}$ denotes the direction pointing away from the central star and $\hat{{\boldsymbol{\phi }}}$ denotes the direction of the disk's rotation, then the drag force causes the small bodies to move with a sub-Keplerian velocity in the $\hat{{\boldsymbol{\phi }}}$ direction, and to drift in the $-\hat{{\boldsymbol{r}}}$ direction. In the above notation, the particle acquires a laminar velocity relative to the gas, given by

$$\begin{eqnarray}&&{v}_{r,\mathrm{gas}}=-2\eta {v}_{k}\left[\displaystyle \frac{{\tau }_{s}}{1+{\tau }_{s}^{2}}\right],\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&{v}_{\phi ,\mathrm{gas}}=-\eta {v}_{k}\left[\displaystyle \frac{1}{1+{\tau }_{s}^{2}}-1\right],\end{eqnarray} \tag{ 15 }$$

where ${\tau }_{s}={t}_{s}{\rm{\Omega }}$. See Nakagawa et al. (1986) for further details.⁷

Because the laminar gas velocity relative to Keplerian is simply ${v}_{\mathrm{gas},k}=-\eta {v}_{k}\hat{{\boldsymbol{\phi }}}$, the velocity of the particle relative to Keplerian is

$$\begin{eqnarray}&&{v}_{r,k}=-2\eta {v}_{k}\left[\displaystyle \frac{{\tau }_{s}}{1+{\tau }_{s}^{2}}\right],\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}&&{v}_{\phi ,k}=-\eta {v}_{k}\left[\displaystyle \frac{1}{1+{\tau }_{s}^{2}}\right].\end{eqnarray} \tag{ 17 }$$

It is straightforward to show that the magnitude of the laminar component of the particle's velocity is

$$\begin{eqnarray}&&{v}_{{pg},{\ell }}=\eta {v}_{k}{\tau }_{s}\displaystyle \frac{\sqrt{4+{\tau }_{s}^{2}}}{1+{\tau }_{s}^{2}},\end{eqnarray} \tag{ 18 }$$

relative to the gas, and

$$\begin{eqnarray}&&{v}_{{pk},{\ell }}=\eta {v}_{k}\displaystyle \frac{\sqrt{1+4{\tau }_{s}^{2}}}{1+{\tau }_{s}^{2}},\end{eqnarray} \tag{ 19 }$$

relative to Keplerian.

In order to describe the strength of the turbulence in the disk, we use the standard Shakura-Sunyaev α parameterization of the effective kinematic viscosity. We employ α simply as a convenient parameterization of the strength of turbulence; for our purposes, α is fundamentally a local quantity and not necessarily connected with the accretion rate onto the star. In terms of α, the effective kinematic viscosity of the turbulent gas, ν_t, is given by (Shakura & Sunyaev 1973):

$$\begin{eqnarray}&&{\nu }_{t}=\alpha {c}_{s}{H}_{g},\end{eqnarray} \tag{ 20 }$$

where ${H}_{g}={c}_{s}/{\rm{\Omega }}$ is the scale height of the gas. If we write the viscosity as the product of the local turbulent velocity and the largest-scale turbulent eddies, ${\nu }_{t}={v}_{t}{{\ell }}_{t}$, and have ${l}_{t}\approx {v}_{t}/{\rm{\Omega }}$, then ${v}_{t}=\sqrt{\alpha }{c}_{s}$.

We use the Kolmogorov theory of turbulence to determine the turbulent energy spectrum. As described in Cuzzi & Hogan (2003), per the Kolmogorov theory, turbulence exists over a range of scales or "eddies," which are characterized by their wave number ${k}_{{\ell }}=1/{\ell }$. The largest-scale eddies, which occur on the length scale l_t, are the scale on which energy is supplied by the turbulence; these large-scale eddies "turn over" and transfer their energy to smaller-scale eddies until energy is finally dissipated by kinematic viscosity on some smallest scale $\tilde{\eta }$. If we assume that the rate of energy transfer between scales is independent of eddy size, we can show that the energy spectrum of the turbulence is given by $E(k)\propto {k}^{-5/3}$, where E(k) has units of energy per mass per wavenumber. The velocity associated with k is then $v{(k)=(2{kE}(k))}^{1/2}\propto {k}^{-1/3}$, and the overturn time for an eddy of wavenumber k is ${t}_{k}=1/({kv}(k))\propto {k}^{-2/3}$. Thus, the larger-scale eddies overturn more slowly and contain more of the turbulent kinetic energy. The size of the smallest-scale eddies can be determined by setting the rate of energy loss from molecular viscosity equal to the eddy turnover time, which gives $\tilde{\eta }={({\nu }^{3}/\epsilon )}^{1/4}$, where $\nu \sim {v}_{\mathrm{th}}\lambda$ is the molecular viscosity and the rate of energy dissipation. The length $\tilde{\eta }$ is known as the "Kolmogorov microscale" of length. We can determine the relation between the smallest- and largest-scale eddies by setting the rate that energy is supplied at the largest scales, ${v}_{t}^{2}/{t}_{L}\sim {v}_{t}^{3}/{l}_{t}$ equal to , the rate of energy dissipation at small scales. Plugging in our expression for in terms of $\tilde{\eta }$ and ν gives $\tilde{\eta }/{l}_{t}\,\sim {({v}_{t}{l}_{t}/\nu )}^{-3/4}={{Re}}_{t}^{-3/4}$, where we have defined the Reynolds number of the turbulence, ${{Re}}_{t}\,\equiv \,{\nu }_{t}/\nu$. In terms of α, we have ${{Re}}_{t}=\alpha {c}_{s}{H}_{g}/({v}_{\mathrm{th}}\lambda )$.

Qualitatively, the behavior of a small body in response to an eddy of wavenumber k depends on the ratio of the particle's stopping time to the eddy turnover time t_k. Particles with ${t}_{s}\lt {t}_{k}$ will come to equilibrium with the eddy before it turns over, and will follow the large-scale motion of the eddy. On the other hand, particles with ${t}_{s}\gt {t}_{k}$ will not come to equilibrium with the eddy before it turns over, and will therefore only receive a small perturbative kick from the eddy. To characterize the small body's response to the turbulence, most authors use the Stokes number, ${St}\equiv {t}_{s}/{t}_{L}$, where t_L is the overturn time of the largest eddies. We make the usual assumption that ${t}_{L}={{\rm{\Omega }}}^{-1}$, in which case ${St}={\tau }_{s}$, the parameter used in describing the particle's interaction with the laminar gas flow. This permits the use of one parameter, which we will hereafter refer to as St, to characterize the particle's interaction with both the laminar and turbulent components of the nebular gas flow. See Youdin & Lithwick (2007) for a more in-depth discussion of the effect of varying t_L.

Due to their interaction with the turbulent gas, small bodies will move relative to inertial space with some non-zero root-mean-square (rms) velocity, ${v}_{{pi},t}$. The rms velocity can be calculated through order-of-magnitude means, as follows (Youdin & Lithwick 2007): For ${St}\gg 1$, we have ${t}_{s}\gg {t}_{L};$ in this regime, particles receive many uncorrelated "kicks" from the largest-scale eddies over a single stopping time, causing the particle to experience a random walk in velocity. Because the particle's velocity damps out over a time t_s, the particle receives approximately $N\sim {t}_{s}/{t}_{L}\sim {St}$ kicks of amplitude $v\sim {v}_{t}\,{t}_{L}/{t}_{s}={v}_{t}/{St}$, resulting in a velocity ${v}_{{pi},t}\sim {v}_{t}/\sqrt{{St}}$. On the other hand, for ${St}\ll 1$, we expect ${v}_{{pi},t}\sim {v}_{t}$. The simplest expression that has the correct behavior in each of these limits is given by

$$\begin{eqnarray}&&{v}_{{pi},t}=\displaystyle \frac{{v}_{t}}{\sqrt{1+{St}}}.\end{eqnarray} \tag{ 21 }$$

In order to better calculate the velocities of smaller particles, i.e., particles with ${St}\lt 1$, we employ expressions from Ormel & Cuzzi (2007, hereafter OC07), who give closed-form equations for the rms velocity of solid particles suspended in a turbulent medium. By following the methodology of Voelk et al. (1980), but also using results that take into account the finite inner scale of eddies at ${k}_{\eta }=1/\tilde{\eta }$,⁸ OC07 arrive at an analytic expression for the rms inertial space velocity of particles suspended in a turbulent medium

$$\begin{eqnarray}&&{v}_{{pi},t}^{2}={v}_{t}^{2}\left(1-\displaystyle \frac{{{St}}^{2}(1-{{Re}}_{t}^{-\tfrac{1}{2}})}{({St}+1)({St}+{{Re}}_{t}^{-\tfrac{1}{2}})}\right).\end{eqnarray} \tag{ 22 }$$

Relative to the gas, the velocity is:

$$\begin{eqnarray}&&{v}_{{pg},t}^{2}={v}_{t}^{2}\left(\displaystyle \frac{{{St}}^{2}(1-{{Re}}_{t}^{-\tfrac{1}{2}})}{({St}+1)({St}+{{Re}}_{t}^{-\tfrac{1}{2}})}\right).\end{eqnarray} \tag{ 23 }$$

Equation (23) was originally derived by Cuzzi & Hogan (2003) for particles with ${St}\ll 1$. The results of numerical calculations and comparison with simulations, however, lead OC07 to argue that Equations (22) and (23) hold to order unity for particles of arbitrary Stokes number. We are therefore justified in applying this equation to determine the rms turbulent velocity of arbitrary-sized accreting particles, with one caveat. We first note that, for large Stokes numbers, we have ${v}_{{pi},t}\to {v}_{t}/{{Re}}_{t}^{1/4}$. This cannot be completely correct because, in the limit ${St}\to \infty$, we expect for the particles to be so massive that they are completely unaffected by turbulence, and thus ${v}_{{pi},t}\to 0$. In addition, we note that, for ${St}\ll {{Re}}_{t}$, as long as we do not have ${St}\ll 1$, Equation (22) reduces to (21). Thus, we can use (21) to determine the velocity of larger St particles; in practice, we use (21) for ${St}\geqslant 10$. This form should connect smoothly to the more precise form given in Equation (22) as long as ${{Re}}_{t}\,\gg \,10$. Using the fiducial values given in Section 5.1, the Reynolds number of the turbulence is given by

$$\begin{eqnarray}&&{{Re}}_{t}\,=\,4.07\times {10}^{10}\,\alpha \,{\left(\displaystyle \frac{a}{\mathrm{au}}\right)}^{-1}.\end{eqnarray} \tag{ 24 }$$

Thus, unless we are at extremely large orbital separations with weak turbulence, we can model the full range of Stokes numbers in this manner with little error.

We have previously discussed how to calculate the laminar and turbulent components of the particle's velocity. However, in order to calculate the total rms velocity, we need to understand how to combine these two components. Furthermore, calculation of drag forces requires knowledge of the velocity of the particle relative to the gas, whereas calculation of ${v}_{\infty }$ will require the velocity of the particle relative to the local Keplerian velocity, because this is the velocity at which small bodies approach the core. Therefore, we also need to understand how to convert our velocities from one frame to another. The methodology necessary for our calculation is summarized here; for a more detailed derivation, see Appendix B.

We let ${\boldsymbol{\delta }}{\boldsymbol{v}}$ represent the "fluctuating" or turbulent component of the gas velocity and $\bar{{\boldsymbol{v}}}$ represent the laminar component—i.e., if we write ${\boldsymbol{v}}=\bar{{\boldsymbol{v}}}+{\boldsymbol{\delta }}{\boldsymbol{v}}$, then $\langle {\boldsymbol{v}}\rangle =\bar{{\boldsymbol{v}}}$, where $\langle \ldots \rangle$ denotes ensemble averaging. We can then show that

$$\begin{eqnarray}&&{v}^{2}={\bar{v}}^{2}+\langle \delta {v}^{2}\rangle ,\end{eqnarray} \tag{ 25 }$$

i.e., the turbulent and laminar components combine in quadrature.

If the subscripts $p,g,$ and k denote the velocity of the small bodies, the gas, and the local Keplerian velocity, respectively, then

$$\begin{eqnarray}&&{\bar{{\boldsymbol{v}}}}_{{pk}}={\bar{{\boldsymbol{v}}}}_{{pg}}+{\bar{{\boldsymbol{v}}}}_{{gk}}.\end{eqnarray} \tag{ 26 }$$

Thus, the laminar component of the particle's velocity can be changed from one frame to another in the usual manner (i.e., the $\hat{{\boldsymbol{\phi }}}$ component is altered by $\eta {v}_{k}$ while the $\hat{{\boldsymbol{r}}}$ component is unchanged), independent of the turbulent velocity.

For the turbulent component of the velocity, one can show

$$\begin{eqnarray}&&\langle \delta {v}_{{gk}}^{2}\rangle =\langle \delta {v}_{{pk}}^{2}\rangle +\langle \delta {v}_{{pg}}^{2}\rangle ,\end{eqnarray} \tag{ 27 }$$

which is given in a number of works (e.g., Csanady 1963; Cuzzi et al. 1993, and OC07), and is usually derived by considering the Fourier components of the turbulent velocities in frequency space. An alternate derivation is given in Appendix B.

Neglecting the effects from the fact that the Keplerian velocity is not truly an inertial reference frame (see Youdin & Lithwick (2007) for a discussion of non-inertial effects), Equations (25)–(27) fully specify how to calculate all the velocities relevant to the problem from the input given by Equations (14), (15), and (22).

Because the Keplerian velocity of the disk varies as ${v}_{k}\propto {a}^{-1/2}$, bodies that are separated substantially in the radial direction will move relative to one another in the azimuthal direction, even in the absence of other effects. For our purposes, we approximate this "shear" velocity as

$$\begin{eqnarray}&&{v}_{\mathrm{shear}}={\rm{\Omega }}r,\end{eqnarray} \tag{ 28 }$$

where r is the separation between the two bodies in the radial direction. For small bodies encountering growing protoplanets, we have r ≈ R_acc. If this velocity is larger than the drift-dispersion velocity of small bodies, v_pk, then v_shear will set ${v}_{\infty }$, the velocity at which small bodies encounter cores; i.e., we set

$$\begin{eqnarray}&&{v}_{\infty }=\max ({v}_{{pk}},{v}_{\mathrm{shear}}).\end{eqnarray} \tag{ 29 }$$

We will refer to particles with ${v}_{\infty }={v}_{{pk}}$ as "drift-dispersion dominated," and particles with ${v}_{\infty }={v}_{\mathrm{shear}}$ as being "shear-dominated."

The relevant velocity for calculating the drag force remains to be determined. Note that this velocity is relative to the gas (in contrast to ${v}_{\infty }$, which is relative to the local Keplerian velocity) and cannot be set by shear, because the gas and the particles will have the same shear velocity.⁹ For particles that approach the core with velocities faster than the local circular orbit velocity about the large body, i.e., for particles with ${v}_{\infty }\gt {v}_{\mathrm{orbit}}=\sqrt{{GM}/{R}_{\mathrm{acc}}}$, the dominant effect of the gravitational force from the core will be to give the particle a small "kick" in the direction perpendicular to its motion. By using the impulse approximation, one can show that the magnitude of the resultant perpendicular velocity is of order ${v}_{\mathrm{kick}}={GM}/({R}_{\mathrm{acc}}{v}_{\infty })$ (see, e.g., Binney & Tremaine 2008). This effect is illustrated in Figure 4 (solid curve). As also shown in Figure 4, particles that approach the core with velocities such that ${v}_{\infty }\lt {v}_{\mathrm{orbit}}$ will experience a substantial change in their total velocity during their interaction with the core (dashed curve). During the encounter, the magnitude of the particle's velocity when the particle is a distance r from the core is approximately $| v| \approx \sqrt{{GM}/r}$ (see the lower right-hand panel of Figure 4). If the magnitude of the work done on the particle can be expressed as $W(r)={F}_{D}(r)r$, then for a particle in a linear drag regime, we have $W(r)\propto {r}^{1/2}$, i.e., the majority of the work done during encounter occurs when the particle is at the largest scales. If the particle is in a quadratic (Ram) drag regime, then F_D is independent of r, so the work done at all scales is approximately equal. In either regime we are therefore justified in approximating the work done during the encounter as $W=2{F}_{D}({v}_{\mathrm{enc}}){R}_{\mathrm{acc}}$, where ${v}_{\mathrm{enc}}={v}_{\mathrm{orbit}}({R}_{\mathrm{acc}})=\sqrt{{GM}/{R}_{\mathrm{acc}}}$.

Zoom In Zoom Out Reset image size

Finally, if the particle's drift velocity relative to the gas, v_pg, is larger than the velocity from the gravitational influence of the core, then v_pg will set the velocity during the encounter.

In summary, we define the relevant velocity for the work calculation, v_enc to be

$$\begin{eqnarray}{v}_{\mathrm{enc}}=\left\{\begin{array}{ll}\max ({v}_{\mathrm{orbit}},{v}_{{pg}}), & {v}_{\infty }\lt {v}_{\mathrm{orbit}}\\ \max ({v}_{\mathrm{kick}},{v}_{{pg}}), & {v}_{\infty }\gt {v}_{\mathrm{orbit}}.\end{array}\right.\end{eqnarray} \tag{ 30 }$$

The capture radius, R_acc, is the radius interior to which a small body will accrete if certain energy criteria are met, as discussed in Section 2.1. The scale R_acc is also used to determine both the incoming kinetic energy of a small body and the work done by gas drag during the encounter. In this work, we assume that particles that cannot dissipate their kinetic energy at R_acc will not be able to accrete on a smaller length scale. To demonstrate this is generally the case, we note that, for a general impact parameter b, we have

$$\begin{eqnarray}&&\displaystyle \frac{{KE}}{W}=\displaystyle \frac{{{mv}}_{\infty }^{2}}{4{F}_{D}({v}_{\mathrm{enc}})b}.\end{eqnarray} \tag{ 31 }$$

If the small body's velocity is shear-dominated, i.e., ${v}_{\infty }\,={v}_{\mathrm{shear}}$, then it can be shown analytically that the only particles that may have ${KE}/W\gt 1$ are those shearing into R_H. In this case, particles with substantially smaller impact parameters will not penetrate into R_H due to the nature of three-body trajectories (see, e.g., Petit & Henon 1986). If instead ${v}_{\infty }={v}_{{pk}}$, then we have

$$\begin{eqnarray}\displaystyle \frac{{KE}}{W}\propto \displaystyle \frac{1}{{F}_{D}({v}_{\mathrm{enc}})b}\propto \left\{\begin{array}{ll}{({v}_{\mathrm{enc}}b)}^{-1}, & \mathrm{Epstein},\mathrm{Stokes}\\ {({v}_{\mathrm{enc}}^{2}b)}^{-1}, & \mathrm{Ram}.\end{array}\right.\end{eqnarray} \tag{ 32 }$$

The maximal possible scaling of v_enc with impact parameter is ${v}_{\mathrm{enc}}={v}_{\mathrm{kick}}\propto {b}^{-1}$. Therefore, Equation (32) implies that KE/W is always minimized at large b when gas drag is in the Epstein or Stokes regime, and our calculation of the energy criterion at the largest possible accretion scale is sufficient. The only case for which gas-assisted growth may be ruled out by the energy criterion at a large scale and yet possible at a smaller impact parameter is in the Ram pressure drag regime. We do not include this possibility in our calculation of the accretion rate. This case only holds over a small amount of parameter space, as it requires particles to have both ${r}_{s}\gt 9\lambda /4$ and Re ≫ 1. Thus, our modeling may rule out accretion of large particles in the inner disk (where the Ram regime is most important) that would be able to accrete at smaller scales.

To determine ${R}_{\mathrm{acc}}$, we note that capture is, in principle, possible either within the planet's atmosphere, where the gas density increases substantially, or within the radius at which a small body can stably orbit the core. Thus, the capture radius is set by the relative size of the atmospheric radius, which extends up to R_b (for ${R}_{b}\lt {R}_{H};$ see Section 7 for ${R}_{b}\gt {R}_{H}$), and the stability radius, R_stab:

$$\begin{eqnarray}&&{R}_{\mathrm{acc}}=\max ({R}_{\mathrm{stab}},{R}_{b}).\end{eqnarray} \tag{ 33 }$$

The stability radius is determined by requiring that the force on the small body is dominated by the gravity of the core. We also consider two other forces that can disrupt accretion, which leads to two different length scales that can set R_stab.

4.1.1. The Hill Radius

The first scale is set by demanding that the small body not be sheared off by the gravity from the host star; this leads to the traditional measure of stability, the Hill radius, where the gravitational acceleration from the core is balanced by the tidal gravity from the star (Hill 1878). For ${M}_{* }\gg M$, where M_* is the mass of the central star,

$$\begin{eqnarray}&&{R}_{H}\approx a{\left(\displaystyle \frac{M}{3{M}_{* }}\right)}^{1/3},\end{eqnarray} \tag{ 34 }$$

where a is the semimajor axis of the core's orbit.

The fact that accretion can occur for impact parameters of order R_H is one of the main enhancements to the growth rate that come from pebble accretion. In the absence of gas, the maximal impact parameter for accretion is smaller than R_H by a factor ${(R/{R}_{H})}^{1/2}$, where R is the radius of the planet (see Appendix C). Both Lambrechts & Johansen (2012, hereafter LJ12) and Ormel & Klahr (2010, hereafter OK10) find regimes where ${R}_{\mathrm{acc}}\sim {R}_{H}$ in their numerical calculations of pebble accretion. In LJ12's framework, accretion at R_H occurs for core masses greater than ${M}_{t}=\sqrt{1/3}{v}_{\mathrm{gas}}^{3}/(G{\rm{\Omega }})$ and particle sizes St ≳ 1.¹⁰

In OK10's framework, impact parameters comparable to R_H occur in what they refer to as the "three-body regime," which occurs for St > 1, and ${St}\gt {\zeta }_{w}\equiv \eta {v}_{k}/{v}_{H}$, where ${v}_{H}\equiv {R}_{H}{\rm{\Omega }}$ is the Hill velocity. In this regime, OK10 give the impact parameter as ${b}_{3b}/{R}_{H}=1.7{p}^{1/2}+1.0/{St}$, where $p\equiv R/{R}_{H}$. This agrees generally with our results that ${R}_{\mathrm{acc}}\approx {R}_{H}$ in the Hill accretion regime.

The $1/{St}$ dependence given by OK10 appears to encapsulate the chaotic nature of trajectories for particles that shear into the Hill sphere: often particles will orbit the core many times before exiting the Hill sphere. This effect can cause particles that would not accrete according to our energy criterion to be captured, because our modeling only includes particles that dissipate their energy over one orbital crossing. If we use the result from Goldreich et al. (2002), who studied binary capture of Kuiper Belt objects by dynamical friction, that the probability of capture, P_cap, for a particle that shears into the Hill radius is equal to the fraction of energy it dissipates over one orbit, then this leads to an effective accretion rate of

$$\begin{eqnarray}&&\displaystyle \frac{{dM}}{{dt}}\sim {R}_{H}{v}_{H}{{\rm{\Sigma }}}_{p}\displaystyle \frac{W}{{KE}}.\end{eqnarray} \tag{ 35 }$$

For high-mass cores, the gravitational force from the core dominates the velocity of the small body during the encounter, and the encounter velocity is ${v}_{\mathrm{enc}}\sim {v}_{H}$. This in turn implies that $W/{KE}\propto 1/{St}$. For 2D accretion (${H}_{p}\lt {R}_{\mathrm{acc}}$), which is the regime modeled by OK10, our expression therefore agrees with their results. For low-mass cores, however, the particle-gas relative velocity will instead dominate the velocity during the encounter, which will lead to scaling that differs from $1/{St}$ in this regime.

In summary, particles that have ${R}_{\mathrm{stab}}={R}_{H}$ and ${v}_{\infty }={v}_{H}$ accrete with timescales given by

$$\begin{eqnarray}&&{t}_{\mathrm{grow}}^{{\prime} }=\displaystyle \frac{{t}_{\mathrm{grow}}}{\min (1,W/{KE})},\end{eqnarray} \tag{ 36 }$$

where ${t}_{\mathrm{grow}}$ is given by Equation (6). In what follows, Equation (36) is used to modify t_grow for particles that shear into R_H.

4.1.2. The WISH Radius and Shearing Radius

In the presence of gas, even if the gravity from the core dominates over the tidal gravity from the star, it is possible that the relative acceleration of the core and the small body will be strong enough to strip the small body away from the core. In this regime, ${R}_{\mathrm{stab}}$ is set not by the Hill radius, but by what Perets & Murray-Clay refer to as the wind-shearing (WISH) radius (Perets & Murray-Clay 2011). At the WISH radius, the differential acceleration between the large and the small body due to gas drag is balanced by gravitational acceleration. That is,

$$\begin{eqnarray}&&{R}_{\mathrm{WS}}^{{\prime} }=\sqrt{\displaystyle \frac{G(M+m)}{{\rm{\Delta }}{a}_{\mathrm{WS}}}},\end{eqnarray} \tag{ 37 }$$

where m is the mass of the small body, and ${\rm{\Delta }}{a}_{\mathrm{WS}}$ is the differential acceleration between the small and the large body, given by

$$\begin{eqnarray}&&{\rm{\Delta }}{a}_{\mathrm{WS}}=\left|\displaystyle \frac{{F}_{D}(M)}{M}-\displaystyle \frac{{F}_{D}(m)}{m}\right|,\end{eqnarray} \tag{ 38 }$$

where F_D is the force exerted on the particle due to gas drag. For the mass ratios considered here, the first term in Equation (38) is generally negligible. To determine the relevant velocity for calculating the drag force F_D, we note that particles that are accreted by the core will have their velocity modified by the core's gravity, increasing their velocity relative to the gas during the encounter. The most extreme case occurs when the velocity of particles match the local Keplerian velocity at which the core is moving. These particles temporarily experience the full gas velocity with respect to the core, which is a combination of the sub-Keplerian and turbulent velocities, as well as the Keplerian shear in the disk. We can encompass this behavior in the limit $M\gg m$ by writing ${R}_{\mathrm{WS}}^{{\prime} }$ as

$$\begin{eqnarray}&&{R}_{\mathrm{WS}}^{{\prime} }\approx \sqrt{\displaystyle \frac{{GMm}}{{F}_{d}(\max [{v}_{\mathrm{gas}},{v}_{\mathrm{shear}}])}},\end{eqnarray} \tag{ 39 }$$

where ${v}_{\mathrm{gas}}=\sqrt{{\eta }^{2}{v}_{k}^{2}+\alpha {c}_{s}^{2}}$ is the rms velocity of the nebular gas; we use the prime to differentiate this radius from the definition of R_WS used below. The case ${v}_{\mathrm{gas}}\gt {v}_{\mathrm{shear}}$ is discussed in detail by Perets & Murray-Clay (2011). In the three drag regimes discussed in Section 3.1, Perets & Murray-Clay (2011) give R_WS as

$$\begin{eqnarray}{R}_{\mathrm{WS}}={({GM}{\rho }_{s})}^{1/2}\times \left\{\begin{array}{ll}{\left(\displaystyle \frac{1}{{\rho }_{g}{v}_{\mathrm{th}}{v}_{\mathrm{gas}}}\right)}^{1/2}{r}_{s}^{1/2}, & \mathrm{Epstein}\\ {\left(\displaystyle \frac{4}{9}\displaystyle \frac{\sigma }{\mu {v}_{\mathrm{th}}{v}_{\mathrm{gas}}}\right)}^{1/2}{r}_{s}, & \mathrm{Stokes}\\ {\left(\displaystyle \frac{1}{0.165}\displaystyle \frac{1}{{\rho }_{g}{v}_{\mathrm{gas}}^{2}}\right)}^{1/2}{r}_{s}^{1/2}, & \mathrm{Ram}\end{array}\right.\end{eqnarray} \tag{ 40 }$$

where ρ_s is the internal density of the small body. We will use the term WISH radius and the symbol R_WS to refer solely to the case ${v}_{\mathrm{gas}}\gt {v}_{\mathrm{shear}}$. For the case ${v}_{\mathrm{shear}}\gt {v}_{\mathrm{gas}}$, on the other hand, the right-hand side of Equation (39) now depends on the impact parameter. In this regime, we will refer to the impact parameter as the "shearing radius," R_shear. In general, R_shear is determined by numerically solving the equation

$$\begin{eqnarray}&&{R}_{\mathrm{shear}}=\sqrt{\displaystyle \frac{{GMm}}{{F}_{d}({R}_{\mathrm{shear}}{\rm{\Omega }})}}.\end{eqnarray} \tag{ 41 }$$

For a particle in a linear drag regime, the particle's Stokes number is independent of velocity, which allows us to solve analytically for R_shear:

$$\begin{eqnarray}&&{R}_{\mathrm{shear}}={(3{St})}^{1/3}{R}_{H}.\end{eqnarray} \tag{ 42 }$$

In Section 5.5, we will compare our results to OK10 and LJ12 in the laminar regime—for now, we translate their impact parameters into our notation.

The WISH radius is equivalent, in the laminar regime, to the effective accretion radius, r_d, used by LJ12, as well as to the settling radius, b_set, used by OK10. LJ12 give the effective accretion radius as ${r}_{d}={({t}_{B}/{t}_{s})}^{-1/2}{r}_{B}$. Here, ${t}_{B}={r}_{B}/{v}_{\infty }$ is the drift time across the "Bondi radius"—${r}_{B}={GM}/{v}_{\infty }^{2}$.¹¹ As LJ12 note, r_d is equal to the WISH radius, as can be seen from that fact that

$$\begin{eqnarray}&&{r}_{d}=\sqrt{{t}_{s}\eta {v}_{k}\displaystyle \frac{{GM}}{{(\eta {v}_{k})}^{2}}}={R}_{\mathrm{WS}}.\end{eqnarray} \tag{ 43 }$$

LJ12 also give a radius equal to R_shear (up to a factor of 3^1/3) for particles with ${St}={10}^{-2}$ and $M\gt {M}_{t}$.

In their "settling" regime, OK10 determine the impact parameter ${b}_{\mathrm{set}}={R}_{\mathrm{acc}}/{R}_{H}$ by solving the cubic equation

$$\begin{eqnarray}&&{b}_{\mathrm{set}}^{3}+\displaystyle \frac{2{\zeta }_{w}}{3}{b}_{\mathrm{set}}^{2}-8{St}=0.\end{eqnarray} \tag{ 44 }$$

If the cubic term in the above equation is negligible, then we get the solution

$$\begin{eqnarray}&&{b}_{\mathrm{set}}=\sqrt{12{St}/{\zeta }_{w}}.\end{eqnarray} \tag{ 45 }$$

In terms of ζ_w, we can rewrite R_WS as

$$\begin{eqnarray}&&{R}_{\mathrm{WS}}=\sqrt{3}{\left(\displaystyle \frac{{St}}{{\zeta }_{w}}\right)}^{1/2}{R}_{H},\end{eqnarray} \tag{ 46 }$$

so ${b}_{\mathrm{set}}=2{R}_{\mathrm{WS}}$ in the laminar regime. If, on the other hand, the middle term is negligible, the solution is

$$\begin{eqnarray}&&{b}_{\mathrm{set}}=2{{St}}^{1/3},\end{eqnarray} \tag{ 47 }$$

which is $\sim {R}_{\mathrm{shear}}$. One can verify by inspection that the solution to Equation (44) is approximately the minimum of Equations (45) and (47).

4.1.3. Combining Length Scales

Once R_H, R_WS, and R_shear are calculated, we take the stability radius to simply be the smallest of these three radii:

$$\begin{eqnarray}&&{R}_{\mathrm{stab}}=\min ({R}_{\mathrm{WS}},{R}_{\mathrm{shear}},{R}_{H}).\end{eqnarray} \tag{ 48 }$$

We now have the pieces necessary to determine R_acc. We first compute R_WS and R_shear; the minimum of these two length scales determines the scale on which gas drag will pull small bodies off the core before they can be accreted. We then calculate R_H, the scale on which the stellar gravity disrupts accretion, and calculate R_stab by taking in the minimum of these scales. Finally, we take R_acc to be:

$$\begin{eqnarray}&&{R}_{\mathrm{acc}}=\max ({R}_{\mathrm{stab}},{R}_{b}).\end{eqnarray} \tag{ 49 }$$

Note that accretion within R_acc only occurs when additional energy criteria are met, as discussed in Section 2.1.

We now turn to the height of the accretion rectangle. The height of this rectangle will be the minimum of the dust scale height and the capture radius, as any particle outside the capture radius will not be able to accrete on to the core:

$$\begin{eqnarray}&&{H}_{\mathrm{acc}}=\min ({R}_{\mathrm{acc}},{H}_{p}).\end{eqnarray} \tag{ 50 }$$

In general, the solid particles will tend to settle to the midplane of the disk, due to the vertical influence of gravity and the lack of pressure support. However, several processes will tend to oppose this settling, giving the particles some finite vertical extent. In our model, we include two different processes that drive particles vertically—the first is turbulent diffusion due to the isotropy of the small body's turbulent velocities.

For St > α, the scale height of the particles due to their interaction with turbulence is given by (Youdin & Lithwick 2007):

$$\begin{eqnarray}&&{H}_{t}=\sqrt{\displaystyle \frac{\alpha }{{St}}}{H}_{g}.\end{eqnarray} \tag{ 51 }$$

We can understand this expression by separately considering larger particles with St ≫ 1 and small, well-coupled particles with St ≪ 1. For the larger particles, as discussed in Section 3.3, we can think of the particle's velocity as resulting from the large number of uncorrelated kicks it receives from the eddies. In this case, we have ${v}_{{pk},t}\sim {v}_{t}/\sqrt{{St}}$, which gives a scale height of

$$\begin{eqnarray}&&{H}_{t}\sim \displaystyle \frac{{v}_{t}}{{\rm{\Omega }}\sqrt{{St}}}\sim \sqrt{\displaystyle \frac{\alpha }{{St}}}{H}_{g},\end{eqnarray} \tag{ 52 }$$

since ${v}_{t}=\sqrt{\alpha }{H}_{g}{\rm{\Omega }}$.

As described in Youdin & Lithwick (2007), for St ≪ 1 we expect the particles to be well-coupled to the gas, and thus expect the particles' diffusion coefficient to be approximately equal to the gas value. We can express the turbulent diffusion coefficient of the gas as ${D}_{g}={\alpha }_{z}{c}_{s}{H}_{g}$. In what follows, we take α_z ≈ α unless otherwise stated. While α and α_z can differ, because α enters our calculation mainly as a parameterization of the turbulent gas velocity, we are more interested in the relative size of α_z and v_z, the vertical turbulent gas velocity. Magnetohydrodynamic simulations find that these two values are similar to order of magnitude (e.g., Xu et al. 2017). Tightly coupled particles should settle to the midplane at approximately their terminal velocity, so their time to settle to the midplane from a height z is approximately ${t}_{\mathrm{set}}\sim z/{v}_{\mathrm{term}}$. At terminal velocity, the drag force balances the vertical component of gravity, ${F}_{g,z}\sim ({GMm}/{r}^{2})(z/r)\sim m{{\rm{\Omega }}}^{2}z$. Setting ${F}_{d}\sim m{{\rm{\Omega }}}^{2}z$ implies that ${v}_{\mathrm{term}}\sim {St}{\rm{\Omega }}z$, which gives ${t}_{\mathrm{set}}\sim 1/({St}\,{\rm{\Omega }})$. Finally, setting ${t}_{\mathrm{set}}$ equal to the particle diffusion time ${H}_{t}^{2}/{D}_{p}$ gives Equation (51). Carballido et al. (2006) derive (51) for St ≫ 1 by more rigorous means, whereas Dubrulle et al. (1995) do the same for St ≪ 1.

For St > α, we have H_t > H_g, which is clearly incorrect; turbulence cannot drive particles to heights larger than the extent of the gas disk. We therefore cap H_t at H_g. Dubrulle et al. (1995) derive (51) as the scale height of the dust-to-gas ratio ${\rho }_{s}/{\rho }_{g}$ and note that this height h is related to the dust scale height by the equation

$$\begin{eqnarray}&&{H}_{t}=h{\left[1+{\left(\displaystyle \frac{h}{{H}_{g}}\right)}^{2}\right]}^{-1/2},\end{eqnarray} \tag{ 53 }$$

which has the expected behavior that ${H}_{t}\to {H}_{g}$ as ${St}\to \infty$. Equation (53) can be used in the place of Equation (51), though this does not substantially change the results. We employ (51) in order to maintain relative simplicity in our analytic results.

Turbulence's modification to the particle scale height is discussed in a number of works, including OK10 and LJ12, though it is often not explicitly included in calculations of growth rate. Instead, authors generally note that, for ${H}_{t}\gt {R}_{\mathrm{acc}}$, the growth timescale is increased by a factor of ${H}_{t}/{R}_{\mathrm{acc}}$. We also note that Equation (53) can be written as

$$\begin{eqnarray}&&{H}_{t}={H}_{g}\sqrt{\displaystyle \frac{\alpha }{\alpha +{St}}},\end{eqnarray} \tag{ 54 }$$

which is employed by Ormel & Kobayashi (2012) for the dust scale height in the presence of turbulence. Equation (54) agrees with Equation (28) of Youdin & Lithwick (2007) without their factor of ${\xi }^{-1/2}$, where $\xi \equiv 1+{St}/(1+{St})$ for the values chosen in this paper.

Even in the absence of strong turbulence, the Kelvin–Helmholtz shear instability prevents small bodies from settling too close to the midplane. For small particles, which are well-coupled to the gas flow, this leads to a scale height of

$$\begin{eqnarray}&&{H}_{\mathrm{KH}}^{{\prime} }=\displaystyle \frac{{H}_{g}^{2}}{a}=\displaystyle \frac{2\eta {v}_{k}}{{\rm{\Omega }}},\end{eqnarray} \tag{ 55 }$$

(e.g., Lee et al. 2010). In order to extend this scale height to include larger particles, we assume that, even when large bodies dominate the mass distribution, a population of small bodies exists that is sufficient to drive the Kelvin–Helmholtz shear instability—and therefore turbulence—close to the midplane. We expect the turbulent velocity in this case to be comparable to the vertical shear rate, which is ∼ηv_k. In this case, we can make arguments analogous to those preceding Equation (52), leading to a scale height ${H}_{p}\sim \eta {v}_{k}/({\rm{\Omega }}\sqrt{{St}})$. Thus, we can describe H_KH over the full range of particle sizes using the expression

$$\begin{eqnarray}&&{H}_{\mathrm{KH}}=\displaystyle \frac{{H}_{g}^{2}}{a}\min (1,{{St}}^{-1/2}).\end{eqnarray} \tag{ 56 }$$

If only large bodies are present, they may need to settle further in order to excite the Kelvin–Helmholtz instability, which would change the dependence on St. We leave a self-consistent calculation to another work. Furthermore, by employing this expression for larger particle sizes, we are neglecting the possibility of dynamical stirring of the large bodies. Mutual scatterings and/or the gravitational force of the core can drive small particles vertically, and may be more important than interactions with the gas for determining the scale height of larger particles, which are decoupled from the gas.

Finally, we take the solid scale height to simply be the maximum of these two heights:

$$\begin{eqnarray}&&{H}_{p}=\max ({H}_{t},{H}_{\mathrm{KH}}).\end{eqnarray} \tag{ 57 }$$

The possibility of 3D accretion (${H}_{p}\gt {R}_{\mathrm{acc}}$) in the non-turbulent regime is generally neglected in works on pebble accretion, but we find that it has non-trivial effects on the growth timescale.

We begin by providing the values of the fiducial parameters used for describing the protoplanetary disk, which are necessary to provide concrete numerical results. A discussion of the effects of varying some of these parameters is given in Section 6.

We assume a surface density profile of

$$\begin{eqnarray}&&{\rm{\Sigma }}(a)=500\,{\rm{g}}\,{\mathrm{cm}}^{-2}\,{\left(\displaystyle \frac{a}{\mathrm{au}}\right)}^{-1},\end{eqnarray} \tag{ 58 }$$

where a is the semimajor axis in the disk. This profile is chosen to be consistent with measurements of the surface density in solids of size 0.1–1 mm, taken from sub-mm observations of protoplanetary disks (Andrews et al. 2009; Andrews 2015), which probe these particle sizes. In order to calculate the surface density profile of solid bodies, Σ_p, we assume a constant dust-to-gas ratio of ${f}_{s}\equiv {{\rm{\Sigma }}}_{p}/{\rm{\Sigma }}=1/100$. The semimajor axis dependence of the temperature profile is taken from Chiang & Youdin (2010)

$$\begin{eqnarray}&&T={T}_{0}{\left(\displaystyle \frac{a}{\mathrm{au}}\right)}^{-3/7}.\end{eqnarray} \tag{ 59 }$$

For most purposes, we will take the prefactor to be ${T}_{0}=200\,{\rm{K}}$, which would place the Earth outside of the water ice line during its formation (e.g., Powell et al. 2017). We also note this prefactor is consistent with the assumption that the disk is heated primarily by irradiation from a central star with luminosity ${L}_{* }\sim 3\,{L}_{\odot }$ (e.g., Ida et al. 2016), which is appropriate for a solar mass star of age 1 Myr (Tognelli et al. 2011). The effects of varying T₀ are discussed in Section 6.2.

From here, we can then calculate the isothermal sound speed and scale height of the disk in the usual manner: ${c}_{s}=\sqrt{{kT}/\mu }$ and ${H}_{g}={c}_{s}/{\rm{\Omega }}$. Here, μ is the mean molecular weight of particles in the disk; we take $\mu =2.35{m}_{H}\approx 3.93\times {10}^{-24}$ g, which assumes a disk composed of 70% ${{\rm{H}}}_{2}$ and 30% He by mass. We can also calculate the gas density ${\rho }_{g}={\rm{\Sigma }}/(2{H}_{g})$, and the mean free path of gas particles, $\lambda =\mu /({\rho }_{g}\sigma )$, where σ is the neutral collision cross section, which we take to be $\sigma \sim {(3\mathring{\rm A} )}^{2}\sim {10}^{-15}\,{\mathrm{cm}}^{2}$. In order to convert from Stokes number, which is the relevant parameter for calculating ${t}_{\mathrm{grow}}$, to physical size, we need to specify the internal density of the pebbles, ρ_s. We take this density to be ${\rho }_{s}=2\,{\rm{g}}\,{\mathrm{cm}}^{-3}$ unless stated otherwise, which is appropriate for rocky bodies. Note, however, that the solids in protoplanetary disks may be fluffy, i.e., their densities may be low, which can have important ramifications for planet formation (see, e.g., Kataoka et al. 2013). While we consider the effect of varying M_* in Section 6.4, unless otherwise noted, all values quoted in the paper are for a solar mass star, ${M}_{* }={M}_{\odot }$. As noted in Section 2.1, the radius of the large body needs to be ≳10 km, which corresponds to M ∼ 10⁻⁹ M_⊕ for a density of 2 g cm⁻³. In practice, we do not consider cores that are close to this lower limit.

A summary of our fiducial disk parameters is given in Table 2.

Table 2.  Fiducial Disk Parameters

Parameter	Symbol	Formula/Valuea
Solid-to-gas ratio	fs	0.01
Mean molecular weight of gas particles	μ	$2.35\,{m}_{H}\approx 3.93\times {10}^{-24}\,{\rm{g}}$
Neutral collision cross section	σ	${10}^{-15}\,{\mathrm{cm}}^{2}$
Gas surface density	Σ	${\rm{\Sigma }}(a)=500{\left(\tfrac{a}{\mathrm{au}}\right)}^{-1}\,{\rm{g}}\ {\mathrm{cm}}^{-2}$
Protoplanetary disk temperature	T	$T(a)=200{\left(\tfrac{a}{\mathrm{au}}\right)}^{-3/7}\,{\rm{K}}$
Isothermal sound speed	cs	${c}_{s}=\sqrt{\tfrac{{kT}}{\mu }}$
Average thermal velocity	vth	${v}_{\mathrm{th}}=\sqrt{\tfrac{8}{\pi }}{c}_{s}$
Gas scale height	Hg	${H}_{g}=\tfrac{{c}_{s}}{{\rm{\Omega }}}$
Gas density	${\rho }_{g}$	${\rho }_{g}=\tfrac{{\rm{\Sigma }}}{2{H}_{g}}$
Gas mean free path	λ	$\lambda =\tfrac{\mu }{{\rho }_{g}\sigma }$

Note.

^aA value indicates that the quantity is a constant input parameter to the code, whereas a formula indicates that the quantity is calculated in the prescribed manner from input parameters and constants.

Download table as:  ASCIITypeset image

For the purpose of understanding the broad features of pebble accretion, we present a typical output of our model in Figure 5. Figure 5 plots the growth timescale (Equation (6)) as a function of the small body radius, r_s, for a core of mass M = 10⁻¹ M_⊕ at a = 1 au. The growth timescales for several α values, ranging from purely laminar flow to extremely strong turbulence, are shown. The timescale is set to $\infty$ for particles that do not satisfy the energy criteria discussed in Section 2.1, so the plotted range for each value of α indicates the range of small body radii the core can accrete via pebble accretion. The Stokes numbers corresponding to the given values of radius are shown for reference (the laminar drift velocity is used to calculate St for particles not in a linear drag regime).

Zoom In Zoom Out Reset image size

We begin by examining the laminar case, where α = 0. For large values of r_s, particles shear into the Hill radius. However, growth is slow because particles are large and cannot fully dissipate their kinetic energy relative to the gas, which increases t_grow by a factor of KE/W. Once particles are small enough that KE ∼ W, they can accrete on rapid timescales. In terms of the small body's Stokes number, we can generally write the criterion for particles to be able to deplete their kinetic energy as

$$\begin{eqnarray}&&{St}\lt \displaystyle \frac{4{v}_{\mathrm{enc}}{R}_{\mathrm{acc}}{\rm{\Omega }}}{{v}_{\infty }^{2}}.\end{eqnarray} \tag{ 60 }$$

For particles of this size, gravity dominates over gas drag effects—therefore, R_H determines R_acc, as opposed to R_WS. The large size of these particles also means that their velocity dispersion and drift velocity are low, so the approach velocity is set by shear, i.e., ${v}_{\infty }={v}_{H}$. Using these values to calculate v_enc and plugging the result into Equation (60) gives

$$\begin{eqnarray}&&{St}\lt 4\sqrt{3},\end{eqnarray} \tag{ 61 }$$

which will be a general upper limit on Stokes numbers that the core is able to accrete rapidly. Note that, if the small body is not in a linear drag regime, then the Stokes number is no longer independent of velocity. In this case, the Stokes number calculated for the small body moving through the disk (which is relative to v_∞) will not necessarily be the same as the Stokes number during the encounter with the core (which is relative to v_enc).

Growth for R_acc = R_H and ${v}_{\infty }={v}_{H}$ is extremely rapid, in comparison to the planetesimal accretion timescale. As small body size continues to decrease, however, gas drag considerations eventually overtake three-body effects, and the shear radius shrinks below the Hill radius. This causes the growth timescale to increase for ${r}_{s}\lesssim 20\,\mathrm{cm}$. From this point, the shear radius continues to shrink with small body size, slowing accretion. The first kink in the slope around r_s ∼ 15 cm comes from small bodies changing from the fluid regime to the diffuse regime, while the second kink around r_s ∼ 0.2 cm stems from R_shear shrinking below H_p, which dilutes the density of small bodies, in turn slowing accretion. Eventually, the shear radius shrinks below the Bondi radius, which signals the point at which particles are so small that they will couple to the gas and flow around R_b. This causes growth to cut off for small values of r_s.

As increasingly strong turbulence is taken into account, the picture becomes more complicated. For the α = 10⁻¹ case, growth for large particle sizes is quite similar to the laminar case, because the growth parameters are solely determined by the core's mass for high St. As the particle size shrinks, we still find a small range of radii where growth occurs at R_acc = R_H, ${v}_{\infty }={v}_{H}$. Rather rapidly however, several new features emerge that were not previously present. First, the strong turbulence greatly increases the particle scale height. As particle size decreases, this scale height increases until it becomes larger than R_H, slowing down growth. As particle size decreases, for St < 1, the drift-dispersion velocity of particles increases as they couple more strongly to the gas. Due to the rapid turbulent velocity of the gas, v_pk actually overtakes v_H around ${r}_{s}\,=100\,\mathrm{cm}$, slightly mitigating the increase in growth timescale. The WISH radius is also decreased due to the strong turbulence, which makes it easier to pull particles off the core. Because of this, the radius at which ${R}_{\mathrm{WS}}\lt {R}_{H}$ is much higher than in the laminar case, and growth is again inhibited around ${r}_{s}=50\,\mathrm{cm}$ as the WISH radius begins determining ${R}_{\mathrm{acc}}$. Another kink in the slope appears around ${r}_{s}=20\,\mathrm{cm}$ as the scale height of particles becomes so large that H_p = H_g, at which point the particle scale height stops increasing. The final kink occurs at r_s ≈ 15 cm, and is again due to the particles switching from the fluid to the diffuse regime. Growth again stops when ${R}_{\mathrm{WS}}\lt {R}_{b}$, but because of the decreased values of WISH radius relative to the laminar case, this occurs at a larger value of r_s. Similar features occur for the smaller values of α, but due to the weaker turbulence, these features occur at smaller values of r_s. The only exception is the domination of the turbulent rms velocity over v_H, which only occurs in the α = 10⁻¹ case for these values of parameters. As can be seen in the figure, for lower values of α, wider ranges of small body sizes can accrete at what appears to be a minimal value of timescale, which we discuss in the next section.

The floor of the growth timescale in Figure 5 corresponds to a minimum value of t_grow, which cores reach when they can accrete over their entire Hill sphere. The appearance of this minimum growth timescale, or maximal growth rate, is a consistent thread throughout the parameter space probed by our model.

We refer to this timescale as the "Hill timescale," or t_Hill. This timescale is reached when all aspects of the accretion processes are set by the core's gravity, i.e., when we have ${R}_{\mathrm{acc}}={R}_{H}\gt {H}_{p}$ and ${v}_{\infty }={v}_{H}$. From Equations (6) and (50), it is easy to see that, if ${H}_{p}\lt {R}_{\mathrm{acc}}$, then ${t}_{\mathrm{grow}}$ is independent of H_p. Using these values for our timescale parameters in (6) gives the value of ${t}_{\mathrm{Hill}}$ as

$$\begin{eqnarray}&&{t}_{\mathrm{Hill}}=\displaystyle \frac{M}{2{f}_{s}{\rm{\Sigma }}{R}_{H}^{2}{\rm{\Omega }}}.\end{eqnarray} \tag{ 62 }$$

A similar regime is identified in Section 3.2 of Lambrechts & Johansen (2012); they refer to it as "Hill Accretion." However, they only use this term to differentiate between the drift-dominated and shear-dominated regimes, which roughly corresponds to where ${R}_{\mathrm{acc}}={R}_{H}$ for their model. Within their Hill regime, the growth timescale is not necessarily equal to ${t}_{\mathrm{Hill}}$. See Section 5.5 for a more in-depth comparison of the two models.

In general, ${t}_{\mathrm{Hill}}$ represents an extremely rapid timescale for growth. Scaled to fiducial values, the minimum timescale can be expressed as

$$\begin{eqnarray}&&{t}_{\mathrm{Hill}}\approx 4300\,\mathrm{years}\,{\left(\displaystyle \frac{a}{\mathrm{au}}\right)}^{1/2}{\left(\displaystyle \frac{M}{{M}_{\oplus }}\right)}^{1/3}.\end{eqnarray} \tag{ 63 }$$

This represents a significant enhancement over the timescale achievable via canonical core accretion, as can be seen by comparing Equations (63) and (130).

This is the substantial decrease in growth timescale promised by gas-assisted growth. The source of this decrease can be mostly attributed to an increase in ${R}_{\mathrm{acc}}$. As previously stated, the minimum timescale indicates accretion in a regime such that essentially all particles that encounter the Hill radius of the growing planet are accreted. By contrast, in the absence of gas, the trajectories of particles that enter the Hill radius are chaotic, and whether a collision occurs is extremely sensitive to the value of the small body's impact parameter; see, for example, Petit & Henon (1986). For this reason, encounters with the core in this regime are often treated probabilistically, with the collision rate of small bodies equal to the product of the rate at which small bodies encounter the radius and the probability that a particle inside the Hill radius will accrete. Performing such a calculation leads to an effective impact parameter for accretion in the so-called "three-body regime" that is, to order of magnitude, equal to the geometric mean of the Hill radius and the planetary radius—$b\sim \sqrt{{{RR}}_{H}}$; see GLS and Ormel & Klahr (2010), Equation (5). From comparison with Equation (129), we see that the maximum increase in accretion rate that can be provided by pebble accretion is of order ${R}_{H}/R$.

To see why ${t}_{\mathrm{Hill}}$ is generally the fastest rate of growth possible, we first note that because ${R}_{\mathrm{stab}}\,=\min ({R}_{H},{R}_{\mathrm{WS}},{R}_{\mathrm{shear}})$, it is clear that the maximal accretion cross section is of order ${\sigma }_{\mathrm{gas},\max }\sim {R}_{H}^{2}$. Given that turbulence can increase the rms velocity of small bodies to values substantially larger than v_H, however, it is conceivable that strong turbulence could increase the growth rate by increasing ${v}_{\infty }$, i.e., by increasing the rate at which small bodies encounter the core. However, increasing the velocity of small bodies drives them vertically as well as horizontally, which decreases the density of small bodies and slows accretion. This trade-off rules out accretion on timescales faster than ${t}_{\mathrm{Hill}}$, which we demonstrate below, to order of magnitude. We first note that

$$\begin{eqnarray}&&{H}_{\mathrm{KH}}=\displaystyle \frac{2\eta {v}_{k}}{{\rm{\Omega }}}\max (1,{{St}}^{-1/2})\gtrsim \displaystyle \frac{{v}_{{pk},{\ell }}}{{\rm{\Omega }}},\end{eqnarray} \tag{ 64 }$$

$$\begin{eqnarray}&&{H}_{t}\approx \displaystyle \frac{{v}_{{pk},t}}{{\rm{\Omega }}}.\end{eqnarray} \tag{ 65 }$$

We also make the approximation that

$$\begin{eqnarray}&&{v}_{{pk}}=\max ({v}_{{pk},{\ell }},{v}_{{pk},t}).\end{eqnarray} \tag{ 66 }$$

rather than the quadrature of ${v}_{{pk},{\ell }}$ and ${v}_{{pk},t}$. We separately consider "2D accretion," defined by ${H}_{p}\lt {R}_{\mathrm{acc}}$, and "3D accretion," where ${H}_{p}\gt {R}_{\mathrm{acc}}$.

For 2D accretion, we can write ${t}_{\mathrm{grow}}$ as

$$\begin{eqnarray}&&{t}_{\mathrm{grow}}=\displaystyle \frac{M}{2{{\rm{\Sigma }}}_{p}{R}_{H}{v}_{\infty }}\end{eqnarray} \tag{ 67 }$$

$$\begin{eqnarray}&&=\,{t}_{\mathrm{Hill}}\left(\displaystyle \frac{{R}_{H}{\rm{\Omega }}}{{v}_{\infty }}\right),\end{eqnarray} \tag{ 68 }$$

where ${t}_{\mathrm{Hill}}\equiv M/(2{{\rm{\Sigma }}}_{p}{R}_{H}^{2}{\rm{\Omega }})$ and ${{\rm{\Sigma }}}_{p}\equiv {f}_{s}{\rm{\Sigma }}$. Thus, ${t}_{\mathrm{grow}}\lt {t}_{\mathrm{Hill}}$ requires ${v}_{\infty }\gt {v}_{H}$ while still having ${R}_{H}\gt {H}_{p}$. However, we see from Equations (64)–(66) that ${R}_{H}\gt {H}_{p}$ implies ${v}_{H}\gt {v}_{{pk}}$, so ${v}_{\infty }={v}_{H}$. Thus, the timescale for growth cannot drop below ${t}_{\mathrm{Hill}}$ for 2D accretion.

We now consider the 3D accretion regime when the Hill radius sets the geometry, i.e., ${H}_{p}\gt {R}_{\mathrm{acc}}={R}_{H}$. In this case, we can write the growth timescale as

$$\begin{eqnarray}&&{t}_{\mathrm{grow}}=\displaystyle \frac{{{MH}}_{p}}{2{{\rm{\Sigma }}}_{p}{R}_{H}^{2}{v}_{\infty }}={t}_{\mathrm{Hill}}\left(\displaystyle \frac{{H}_{p}{\rm{\Omega }}}{{v}_{\infty }}\right).\end{eqnarray} \tag{ 69 }$$

If the particle is shear-dominated, then ${t}_{\mathrm{grow}}\gt {t}_{\mathrm{Hill}}$ because ${R}_{H}\lt {H}_{p}$. Furthermore, because Equations (64)–(66) imply that ${H}_{p}{\rm{\Omega }}\gtrsim {v}_{{pk}}$, if the particle is drift-dispersion dominated then we again have ${t}_{\mathrm{grow}}\gt {t}_{\mathrm{Hill}}$. Thus, in the 3D accretion regime, ${t}_{\mathrm{Hill}}$ also represents the smallest possible timescale.

While growth at ${t}_{\mathrm{Hill}}$ is extremely rapid, it is important to realize that, in practice, the minimum timescale can only be achieved for a certain range of small body sizes. Maximal accretion efficiency occurs for particles where the dominant effect of the gas is to dissipate particles' kinetic energy as they encounter the growing core, and all other aspects of the accretion process are set by the gravitational influence of the core. When the mass of small or large bodies is low enough that gas determines aspects of the accretion process, the timescale for growth will increase—in many cases, quite substantially.

To more thoroughly examine what causes the timescale to increase above ${t}_{\mathrm{Hill}}$ as the interaction between the gas and the small bodies increases in importance, we plot the four quantities that go into the calculation of ${t}_{\mathrm{grow}}$, along with timescale itself, in Figure 6. As α increases, the gas-particle interaction will dominate over the gravitational interaction between the small body and the core. For low values of α, we have ${R}_{\mathrm{WS}}\gt {R}_{H}\gt {H}_{p}$ and ${v}_{\infty }={v}_{\mathrm{shear}}$ (as expected for ${R}_{H}\,\gt {H}_{t}$), so the particle is in the 2D accretion regime and is able to accrete at ${t}_{\mathrm{Hill}}$. As α increases, the particle scale height begins to increase, until eventually ${H}_{p}\gt {R}_{H}$. At this point, accretion becomes less efficient and the timescale correspondingly increases. There are several other kinks in the timescale graph, which are caused by, in order of increasing α:

• 1.  
  R_WS decreasing to the point that ${R}_{\mathrm{WS}}\lt {R}_{H}$.
• 2.  
  ${v}_{\infty }$ being set by dispersion instead of shear, i.e., switching to the drift-dispersion dominated regime.
• 3.  
  Reaching H_p = H_g, at which point H_p stops increasing.

Zoom In Zoom Out Reset image size

For large α values, the growth timescale is almost an order of magnitude larger than ${t}_{\mathrm{Hill}}$, indicating a substantial slowdown in growth rate.

In this section, we compare our model to the analytic results from OK10 and Johansen & Lambrechts (2017, hereafter JL17)¹² in the laminar (α = 0 regime), as well as the extension of the OK10 model by Chambers (2014) to include turbulence. We discuss why our framework is more useful for incorporating the results of turbulence.

In order to facilitate comparison between the various models, we begin by highlighting some important features that emerge in our modeling of the growth rate. Figure 7 shows a plot of ${P}_{\mathrm{col}}\equiv 2{R}_{\mathrm{acc}}{v}_{\infty }$ for α = 0 in our model as a function of r_s and M. Several features, along with their analytic formulae, are marked on the plot. We discuss these features in detail below. In calculating the analytic expressions for these features, we implicitly assume a linear drag regime so that we can employ the analytic expressions for R_WS and R_shear (see Appendix A).

Zoom In Zoom Out Reset image size

In the upper left-hand corner of the plot, we see that accretion is cut off once ${R}_{b}\gt {R}_{\mathrm{stab}}$. At low core mass, ${R}_{\mathrm{stab}}={R}_{\mathrm{WS}}$ when this cutoff occurs, but R_shear dominates at high core masses, which causes the kink in the slope seen in the figure. The combination of these two scales causes a cutoff in accretion for

$$\begin{eqnarray}&&{St}\lt \min \left(3\displaystyle \frac{{v}_{H}^{3}{v}_{\mathrm{gas}}}{{c}_{s}^{4}},9\displaystyle \frac{{v}_{H}^{6}}{{c}_{s}^{6}}\right).\end{eqnarray} \tag{ 70 }$$

In the bottom of the plot, we see that accretion shuts off for particles that pass a certain maximum size at low core masses. In this regime, we have ${R}_{\mathrm{acc}}={R}_{\mathrm{WS}}$, ${v}_{\infty }\approx {v}_{\mathrm{gas}}$, and ${v}_{\mathrm{enc}}={v}_{\mathrm{kick}}$, which leads to a maximum size of

$$\begin{eqnarray}&&{St}=12{\left(\displaystyle \frac{{v}_{H}}{{v}_{\mathrm{gas}}}\right)}^{3}.\end{eqnarray} \tag{ 71 }$$

However, as particle size continues to increase, we see that accretion actually resumes once particles become large enough. This is caused by these large particles decoupling from the gas, which lowers ${v}_{\infty }={v}_{{pk}}$ and raises ${v}_{\mathrm{enc}}={v}_{{pg}}$ enough to overcome the increased mass of these particles. If we set ${R}_{\mathrm{acc}}={R}_{H}$, ${v}_{{pg}}\approx {v}_{\mathrm{gas}}$, and approximate v_pk by its Taylor Expansion at $\infty$ in the laminar regime, ${v}_{{pk}}\approx 2{v}_{\mathrm{gas}}/{St}$, then we can derive an approximate expression for the Stokes number at which accretion commences again:

$$\begin{eqnarray}&&{St}=\displaystyle \frac{{v}_{\mathrm{gas}}}{{v}_{H}}.\end{eqnarray} \tag{ 72 }$$

Note that, as particle size continues to increase, the energy criteria are only satisfied for a small range of particle sizes. However accretion continues probabilistically in this regime, in accordance with Equation (36).

The gap in particle size where accretion is not possible disappears once the core surpasses a certain mass. If particles that have their encounter velocity dominated by the particle-gas relative velocity, ${v}_{\mathrm{enc}}={v}_{{pg}}$, become available for accretion, then the range of particle sizes the core can accrete will be extended because these particles can dissipate more of their kinetic energy due to their larger encounter velocities. See Section 6.1 for more discussion of this effect. For these particles to be available, the previously derived upper limit on particle size, ${St}=12{v}_{H}^{3}/{v}_{\mathrm{gas}}^{3}$, must occur after the transition where ${v}_{{pg}}={v}_{\mathrm{kick}}$. Using the Taylor expansion of v_pg about 0 in the laminar regime, ${v}_{{pg}}\approx 2{v}_{\mathrm{gas}}{St}$, we see that the latter transition occurs at ${St}\approx {(3/4)}^{1/3}({v}_{H}/{v}_{\mathrm{gas}})$. Therefore, the mass scale where this transition occurs is given by

$$\begin{eqnarray}&&\displaystyle \frac{{v}_{H}}{{v}_{\mathrm{gas}}}={48}^{-1/3}.\end{eqnarray} \tag{ 73 }$$

Finally, at high core masses, the transition where KE = W occurs at a fixed Stokes number of

$$\begin{eqnarray}&&{St}=4\sqrt{3},\end{eqnarray} \tag{ 74 }$$

which is noted on the plot as well (see Section 5.2 for further discussion).

We now contrast the output from our model in the laminar regime to that of the analytic models of OK10 and JL17. Figure 8 below shows plots of both ${P}_{\mathrm{col}}\equiv 2{R}_{\mathrm{acc}}{v}_{\infty }$ and ${t}_{\mathrm{grow}}$ for a core at 30 au, using our model, OK10's, and JL17's¹³ (for a discussion of how the parameters of the protoplanetary disk are modeled, see Section 5.1). We plot P_col to disentangle the effects of particle scale height on the growth timescale, as we discuss below.

Zoom In Zoom Out Reset image size

In all three figures, we see the same broad features: while pebble accretion operates for a broad range of particle sizes, there exists an "optimal" range near St = 1 (${r}_{s}\approx 10\,\mathrm{cm}$) at which accretion reaches its maximal possible rate, and the growth timescale becomes comparable to ${t}_{\mathrm{Hill}}=M/(2{R}_{H}^{2}{\rm{\Omega }}{{\rm{\Sigma }}}_{p})$. In all three models, optimum accretion at ${St}\approx 1$ does not begin until the core has reached a "large enough" mass. Furthermore, as the core mass increases, the core's growth rate grows as well, in large part due to the increasing size of the core's Hill radius. The growth timescale, however, increases as the core grows, because the core's growth rate scales with M to a power less than one.

In our and OK10's modeling, we see two other notable features: first, there exists a gap in the range of particle sizes that can be accreted for low core masses. These particles cannot accrete in our framework because their incoming kinetic energies are too high. Second, on the right-hand side of the plots (${r}_{s}\gtrsim 100$ cm), particles have ${R}_{H}\lt {R}_{\mathrm{WS}},{R}_{\mathrm{shear}}$ (i.e., ${R}_{\mathrm{acc}}={R}_{H}$), but accretion is inefficient. This is due to the fact that these particles cannot fully dissipate their kinetic energy during a single orbital crossing; however, because of the chaotic nature of trajectories for particles with impact parameter $\sim {R}_{H}$, their growth timescale is increased by a factor of KE/W (see Section 4.1.1 for more detail.)

As can be verified from the plots, in the regions of parameter space where our model predicts that pebble accretion should operate, our results agree (within factors of a few) with both of these analytic models. The main disagreements are:

• 1.  
  The mass scale above which accretion occurs on timescales comparable to ${t}_{\mathrm{Hill}}$ is not exactly the same in the models. As discussed previously, our model has a transition in behavior past ${v}_{H}/{v}_{\mathrm{gas}}\approx {48}^{-1/3}$. In JL17's model, the parameters that set ${t}_{\mathrm{grow}}$ change for a mass $M\gt {M}_{t}$, where ${M}_{t}=\sqrt{1/3}{v}_{\mathrm{gas}}^{3}/(G{\rm{\Omega }})$. This is equivalent to ${v}_{H}/{v}_{\mathrm{gas}}=\sqrt{1/3}$, which is clearly comparable to the cut in our model. In OK10's modeling, the size of the gap progressively narrows, as the borders between the two regimes, ${St}=12{v}_{H}^{3}/{v}_{\mathrm{gas}}^{3}$ and ${St}={v}_{\mathrm{gas}}/{v}_{H}$, become comparable. Eventually, the gap disappears for ${v}_{H}/{v}_{\mathrm{gas}}=1$.
• 2.  
  As discussed above, in our modeling, accretion shuts off for particles below a certain radius, which stems from the fact that, once ${R}_{\mathrm{stab}}$ shrinks below R_b, we expect particles to flow around the core's atmosphere (see Figure 1, particularly the lower left-hand panel). The regions where this effect occurs are denoted by red hatching. In these regions, we expect growth by all mechanisms to be inhibited, e.g., capture by gravitational focusing or even physical collisions with the core will not occur. (In contrast, in the white regions of the plot, we still expect capture mechanisms other than gas-assisted growth to operate.) Because this process is set by the modification of the flow by the core's gravity, we would not expect OK10's integrations to capture it, because they assume a constant headwind velocity in their integrations. LJ12's simulations do not appear to reach high enough core masses with small enough particle sizes to capture this effect.¹⁴ This effect does, however, appear to be consistent with the results of Ormel (2013), who find minimal accretion for small particle sizes at high core masses (see their Figure 12).¹⁵
• 3.  
  JL17 define a "weak coupling" regime, where the stopping time of the particle exceeds the time to pass the protoplanet, ${t}_{p}={GM}/{({v}_{\mathrm{gas}}+{v}_{H})}^{3}$. For ${v}_{\mathrm{gas}}\gg {v}_{H}$, this criterion is equivalent to ${St}\gt 3{v}_{H}^{3}/{v}_{\mathrm{gas}}^{3}$, which is in line with our and OK10's criterion ${St}\gt 12{v}_{H}^{3}/{v}_{\mathrm{gas}}^{3}$. For ${v}_{H}\gg {v}_{\mathrm{gas}}$, the passing time criterion reduces to ${St}\gt 3$, which is again similar to our ${St}\gt 4\sqrt{3}$ limit. For particle sizes exceeding this critical Stokes number, JL17 exponentially smooth R_acc to zero, which differs from our and OK10's modeling.
• 4.  
  In the lower right-hand region, OK10's results and ours diverge. In this regime, particles that shear into R_H cannot fully dissipate their kinetic energy over one orbital crossing. This makes them accrete probabilistically, in accordance with Equation (36); i.e., instead of shutting off accretion because $W\lt {KE}$ in this regime, we increase the growth timescale by a factor KE/W. For particles with ${v}_{\mathrm{enc}}\sim {v}_{H}$, which holds for larger cores, this increases the growth timescale by a factor of St, which agrees with OK10's expression for ${R}_{\mathrm{acc}}$ in this regime: ${R}_{\mathrm{acc}}={R}_{H}/{St}$. However, for the low-mass cores in the bottom right-hand corner of the plot, the particle-gas relative velocity dominates over the core's gravity, changing the Stokes number dependence.

Thus, while our results are generally in order-of-magnitude agreement with other modeling in the laminar regime, there are few notable regions where our results diverge. It can be shown analytically that the agreement between our model and the other two essentially stems from the fact that, to order of magnitude, many of the transitions between length and velocity regimes occur in the same region of parameter space. For example, the transitions where ${R}_{H}={R}_{\mathrm{WS}}$ and ${v}_{{pk},{\ell }}={v}_{H}$ both occur for ${St}\sim {\zeta }_{w}$. Therefore, in the laminar case, there are essentially fewer regimes that need to be covered.

When turbulence is included, however, this is no longer the case. As an example, Figure 9 compares our modeling to the pebble accretion modeling of Chambers (2014). The expressions used in Chambers are analogous to the OK10 expressions, with the replacement of OK10's ${\zeta }_{w}\equiv \eta {v}_{k}/{v}_{H}$ parameter with ${v}_{\mathrm{gas}}/{v}_{H}$, where ${v}_{\mathrm{gas}}=\max (\eta {v}_{k},\sqrt{\alpha }{c}_{s})$. In our comparison, we neglect the exponential smoothing term they carry over from OK10, as well as terms involving eccentricity and inclination, as they are not included in our model.¹⁶

Zoom In Zoom Out Reset image size

As can be seen from the figure, for low core masses and small body sizes, our models are in order-of-magnitude agreement. Again, analytic calculations can be used to demonstrate that, for low core masses and small particle sizes, the expressions between the two models agree to order of magnitude. For example, the cutoff in Stokes number employed by Chambers—${{St}}_{\mathrm{crit}}=12{({v}_{H}/{v}_{\mathrm{gas}})}^{3}$ is replicated in our model for small cores and small particle radii. However, as can also be seen in the figure, there are numerous regimes covered in our modeling that are not captured by extending the OK10 expressions. A prominent example is the feature in the lower right-hand corner, where low-mass cores can accrete "large" particles ($\sim {10}^{2}\mbox{--}{10}^{3}$ cm) on rapid timescales. This feature results from the fact that, for low-core masses, the velocities of these larger particles are set by their interactions with the nebular gas. Because these large particles are not well-coupled to the gas flow, they have low kinetic energies relative to the core, but high velocities relative to the gas, allowing a certain range of particle radii to dissipate an amount of energy during the encounter that is comparable to their kinetic energy. The shape of this feature depends on the transitions where ${KE}\sim W$ as well as where ${R}_{\mathrm{WS}}\sim {R}_{H}$ and ${v}_{{pk}}\sim {v}_{H}$. In the presence of turbulence, these transitions no longer need occur at similar regions of parameter space, leading to a more complex model of pebble accretion that is not well-captured by extending the OK10 expressions.

We close this section by comparing the results of our model to numerical simulations of pebble accretion in the presence of turbulence. Xu et al. (2017) performed magnetohydrodynamic (MHD) simulations of gas-assisted growth of turbulence due to the magnetorotational instability (MRI). These simulations provide an excellent point of comparison for our order-of-magnitude model, because they use much more detailed physics but apply over a narrower range of parameter space. As we will show below, our model agrees with the Xu et al. results to order of magnitude, as we would expect.

Xu et al. perform six sets of simulations; they use two different values of core mass, at three different levels of turbulence for each core. The core masses are parameterized in terms of the "thermal mass," which Xu et al. give as

$$\begin{eqnarray}&&{M}_{T}\approx 160{M}_{\oplus }{\left(\displaystyle \frac{a}{30\mathrm{au}}\right)}^{3/4}.\end{eqnarray} \tag{ 75 }$$

Note that this is the same as the flow isolation mass (see Section 7) for the values used in their work. The simulations are performed at core masses of ${\mu }_{c}\equiv M/{M}_{T}=3\times {10}^{-3}$ and ${\mu }_{c}=3\times {10}^{-2}$. The three levels of turbulence are achieved by performing a pure hydrodynamic simulation, a non-ideal MHD simulation with ambipolar diffusion, and an ideal MHD simulation. Thus far in our work, we have assumed that $\alpha \approx {\alpha }_{z}$, where ${\alpha }_{z}$ is the diffusion coefficient for the particle scale height, but Xu et al. are able to separately calculate α and ${\alpha }_{z}$. We also use their value for v_z to calculate our α value, as we use α to parameterize the turbulent gas velocity. For the pure hydrodynamic run, they have ${v}_{z}={\alpha }_{z}=0$. For the ambipolar diffusion run, they give $\langle {v}_{z}^{2}/{c}_{s}^{2}\rangle =3.0\times {10}^{-4}$ and ${\alpha }_{z}=7.8\times {10}^{-4}$, whereas they give $\langle {v}_{z}^{2}/{c}_{s}^{2}\rangle =1.21\times {10}^{-2}$ and ${\alpha }_{z}=4.4\times {10}^{-3}$ for the ideal MHD simulation.

A comparison between our analytic model and their simulations is shown in Figure 10. This figure plots ${k}_{\mathrm{abs}}\equiv \dot{M}/{\dot{M}}_{\mathrm{Hill}}$ as a function of the particle Stokes number St. Here, $\dot{M}=M/{t}_{\mathrm{grow}}$ is the growth rate of the core, and ${\dot{M}}_{\mathrm{Hill}}=M/{t}_{\mathrm{Hill}}$ is the growth rate for Hill accretion. The solid lines depict the model presented in this paper, while the data points show the results from the Xu et al. paper. For the purposes of matching their model, we have used a temperature profile consistent with their results, as well as their values of ${\alpha }_{z}$ for the calculation of H_p (as opposed to using α), and set the scale height in the laminar case to be ${H}_{p,\mathrm{lam}}=0.01{H}_{g}$ to be consistent with their choice.

Zoom In Zoom Out Reset image size

As can be seen in Figure 10, our model nicely achieves the intended goal of reproducing the trends found by the numeric results, as well as being accurate to within a factor of two for all of the values found in the numeric simulations. It is worth noting here that the simulations are carried out at quite large values of core mass, and serve mostly to confirm our prediction that the efficiency of pebble accretion is not hugely reduced by increasing turbulence. It would be interesting to investigate in future simulations whether the drop-off in efficiency for lower core masses predicted in our model (see Section 6.1) is reproduced in the numeric simulations.

In Section 5.3, we identify the minimal timescale ${t}_{\mathrm{Hill}}$ on which pebble accretion can operate. We also emphasize, however, that pebble accretion can operate on timescales that are substantially slower—in many cases, these timescales can exceed the lifetime of the disk, ${\tau }_{\mathrm{disk}}\sim 2.5\,\mathrm{Myr}$, which implies that pebble accretion essentially will not occur. In this section, we highlight how, at low core masses and wide orbital separations, only a small range of particle sizes, if any, accrete on timescales comparable to ${t}_{\mathrm{Hill}}$.

To illustrate how core mass and orbital separation affect the growth timescale, we plot ${t}_{\mathrm{grow}}$ as a function of both r_s and M for a core at a = 5 au in Figure 11, and for a core at a = 50 au in Figure 12. To begin, we note that growth is generally slower at wide orbital separation, because the dynamical time ${t}_{\mathrm{dyn}}\sim {{\rm{\Omega }}}^{-1}$ is larger and the solid surface density is smaller. For low core mass, there is also a gap in the particle sizes that can be accreted. This gap is small at low turbulence, but as α increases, the width of this gap increases as well. This effect is enhanced at wide orbital separation, where the gap extends to higher core masses and is overall wider. Thus, low-mass cores can often only accrete at small pebble radii. These small particles are accreted on slow timescales that are $\gtrsim {\tau }_{\mathrm{disk}}$, as small pebbles can only be captured at low impact parameters and have low densities because they can be easily lofted by turbulence. Thus, we see that growth is much slower at low core mass and wide orbital separation, particularly as the strength of turbulence is increased. If there exists a population of ≳10 m "boulders," then we see from the figures that there exists a narrow range of these larger particle sizes for which accretion can proceed on rapid timescales, even if accretion of pebble-sized particles is slow.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Once the core exceeds a certain critical mass, the gap disappears and accretion proceeds much more rapidly, allowing pebbles to accrete on timescales comparable to ${t}_{\mathrm{Hill}}$. To emphasize this effect, in Figure 13 we plot the range of particle sizes that can accrete on timescales within a factor of two of ${t}_{\mathrm{Hill}}$. Each panel depicts a different core mass, while different hatching patterns depict different levels of turbulence. As expected from the above discussion, for low core mass, only close-in cores with low levels of nebular turbulence can accrete particles on timescales comparable to ${t}_{\mathrm{Hill}}$. As core mass increases, however, particles farther out in the disk can be accreted with growth timescales $\sim {t}_{\mathrm{Hill}}$, even when turbulence is strong.

Zoom In Zoom Out Reset image size

We shall now discuss more quantitatively what sets the range of particle sizes that can be accreted, as well as how that range is affected by changing the core mass and orbital separation. For low-mass cores accreting small particles, we expect R_WS to set the impact parameter. Because the particles are so small, and therefore well-coupled to the gas, we will have ${v}_{{pk}}\approx {v}_{\mathrm{gas}}$ and ${v}_{{pg}}\ll {v}_{\mathrm{kick}}$. Thus, in this regime, we expect ${v}_{\infty }\approx {v}_{\mathrm{gas}}$ (because shear is negligible for low-mass cores) and ${v}_{\mathrm{enc}}={v}_{\mathrm{kick}}$ (because ${v}_{\mathrm{gas}}\lt {v}_{\mathrm{orbit}}$ in this regime). Plugging these values into our expressions for KE and W, we can show that having ${KE}\lt W$ requires

$$\begin{eqnarray}&&{St}\lt 12\displaystyle \frac{{v}_{H}^{3}}{{v}_{\mathrm{gas}}^{3}},\end{eqnarray} \tag{ 76 }$$

which was given previously in Equation (71). This is one facet that makes accretion less efficient for low-mass cores, high turbulence, and wide orbital separation—because ${v}_{H}/{v}_{\mathrm{gas}}\propto {M}^{1/3}{a}^{-4/7}$ in the laminar regime (i.e., approximating ${v}_{\mathrm{gas}}\approx \eta {v}_{k}$) and $\propto \,{M}^{1/3}{a}^{-2/7}{\alpha }^{-1/2}$ in the turbulent regime (${v}_{\mathrm{gas}}\approx \sqrt{\alpha }{c}_{s}$), decreasing core mass, increasing α, or increasing orbital distance will all make the limit on St more stringent.

As small body radius increases, ${v}_{\mathrm{kick}}$ decreases due to the increasing size of the core's WISH radius, while v_pg rises as particles decouple from the gas. Eventually, we reach the point where ${v}_{{pg}}\gt {v}_{\mathrm{kick}}$, meaning that now ${v}_{\mathrm{enc}}={v}_{{pg}}$. In this regime, the dependence of the energy on small body radius becomes more complex: throughout a large amount of parameter space, KE/W increases with particle size, which encapsulates the fact that it is more difficult for heavier particles to dissipate their kinetic energy. If ${v}_{\mathrm{enc}}={v}_{{pg}}$, however, the work done during the encounter increases strongly with small body radius, which actually makes it easier for large particles to be accreted. Furthermore, if ${v}_{\infty }={v}_{{pk}}$, then the incoming KE of particles can also decrease with small body radius, making it even easier for heavier particles to be accreted.

This effect is illustrated in Figure 14, which plots the ratio KE/W as a function of r_s for different core masses. The two different panels show cores at different orbital separations. For the low-mass cores in the right-hand panel, we see that once r_s increases past a certain value, the qualitative behavior of KE/W changes—instead of monotonically increasing, the slope flattens out or even decreases. However, for the very low-mass cores, this behavior is inconsequential because particles are ruled out from accreting before we reach a large enough particle size that ${v}_{\mathrm{enc}}={v}_{{pg}}$ and this more complex behavior starts. For larger core masses, however, this is no longer the case, and the range of particle sizes that can be accreted is greatly extended. This is what causes the gap seen in the range of accreted particle sizes to disappear—once particles with ${v}_{\mathrm{enc}}={v}_{{pg}}$ become available, accretion generally continues until the limit ${St}=4\sqrt{3}$ discussed in Equation (61) is reached. From comparison of the two panels, we also see that this effect is much more prominent at wide orbital separations—at 1 au, none of the cores exhibit the change in slope seen at 50 au.

Zoom In Zoom Out Reset image size

Thus, we see that growth changes qualitatively once particles with ${v}_{{pg}}\gt {v}_{\mathrm{kick}}$ can be accreted. The critical value of mass where this occurs can be approximated by calculating the mass at which ${v}_{{pg}}\gt {v}_{\mathrm{kick}}$ is reached before the critical Stokes number given in Equation (76) is reached. If we keep the approximations that led to Equation (76) but take ${v}_{{pg}}\approx 2{v}_{\mathrm{gas}}{St}$ in the laminar regime and ${v}_{{pg}}\approx {v}_{\mathrm{gas}}\sqrt{{St}}$ in the turbulent regime, we can solve for the Stokes number past which ${v}_{{pg}}\gt {v}_{\mathrm{kick}}$. Doing so and setting the resulting Stokes number equal to (76) yields a mass limit that can, in both regimes, be approximated by

$$\begin{eqnarray}&&\displaystyle \frac{{v}_{H}}{{v}_{\mathrm{gas}}}\approx {48}^{-1/3}.\end{eqnarray} \tag{ 77 }$$

This inefficiency of pebble accretion at lower core masses is also identified by Visser & Ormel (2016), who calculate the growth timescale of planetesimals by numerically integrating the pebble equation of motion in the presence of two different laminar gas flow patterns. Visser & Ormel find that pebble accretion is only faster than gravitational focusing (which they refer to as the "Safronov regime") once the growing planetesimal exceeds a certain radius, R_PA. They also find that this critical radius increases further out in the disk, which is also consistent with our results. It can be seen from our analytic expressions, however, that while these effects are present in our model even in the laminar case, they are strongly amplified by the presence of turbulence.

In this section, we consider the effects of varying the prefactor to the temperature profile, T₀. While changing the disk temperature will affect the properties of the small bodies, such as number density and composition, by changing what volatile species are present in solids at a given location, we neglect such effects in what follows. The consequences of changing the temperature profile are complex and depend on the local disk parameters as well as the core mass and strength of turbulence. Nevertheless, in general, increasing the temperature is detrimental to the accretion rate when gas-dominated processes set the scales relevant to growth. For growth at ${t}_{\mathrm{Hill}}$, however, none of the scales depend on temperature, so the most rapid growth timescales are unaffected by increasing the disk temperature.

The primary effect of varying T₀ on growth timescales is due to the dependence of c_s on T: ${c}_{s}\propto {T}^{1/2}$. Because ${H}_{\mathrm{KH}}\propto {c}_{s}^{2}$, and ${H}_{t}\propto {c}_{s}$, increasing T₀ will increase H_p. Similarly, ${v}_{{pk},{\ell }}\propto {c}_{s}^{2}$ and ${v}_{{pk},t}\propto {c}_{s}$, so increasing T₀ will also increase ${v}_{\infty }$ in the drift-dominated regime. Finally, from inspection of Equation (40), we see that T₀ affects R_WS by changing ${\rho }_{g}\,(\propto {c}_{s}^{-1}$), v_th, and ${v}_{\mathrm{gas}}$. By inspecting the scaling of R_WS one can show that R_WS is a decreasing function of T₀.

While some of the above effects increase ${t}_{\mathrm{grow}}$, while other decrease it, the general effect of increasing T₀ is to slow down growth. To see this, we can make the same approximations described in Equations (64)–(66). As in Section 5.3, we consider the 2D and 3D regimes separately.

In the 2D regime where ${H}_{p}\lt {R}_{\mathrm{acc}}$, ${H}_{\mathrm{acc}}={R}_{\mathrm{acc}}$, so we have ${t}_{\mathrm{grow}}\propto {({R}_{\mathrm{acc}}{v}_{\infty })}^{-1}$. Having ${R}_{\mathrm{acc}}\gt {H}_{p}$ implies that ${v}_{\mathrm{shear}}={R}_{\mathrm{acc}}{\rm{\Omega }}\gt {v}_{{pk}}$, which in turn implies that ${R}_{\mathrm{shear}}\lt {R}_{\mathrm{WS}}$. Thus, we expect that ${t}_{\mathrm{grow}}$ is independent of T₀ in this regime.

For 3D accretion, we have ${H}_{p}\gt {R}_{\mathrm{acc}}$, which implies that ${v}_{\infty }={v}_{{pk}}$. Therefore we have ${t}_{\mathrm{grow}}\propto {({R}_{\mathrm{WS}}^{2}{\rm{\Omega }})}^{-1}$, which is an increasing function of T₀. We therefore expect that (approximately) higher T₀ should lead to slower growth. An example of the difference caused by increasing T₀ is shown in Figure 15.

Zoom In Zoom Out Reset image size

A secondary effect of increasing T₀ is that t_s and St, which depend on ρ_g, v_th, and λ, are also affected. However, inspection of Equation (13) shows that the effects of T₀ on t_s cancel out in the Epstein regime. Therefore, T₀ only affects St in the fluid regime. This substantially diminishes the importance of this dependence, because the small body sizes that grow the most efficiently are in the Epstein regime for a large amount of parameter space. For particles that are in the fluid regime, higher values of T₀ will shift values of St to lower values of r_s.

Because of the approximate cancellation between ${v}_{\infty }$ and H_p in the 3D regime, the maximal increase in timescale that can be provided by increasing T₀ is generally of order ${t}_{\mathrm{hot}}/{t}_{\mathrm{cold}}\sim {R}_{\mathrm{acc},\mathrm{cold}}^{2}/{R}_{\mathrm{acc},\mathrm{hot}}^{2}$, where hot and cold denote the higher and lower values of T₀ respectively. Because R_H generally represents an upper limit on R_acc (ignoring accretion in the ${R}_{\mathrm{acc}}={R}_{b}$ regime), the maximal increase is of order ${({R}_{{WS},\mathrm{cold}}/{R}_{{WS},\mathrm{hot}})}^{2}$, which can be of order ${({c}_{s,\mathrm{hot}}/{c}_{s,\mathrm{cold}})}^{3}\,={({T}_{0,\mathrm{hot}}/{T}_{0,\mathrm{cold}})}^{3/2}$ in the laminar Stokes and Ram regimes. In practice, this maximal value is rarely reached because it requires satisfying ${r}_{s}\gt 9\lambda /4$ to be in the fluid regime, while simultaneously having ${St}\lt {v}_{\mathrm{gas}}^{3}/(3{v}_{H}^{3})$ in order to have ${R}_{\mathrm{WS}}\lt {R}_{\mathrm{shear}}$. In practice, therefore, the increase in growth timescale generally goes as ${T}_{0,\mathrm{hot}}/{T}_{0,\mathrm{cold}}$. Because accretion timescales at t_Hill are not dependent on temperature, increasing T₀ does not affect the rapid growth rates supplied by gas-assisted growth.

Changing T₀ will also have an effect on the range of particle sizes that the growing core is able to accrete. For small particle sizes, the accretion cutoff is determined by the point where ${R}_{b}={R}_{\mathrm{WS}}$. Because ${R}_{b}\propto {c}_{s}^{-2}$, which is stronger than the dependence of R_WS on c_s regardless of drag regime, increasing T₀ will decreases the size of the Bondi radius relative to the WISH radius. This, in turn, will allow the core to accrete smaller-sized bodies. For larger particle sizes, the effect of increasing T₀ is not as clear-cut. For particles accreting at t_Hill, the only effect of increasing T₀ is to increase the drag force on accreting small bodies, which allow the core to accrete larger sizes. In other regimes, however, increasing T₀ can also increase the approach velocity of small bodies. In these cases, the increase in kinetic energy of accreting particles outweighs the larger amount of work done, and the maximal size of small bodies accreted is reduced.

In this section, we explore growth in a disk where the gas density has been reduced by a factor of 100, but the solid surface density is unchanged; i.e., we have ${\rm{\Sigma }}={{\rm{\Sigma }}}_{g,0}{(a/\mathrm{au})}^{-1}$, where ${{\rm{\Sigma }}}_{g,0}={{\rm{\Sigma }}}_{0}/100$, and ${{\rm{\Sigma }}}_{0}=500\,{\rm{g}}\,{\mathrm{cm}}^{-2}$ is the prefactor employed elsewhere in this paper for the gas surface density. However, we keep ${{\rm{\Sigma }}}_{p,0}=5\,{\rm{g}}\,{\mathrm{cm}}^{-2}$. We note that this may affect some of the expressions used in this work, which implicitly assume ${\rho }_{g}\gg {\rho }_{p}$, where ${\rho }_{p}$ is the volumetric mass density of the small bodies. We neglect any such effects in what follows. These choices for gas and solid densities are made to emulate the conditions in the disk when the gas component of the disk is in the process of photoevaporating, which is an important stage in some theories of planet formation (e.g., Lee & Chiang 2016; Frelikh & Murray-Clay 2017). As we shall show below, the predominant effect of reducing the gas surface density is to shift the range of small body sizes where accretion occurs to lower values.

An example of growth in a depleted disk is shown in Figure 16. As can be seen in the figure, the predominant effect of changing Σ is to shift growth down to lower values of r_s. This is due to the fact that the quantities that go into calculating t_grow, even the energy criteria, are functions of r_s through their dependence on St alone. In the Epstein regime, ${t}_{s}\propto {r}_{s}/{\rho }_{g}$, and in the Stokes regime, ${t}_{s}\propto {r}_{s}^{2}/{\rho }_{g}$, so the radius corresponding to a given Stokes number decreases when the surface density, and correspondingly the volumetric density, are decreased. This is what causes the shift to lower radii seen in Figure 16. Other than this shift, however, there is essentially no change in the growth timescale. In other words, when the timescale is viewed as a function of Stokes number, i.e., if we consider ${t}_{\mathrm{grow}}({St})$ as opposed to ${t}_{\mathrm{grow}}({r}_{s})$, then this function is independent of Σ.

Zoom In Zoom Out Reset image size

The sole caveat is that, for particles not in a linear drag regime, the Stokes number of a particle is now dependent on the particle's velocity, which means that the Stokes numbers used for calculating different quantities might not be the same. For example, we can express the WISH radius as ${R}_{\mathrm{WS}}=\sqrt{{GMSt}/({v}_{\mathrm{gas}}{\rm{\Omega }})}$, but the St value in that expression is defined with respect to ${v}_{\mathrm{rel}}={v}_{\mathrm{gas}}$, meaning it will not be the same as the Stokes number for, e.g., v_pg. In this case, the simple argument that all quantities are solely dependent on St breaks down. In practice, if we use, e.g., the Stokes number defined with respect to laminar drift velocities for comparison purposes, the discrepancy between the two surface densities is minor. Furthermore, if we write the critical radius dividing the fluid and diffuse drag regimes, r_s = 9λ/4, in terms of Stokes number, it is straightforward to show that the particle is in the fluid regime for

$$\begin{eqnarray}&&{{St}}_{\mathrm{crit}}\gt 7.4\times {10}^{-2}{\left(\displaystyle \frac{a}{\mathrm{au}}\right)}^{23/7}{\left(\displaystyle \frac{{{\rm{\Sigma }}}_{g,0}}{500{\rm{g}}{\mathrm{cm}}^{-2}}\right)}^{-2}.\end{eqnarray} \tag{ 78 }$$

Because of the strong scaling of St_crit with a, the differences between the two surface densities disappear entirely for a ≳ 5 au.

Thus, in general, the range of small body sizes where accretion is effective will shift to lower values as the gaseous component of the disk dissipates. The subsequent effects of such a shift are quite sensitive to the underlying size distribution of the small bodies. A particularly salient issue here is whether radius or Stokes number controls the processes that produce the size distribution, as could be the case if, e.g., fragmentation near St ∼ 0.1 generates an upper cutoff to growth (Blum & Wurm 2008). If the Stokes number is what matters, then the effects of dissipating the surface density would be rather minimal, provided the surface density of the disk evolves faster than the disk dissipation timescale. If the size distribution is determined by radius, however, and there exists a small range of sizes where most of the mass of the disk is located, then the shift in effective accretion range could have substantial effects, either positive or negative, on the accretion rate of small bodies. In the outer disk, for example, the combination of growth and radial drift may set an important particle scale (e.g., Birnstiel et al. 2012; Powell et al. 2017). The combination of these effects merits further investigation. See Lambrechts & Johansen (2014) for an example of pebble accretion in the presence of particle sizes determined by the interplay of growth and drift (note that these authors use a full gas disk, not one depleted in gas surface density).

In this section, we discuss the dynamical effect of varying M_*. Clearly, higher-mass stars will have a higher overall temperature, and should on average have higher surface densities as well. Because we have discussed those effects in previous sections, in this section we consider only the effect of varying M_*, leaving the other properties of the disk unchanged.

Changing M_* has an effect on the growth process solely though its effect on the local Keplerian orbital frequency Ω and the size of the Hill radius R_H. Higher-mass stars have lower dynamical times, which tends to speed up growth. On the other hand, in the presence of more massive stars, the Hill radius of a growing planet will shrink as it becomes more difficult to hold on to accreting material, which clearly inhibits growth. The interplay between these two factors will determine the net effect of varying the stellar mass—which effect dominates, and therefore whether changing stellar mass is beneficial or detrimental to growth, depends on the other input parameters.

An example is shown in Figure 17. The figure shows a plot of t_grow versus r_s for five different stellar masses; the panels shows the growth timescale for two different levels of turbulence. Many of the main features of varying M_* are visible in the figure. For the α = 0 case, small particles (r_s ≲ 0.1 mm) are accreted much more rapidly by the more massive stars. In this regime, particles have R_acc = R_WS and accrete in 3D, with scale height ${H}_{p}={H}_{\mathrm{KH}}.$ Thus, the primary effect here of varying stellar mass is to change $\eta {v}_{k}$ and ${\rho }_{g}$ along with Ω. In the Epstein regime, however, these effects on both R_WS and St cancel out, and because ${H}_{\mathrm{KH}}\sim 2\eta {v}_{k}/{\rm{\Omega }}$, the increase in growth rate in this regime scales roughly as ${\rm{\Omega }}\propto {M}_{* }^{1/2}$. Eventually, the particle size increases to the point that ${R}_{\mathrm{shear}}\lt {R}_{\mathrm{WS}}$, which also implies that particles are now shear-dominated. Because this change roughly corresponds to where we shift from 3D to 2D accretion in the laminar regime, the effect of increasing stellar mass in this regime is actually to increase the growth timescale, ${t}_{\mathrm{grow}}\propto {({R}_{\mathrm{shear}}^{2}{\rm{\Omega }})}^{-1}\propto {{St}}^{-2/3}{R}_{H}^{-2}{{\rm{\Omega }}}^{-1}\propto {M}_{* }^{1/6}$, because St is independent of M_* in the Epstein regime. This is what causes the overlap in ${t}_{\mathrm{grow}}$ curves as r_s increases. Due to the fact that changing M_* can switch the regime that determines one of our growth parameters (e.g., shear- versus dispersion-dominated) the maximal change in growth rate for two stellar masses ${M}_{* ,1}$ and ${M}_{* ,2}$ in this regime can be complicated. In general, however, the change is of order ${{\rm{\Omega }}}_{1}/{{\rm{\Omega }}}_{2}={({M}_{* ,1}/{M}_{* ,2})}^{1/2}$, as can be seen in Figure 17.

Zoom In Zoom Out Reset image size

As particle size increases, the cores reach the point where they can accrete particles at ${R}_{\mathrm{acc}}={R}_{H}$. As can be seen in Figure 17, the value of r_s where this transition occurs is a decreasing function of M_* due to the decreased size of the Hill radius at higher stellar mass. Figure 17 also shows that increasing stellar mass initially decreases the growth timescale, but further increasing the mass of the star actually increases the growth timescale. This is due to the fact that decreasing M_* increases H_KH, affecting whether the largest particles are accreting in the 2D or 3D regime. Low-mass stars will accrete in 3D, in which case increasing M_* will decrease the growth timescale by bringing the particle density down and the rate of shear up. Once the stellar mass becomes sufficiently large, however, we will have ${R}_{H}\gt {H}_{\mathrm{KH}}$, and the core will accrete in 2D. In this case, increasing M_* will slow down growth by decreasing the size of the Hill radius.

As turbulence increases, the difference between the various stellar masses for small particle size disappears. This effect occurs due to a balance between the decreased growth rate from decreasing encounter rate and an increased growth rate due to an increased particle density, because we generally accrete in 3D in this regime. If particles are drift-dispersion dominated, then ${v}_{\infty }={v}_{{pk}}\approx {v}_{\mathrm{gas}}$ and ${R}_{\mathrm{acc}}={R}_{\mathrm{WS}}$. When $\alpha \ne 0$, R_WS is no longer independent of M_*. If $\alpha \gt \eta$, then ${v}_{\mathrm{gas}}$ is approximately constant with respect to varying Ω, meaning that the quantity in the denominator of ${t}_{\mathrm{grow}}$—${R}_{\mathrm{WS}}^{2}{v}_{{pk}}\propto {M}_{* }^{-1/2}$. Furthermore, if ${v}_{\mathrm{shear}}\gt {v}_{{pk}}$, then ${R}_{\mathrm{acc}}={R}_{\mathrm{shear}}$, and the denominator of ${t}_{\mathrm{grow}}$ is $\propto {R}_{\mathrm{shear}}^{3}{\rm{\Omega }}\propto {M}_{* }^{-1/2}$. For sufficiently strong turbulence that ${H}_{\mathrm{turb}}\gt {H}_{\mathrm{KH}}$, we have ${H}_{p}\propto {{\rm{\Omega }}}^{-1}\propto {M}_{* }^{-1/2}$. Thus, these effects approximately cancel out for 3D accretion in the strong turbulence regime, causing the growth curves to merge.

For high core masses, the transition on the right side of the graphs where ${t}_{\mathrm{grow}}$ increases as a function of r_s occurs at essentially the same value of small body radius, independent of M_*. As noted in Equation (61), the energy criterion ${KE}\lt W$ can be rewritten as ${St}\lt 4\sqrt{3}$ for particles in this growth regime. For particles in the Epstein regime the Stokes number is independent of M_*, causing the cutoff radius to be approximately the same for all stellar masses. The timescale in this regime is weakly dependent on M_*—these large particles accrete at ${t}_{\mathrm{Hill}}$, with an increase in growth timescale $\propto {St}$, which again is essentially independent of M_*. Thus, the growth timescale goes as ${t}_{\mathrm{grow}}\propto {({R}_{H}{v}_{H})}^{-1}\propto {M}_{* }^{1/6}$.

In the laminar regime, the cutoff at small radii is also independent of r_s. This cutoff occurs when ${R}_{\mathrm{WS}}\lt {R}_{b}$, and because, as previously noted, R_WS is independent of M_* in the laminar regime, this cutoff is independent of M_* as well. As turbulence increases, R_WS becomes a decreasing function of M_*, causing the cutoff in growth to shift to higher values of r_s.

For low core masses, the low-end cutoff still has the same dependence on M_*. However, growth will now cut off at the limit given by Equation (76), which is an increasing function of M_* in both the laminar and turbulent regimes.