Highlight
Free Access
Issue
A&A
Volume 520, September-October 2010
Article Number A43
Number of page(s) 15
Section Planets and planetary systems
DOI https://doi.org/10.1051/0004-6361/201014903
Published online 29 September 2010
A&A 520, A43 (2010)

The effect of gas drag on the growth of protoplanets

Analytical expressions for the accretion of small bodies in laminar disks

C. W. Ormel - H. H. Klahr

Max-Planck-Institut für Astronomie, Königstuhl 17, 69117 Heidelberg, Germany

Received 30 April 2010 / Accepted 1 July 2010

Abstract
Planetary bodies form by accretion of smaller bodies. It has been suggested that a very efficient way to grow protoplanets is by accreting particles of size $\ll$km (e.g., chondrules, boulders, or fragments of larger bodies) as they can be kept dynamically cold. We investigate the effects of gas drag on the impact radii and the accretion rates of these particles. As simplifying assumptions we restrict our analysis to 2D settings, a gas drag law linear in velocity, and a laminar disk characterized by a smooth (global) pressure gradient that causes particles to drift in radially. These approximations, however, enable us to cover an arbitrary large parameter space. The framework of the circularly restricted three body problem is used to numerically integrate particle trajectories and to derive their impact parameters. Three accretion modes can be distinguished: hyperbolic encounters, where the 2-body gravitational focusing enhances the impact parameter; three-body encounters, where gas drag enhances the capture probability; and settling encounters, where particles settle towards the protoplanet. An analysis of the observed behavior is presented; and we provide a recipe to analytically calculate the impact radius, which confirms the numerical findings. We apply our results to the sweepup of fragments by a protoplanet at a distance of 5 AU. Accretion of debris on small protoplanets ($\la$50 km) is found to be slow, because the fragments are distributed over a rather thick layer. However, the newly found settling mechanism, which is characterized by much larger impact radii, becomes relevant for protoplanets of $\sim$103 km in size and provides a much faster channel for growth.

Key words: planets and satellites: formation - protoplanetary disks - minor planets, asteroids: general

1 Introduction

We consider how gas drag affects the collision rates between a big body - a planetesimal or protoplanet - and small particles, e.g., dust, chondrules, or boulders. Although the core accretion model (Hubickyj et al. 2005; Pollack et al. 1996) in its initial stages, i.e., before the formation of a $\sim$10 Earth mass ($M_\oplus$) core, concerns the accumulation of solid bodies, the role of the gas cannot be overstated. In the early phases of planet formation - the growth of dust to planetesimals - the gas damps the velocities of small particles. Initially, the (relative) velocities between particles are tiny and this is the reason why dust grains can coagulate due to intermolecular forces (Blum & Wurm 2000; Dominik & Tielens 1997) - an effect much harder to envision in the diffuse interstellar medium or even in molecular clouds. In this stage mechanisms that induce a relative velocity among the dust particles include Brownian motion, settling, radial drift, and turbulent motions. The latter three are all functions of the particle's stopping time, a measure of how well particles couple to the gas. With increasing size (or, more correctly, increasing mass-to-surface area) particles couple less well to the gas and relative velocities increase, culminating in the so called meter-size barrier, which, at our current level of understanding, can best be overcome by the combined efforts of turbulent concentration and gravitational collapse (Cuzzi et al. 2010; Johansen et al. 2007,2009).

Gas drag also affects the collisional behavior at a much later stage, when protoplanets accrete planetesimals of perhaps $\sim$1-102 km in size. The collisional cross section between these big bodies is increased by gravitational focusing, i.e., the body can accrete particles at a cross section larger than its geometrical cross section due to gravitational deflection (Wetherill & Stewart 1989; Safronov 1969; Greenzweig & Lissauer 1992,1990). This effect, however, is very sensitive to the velocity $v_{\rm a}$ at which the bodies approach: if $v_{\rm a}$ is too large, the focusing vanishes. In planetesimal accretion theory it is expected that a protoplanet will excite the random motions (eccentricities and inclinations) of the bodies it is accreting from, leading to a self-regulated accretion behavior, which slows down the growth (Ida & Makino 1993; Kokubo & Ida 1998; Ormel et al. 2010a). Gas drag can provide some relief since, by damping the random motions of the planetesimals, the gravitational focusing is kept large. Moreover, the capture probability of planetesimals is also significantly increased when (proto)planets are surrounded by atmospheres (Tanigawa & Ohtsuki 2010; Inaba & Ikoma 2003) - again, gas drag is the mechanism that facilitates their accretion. Still, it is unclear if these effects are sufficient to overcome the timescale problem, i.e., to grow protoplanets to $\sim$10 $M_\oplus$ within the time the gas disk dissipates ($\sim$106 yr); see Levison et al. (2010) for a recent review.

Due to the dynamical heating of planetesimals, planetesimal-planetesimal collisions may become disruptive, producing smaller planetesimals or even fragments (Leinhardt et al. 2009; Wetherill & Stewart 1993). These fragments can be kept dynamically cold, e.g., by mutual collisions or by gas drag. The accretion then takes place at low $v_{\rm a}$ - the shear-dominated regime - which is very favorable for growth (Goldreich et al. 2004). The generation of large amounts of fragments therefore can significantly boost accretion. In particular, the accretion rate in the two dimensional (interactions are confined to a plane), gas-free, three-body regime (including the gravity of the central star) is derived by a number of studies to be

\begin{displaymath}\left(\frac{{\rm d}M}{{\rm d}t}\right)_{\rm gf} \approx 11 \alpha_{\rm p}^{1/2} R_{\rm h} v_{\rm h} \Sigma
\end{displaymath} (1)

(e.g., Rafikov 2004; Ida & Nakazawa 1989; Greenberg et al. 1991; Ormel et al. 2010b; Weidenschilling 2005; the numerical constant is adopted from Inaba et al. 2001), where $\Sigma$ is the density in solids, $R_{\rm h}$ the Hill radius,

\begin{displaymath}R_{\rm h} = a_0 \left( \frac{M_{\rm p}}{3M_\star} \right)^{1/3},
\end{displaymath} (2)

$v_{\rm h}$ the Hill velocity, $v_{\rm h} = R_{\rm h} \Omega$, a0 the semi-major axis, $\Omega$ the corresponding orbital frequency, $M_{\rm p}/M_\star$ the ratio between the mass of the protoplanet and the central star, and $\alpha_{\rm p}$ the ratio between the protoplanet radius and the Hill radius, $\alpha_{\rm p} = R_{\rm p}/R_{\rm h}$. Equation (1) is often used in statistical models for the accretion rate (e.g., Kobayashi et al. 2010; Chambers 2008,2006; Brunini & Benvenuto 2008; Inaba et al. 2001). It represents a fast accretion rate. Kenyon & Bromley (2009), applying such a fragmentation-driven accretion scenario, calculate that the core formation process can be completed within 106 yr.

How would gas drag affect these conclusions; i.e., does the rather large accretion rate of Eq. (1) also materialize in the presence of gas drag? Qualitatively, two directions can be envisioned. On the one hand, the dissipative nature of the drag will enhance the collision (impact) radius, like in the case of a dense atmosphere. Conversely, strong particle-gas coupling will suppress the accretion efficiency since the gas after all is not accreted but flows past the object (until the point where it has become more massive than 10 $M_\oplus$ and gas runaway accretion kicks in). It is a priori unclear which aspect of the drag - the coupling or the dissipation - will turn out to be the more important.

To address these questions we include gas drag as an additional force to the restricted 3-body problem that has been previously used in calculating accretion rates in gas-free systems (or in systems where gas can be neglected; Ida & Nakazawa 1989; Petit & Henon 1986). Using appropriate scaling behavior, we show, in Sect. 2, that the system of equations containing all the physics can be restated into two dimensionless parameters: the dimensionless headwind velocity $\zeta _{\rm w}$ that the protoplanet experiences and the dimensionless stopping time (Stokes number, ${\rm St}$) of the particle. Our setup is idealized in the sense that we assume a steady gas flow of constant density (i.e., no pressure fluctuations or atmospheres), a drag law linear in velocity (applicable to small particles), and only consider drift motions of particles.

After having outlined our setup in Sect. 2, Sect. 3 considers the geometrical limit, in which the 2 body interaction is absent or can be ignored. In Sect. 4 we perform an extensive parameter study to obtain the impact parameters as function of the relevant dimensionless quantities. Section 5 presents an analytic model to obtain the impact radii and accretion rates from first principles, which we compare to our measured values. Section 6 illustrates the significance of our result by calculating the protoplanet growth timescale in which we apply a correction to account for the scaleheight of the particle layer. We discuss limitations of our results and summarize in Sect. 7.

2 Sketch of problem and approach

2.1 Definition of impact radius

In this study we will calculate both numerically and analytically impact radii, $b_\sigma $. In 3D systems, the collision rate ${\rm d}M/{\rm d}t$ is the product of the velocity at which the bodies approach each other, the approach velocity, $v_{\rm a}$, the cross section for collisions, $\sigma$, and the volume density in solids $\rho$ that are accreted, ${\rm d}M/{\rm d}t = \rho \sigma v_{\rm a}$. In 2D configurations the vertical dimension is lacking and we define

\begin{displaymath}\left(\frac{{\rm d}M}{{\rm d}t}\right)_{\rm 2D} \equiv 2b_\sigma v_{\rm a} \Sigma = P_{\rm col} \Sigma,
\end{displaymath} (3)

where $P_{\rm col}\equiv2b_\sigma v_{\rm a}$ is the specific collision rate (cf. Nakazawa et al. 1989). In the drag-free regime we indicate the impact radius $b_\sigma $ by  $b_{\rm gf}$. Although we primarily focus on 2D-configurations, Sect. 6 considers a 3D extension in which we apply the derived $b_\sigma $ also for the vertical dimension.

In the gas-free regime particles enter the Hill sphere from orbits both interior and exterior to that of the planet, see Fig. 1. Therefore, $b_{\rm gf}$ is associated with the lengthscale over which particles impact for one of these branches. However, particles can only enter at specific intervals, $1.7R_{\rm h}<\vert b\vert<2.5R_{\rm h}$ (e.g., Greenberg et al. 1991); particles on impact parameters $\vert b\vert<1.7R_{\rm h}$ move on horseshoe orbits that do not enter the Hill sphere. The approach velocity $v_{\rm a}$ for the 3-body regime is defined as the average shear velocity ( $3b\Omega/2$) over the above interval, i.e.,

\begin{displaymath}v_{\rm a} \equiv \frac{1}{2.5R_{\rm h} - 1.7R_{\rm h}}\int_{1...
...^{2.5R_{\rm h}} \frac{3\Omega b}{2}\ {\rm d}b = 3.2 v_{\rm h}.
\end{displaymath} (4)

Using $v_{\rm a} = 3.2v_{\rm h}$ and equating Eq. (1) with Eq. (3) gives

\begin{displaymath}b_{\rm gf} = 1.7 \alpha_{\rm p}^{1/2} R_{\rm h}
\end{displaymath} (5)

as the (effective) impact radius for accretion in the 2D gas-free regime. Note that since $\alpha_{\rm p} \ll 1$, $b_{\rm gf} \ll R_{\rm h}$, which signifies that not every particle that enters the Hill sphere will collide. For this reason we distinguish between $b_{\rm app}=2.5$, the impact parameter at which particles approach (which is related to $v_{\rm a}$), and $b_\sigma $, the impact radius that enters in the expression for the collision rate $P_{\rm col}$. The fact that $b_\sigma\neq b_{\rm app}$ is peculiar to the three-body regime, where the gravity of the central star becomes important.
\begin{figure}
\par\includegraphics[width=8cm]{14903fg1}
\end{figure} Figure 1:

Sketch of particle trajectories in the comoving frame. We consider the motion of the third (test) particle m in the comoving frame of the second body ($M_{\rm p}$, the planet) while including the gravity of the central star. In the gas-free limit zero-eccentricity particles can enter the Hill sphere from both the first and the third quadrant (black curves) but only from specific impact parameters indicated by the hatched regions. Particles arriving at closer impact parameters move on horseshoe orbits. The magnitude and direction of the gas velocity $\vec{v}_{\rm gas}$ as seen from the comoving frame is indicated by the dashed arrows. Particle trajectories including gas drag (solid gray arrows) can be anything depending on the properties of the particle and the gas.

Open with DEXTER

2.2 The circularly restricted three body problem modified by gas drag

We briefly review the circularly restricted three body problem using the framework of Hill's equations (Hill 1878) and include a drag term. The restricted three body problem assumes that the mass of the third body (M3) can be neglected with respect to the masses of the other two bodies, M1, M2. Furthermore, it is assumed that the orbits of the bodies are confined to a single plane and that these are circular for the two massive bodies. In our case the first body is the central star ($M_\star$), the second the (proto)planet ($M_{\rm p}$), and the third the (test) particle m. We then consider the motion of the test particle in a coordinate system centered on and rotating with the motion of the planet. The resulting equations of motions for m in such a frame rotating with angular frequency $\Omega_0$ read

\begin{displaymath}\frac{{\rm d}\vec{v}}{{\rm d}t} = \vec{F} -2\Omega_0 \times \vec{v} + \Omega_0^2 \vec{r},
\end{displaymath} (6)

where $\vec{r} = (x,y,z)$ are the coordinates in the comoving frame, $2\Omega_0 \times \vec{v}$ is the Coriolis acceleration and $\Omega_0^2 \vec{r}$ the centrifugal acceleration. The force per unit mass, $\vec{F}$, acting on the third body consist of the solar gravity, $\vec{F}_{\rm sun} = \Omega^2 \vec{a}$, the 2-body force with the protoplanet, $\vec{F}_{\rm 2b} = GM\vec{r}/r^3$, and the drag force with the gas  $\vec{F}_{\rm drag}$. Expanding the $\vec{F}_{\rm sun}$ term around a0 enables us to linearize Eq. (6) to obtain

\begin{displaymath}\frac{{\rm d}\vec{v}}{{\rm d}t} = \left( \begin{array}{c} 2\O...
...y}{c} x \\ y \\ z \\ \end{array} \right) + \vec{F}_{\rm drag}.
\end{displaymath} (7)

Next, we rewrite Eq. (7) in dimensionless form by normalizing lengths to Hill radii $R_{\rm h}$ (see Eq. (2)) and times to $\Omega_0^{-1}$. The unit of velocity is then the Hill velocity, $v_{\rm h} = R_{\rm h}\Omega_0$. It can be shown that in Hill units $GM = 3\Omega_0^2 R_{\rm h}^3 = 3$. Dropping the z-term as we will treat planar configurations only, Eq. (7), in Hill units, reads

\begin{displaymath}\frac{{\rm d}\vec{v'}}{{\rm d}t'} = \left( \begin{array}{c} 2...
...r'^3 \\ 3y'/r'^3 \\ \end{array} \right) + \vec{F}'_{\rm drag},
\end{displaymath} (8)

where $\vec{F}'_{\rm drag}$ is related to $\vec{F}_{\rm drag}$ as $\vec{F}'_{\rm drag} = \vec{F}_{\rm drag}/R_{\rm h}^2 \Omega_0$.

2.3 The gas drag force

The drag force, $\vec{F}_{\rm drag}$, can be expressed in terms of a stopping time $t_{\rm s}$,

\begin{displaymath}\vec{F}_{\rm drag} = -\frac{\Delta \vec{v}}{t_{\rm s}} = -\frac{(\vec{v}-\vec{v}_{\rm gas})}{t_{\rm s}},
\end{displaymath} (9)

where $\vec{v}_{\rm gas}$ is the velocity of the gas in the comoving frame and $\Delta \vec{v}$ the velocity difference between that of the particle and the gas, see Fig. 1. Due to pressure support, the gas rotates slower than Keplerian by a magnitude $v_{\rm hw} = \eta v_{\rm K}$, where $v_{\rm K}$ is the Keplerian velocity at disk radius a and $\eta$ a dimensionless quantity that gives the fractional deviation from the Keplerian motion (Nakagawa et al. 1986):

\begin{displaymath}\eta = \frac{{\rm d}P/{\rm d}a}{2a\Omega^2 \rho_{\rm g}} \sim \left( \frac{c_{\rm g}}{v_{\rm K}} \right)^2,
\end{displaymath} (10)

with $c_{\rm g}$ the sound speed. In the comoving frame the headwind is directed towards negative y. However, we should correct for the Keplerian shear; thus,

\begin{displaymath}\vec{v}_{\rm gas} = \left(-v_{\rm hw} - \frac{3}{2}\Omega x\right) \vec{e}_y,
\end{displaymath} (11)

where $\vec{e}_y$ is the unit vector in the y direction.

For the drag force we consider several regimes. The stopping time for solid spheres of internal density $\rho_{\rm s}$ for particles of increasing size s reads (Weidenschilling 1977a):

\begin{displaymath}t_{\rm s} = \left\{\begin{array}{ll}
\displaystyle
\frac{\r...
...rm gas}\vert} & {\rm (Quadratic~drag)} \\
\end{array}\right.
\end{displaymath} (12)

where $\rho_{\rm g}$ is the density of the gas, and $\ell_{\rm mfp}$ the mean free path of the gas. For small particles the Epstein regime holds. The Stokes regime supersedes the Epstein regime for particle sizes $s>9\ell_{\rm mfp}/4$. In both the Epstein and the Stokes regime the gas drag is linear with velocity and the stopping time reflects a particle property. These are the regimes for which our study is applicable. In the quadratic regime the stopping time becomes a function of the particle velocity since here $F_{\rm drag}\propto \vert\Delta\vec{v}\vert^2$. In fact, there is a transition regime between the Stokes and quadratic drag regimes where stopping times are proportional to $\vert\Delta \vec{v}\vert^{0.4}$, which we have, for reasons of simplicity, ignored here (following Rafikov 2004).

As an (approximate) upper limit for $\vert\Delta \vec{v}\vert$ we can take the headwind velocity, $v_{\rm hw}$. The transition between the Stokes and the quadratic drag regimes then occurs at a size of

                                     s $\textstyle \approx$ $\displaystyle s_{\rm max} = \frac{27 \ell_{\rm mfp} c_{\rm g}}{2 v_{\rm hw}}$ (13)
  = $\displaystyle \ 90\ {\rm m} \left( \frac{c_{\rm g}}{10^5\ {\rm cm\ s^{-1}}} \ri...
...{-3}}} \right)^{-1} \left( \frac{v_{\rm hw}}{30\ {\rm m\ s^{-1}}} \right)^{-1},$  

where we used $\ell_{\rm mfp} = 2\times10^{-9}/\rho_{\rm g}$ (in cgs units; Nakagawa et al. 1986). Since we consider a drag law that is linear in velocity, our results are only applicable for particle sizes less than $s_{\rm max}$. In the inner disk (where the gas density is large) the results should be applicable to the sweepup of chondrule-like particles and m-size boulders. In the outer disk, $\rho_{\rm g}$ is much lower and the particles for which our results are applicable include (small) planetesimals.

Expressed in dimensionless form the drag law reads

\begin{displaymath}\vec{F}_{\rm drag}' = \frac{\vec{F}_{\rm drag}}{v_{\rm h}\Ome...
...c} -v_x' \\ -v_y' -\zeta_{\rm w} -3x'/2 \\ \end{array} \right)
\end{displaymath} (14)

where we used Eq. (9) for $\vec{F}_{\rm drag}$ and Eq. (11) for $\vec{v}_{\rm gas}$, normalized velocities to $v_{\rm h}$, and have introduced the Stokes number,  ${\rm St}=t_{\rm s} \Omega_0$[*], and the dimensionless headwind velocity,  $\zeta _{\rm w}$:
                            $\displaystyle \zeta_{\rm w}$ $\textstyle \equiv$ $\displaystyle \frac{v_{\rm hw}}{\Omega R_{\rm h}} \approx 12.5 \left( \frac{\rh...
...\ cm^{-3}}}\right)^{-1/3} \left( \frac{v_{\rm hw}}{30\ {\rm m\ s}^{-1}} \right)$  
    $\displaystyle \times\left( \frac{R_{\rm p}}{100\ {\rm km}} \right)^{-1} \left( \frac{a}{1\ {\rm AU}} \right)^{1/2},$ (15)

with $R_{\rm p}$ the radius of the (proto)planet. Note that due to the normalization to $v_{\rm h}$, $\zeta _{\rm w}$ is primarily an indicator of the size of the (proto)planet rather than of the strength of the headwind $v_{\rm hw}$ as the latter is approximately constant throughout the disk.

2.4 Dimensionless quantities

Table 1:   Dimensional and dimensionless parameters.

Table 1 compiles some key quantities in both dimensional and dimensionless form. These include the dimensionless headwind velocity $\zeta _{\rm w}$ (Eq. (15)), the particle Stokes number ${\rm St}$, and the protoplanet radius

\begin{displaymath}\alpha_{\rm p} = \frac{R_{\rm p}}{R_{\rm h}} = 5.7\times10^{-...
...-3}}} \right)^{-1/3} \left( \frac{a_0}{{\rm AU}} \right)^{-1},
\end{displaymath} (16)

which mainly depends on semi-major axis a0. In Fig. 2 we give the relation between the dimensionless $\zeta _{\rm w}$ and ${\rm St}$ to the physical protoplanet size $R_{\rm p}$ and particle size s for two disk radii a. We have adopted $v_{\rm hw} = 30\ {\rm m\ s^{-1}}$ and disk parameters that correspond (approximately) to a typical minimum mass solar nebula model (Weidenschilling 1977b; Hayashi et al. 1985).

The full set of equations of motions in dimensionless form, dropping the primes, reads

\begin{displaymath}\frac{{\rm d}}{{\rm d}t}\left( \begin{array}{c} v_x \\ v_y \\...
...rray}{c} v_x \\ v_y+\zeta_{\rm w}+3x/2 \\ \end{array} \right).
\end{displaymath} (17)

The drag-free equations of motions are retrieved when ${\rm St}\rightarrow \infty$, which signifies that particles are not coupled to the gas. However, for Stokes number ${\rm St} \la 1$ particles are coupled to the gas and the importance of the drag terms becomes relevant or even dominant. Petit & Henon (1986) and Ida & Nakazawa (1989), working in the gas-free regime, only had to care about the first terms on the RHS of Eq. (17) and the equations of motions did not include any parameter. The addition of gas drag introduces two parameters: the velocity of the gas flow $\zeta _{\rm w}$ and the coupling parameter ${\rm St}$. Together with the size of the planet, $\alpha_{\rm p}$, these fully specify the problem; i.e., impact parameters $b_\sigma $ depend on these three dimensionless quantities only. Although not as clean as the drag-free equations, Eq. (17) still represent a significant reduction of the parameters involved (semi-major axis, particle size s, protoplanet size $R_{\rm p}$, headwind velocity, gas density, etc.). In our parameter study we only have to care about these three parameters.

We do not include the eccentricity in our prescription (and also not the inclination since the interaction is assumed to be 2D). Rather, the initial velocity of the approaching particle is given by the radial drift equations for individual particles, neglecting the 2-body interaction term. These we will now review.

3 The geometrical limit[*]

Ignoring the 2-body interaction force, we will analytically solve for the particle's trajectory in the comoving frame. As we will soon see, impact parameters along a particle trajectory are generally not conserved. We provide a general relation between the impact parameter at the interaction point ($b_\sigma $) and its projected value on the x and y axes at any arbitrary point (Eq. (21)). This relation will be used later in Sect. 4 to obtain the impact radii $b_\sigma $.

3.1 Steady-state velocities

Without the two-body interaction term, the motion of the particle fulfills the well-known drift equations (Weidenschilling 1977a; Nakagawa et al. 1986; Brauer et al. 2007):

\begin{displaymath}
v_{\rm r} = - \frac{2v_{\rm hw} {\rm St}}{1+{\rm St}^2};
\end{displaymath} (18a)

\begin{displaymath}v_\phi = - \frac{v_{\rm hw}}{1+{\rm St}^2},
\end{displaymath} (18b)

where $v_{\rm r}$ is the radial velocity and $v_\phi$ the azimuthal velocity with respect to the local Keplerian rotation. Thus, in the context of a fixed-rotating orbital frame we have $v_x = v_{\rm r}$ and $v_y = v_\phi - \frac{3}{2}\Omega_0 x$, or in dimensionless units (divide by $v_{\rm h}$)

\begin{displaymath}v_x = - \frac{2\zeta_{\rm w} {\rm St}}{1+{\rm St}^2}
\end{displaymath} (19a)

\begin{displaymath}v_y = - \frac{\zeta_{\rm w}}{1+{\rm St}^2} - \frac{3}{2}x
\end{displaymath} (19b)

and it can be verified that with these expressions the RHS of Eq. (17) vanishes when the two body interaction terms (-3x/r3 and -3y/r3) are omitted.

3.2 The parabola solution

\begin{figure}
\par\includegraphics[width=8.5cm]{14903fg2}
\end{figure} Figure 2:

Relation between dimensionless and physical quantities. The Stokes number (${\rm St}$, solid lines) and the dimensionless headwind velocity $\zeta _{\rm w}$ (dashed lines) are plotted on the y-axis as function of the radius s of the particle and the radius $R_{\rm p}$ of the protoplanet. Note the different units of s and $R_{\rm p}$ on the x-axis. The black dot denotes $s_{\rm max}$ (Eq. (13)). Lines are shown for: (i) $c_{\rm g}=10^5~{\rm cm\ s^{-1}}$ and $\rho_{\rm g}=10^{-9}~{\rm g~cm^{-3}}$ at a disk radius of 1 AU (black lines) and (ii) $c_{\rm g} = 6\times10^4~{\rm cm~s^{-1}}$ and $\rho_{\rm g} = 10^{-11}~{\rm g~cm^{-3}}$ at a position of 10 AU (gray lines). The internal density of solids is fixed at $\rho_{\rm s} = 3~{\rm g~cm^{-3}}$, the nebula headwind is $v_{\rm hw}=30~{\rm m~s^{-1}}$, and the mass of the central object is solar.

Open with DEXTER
\begin{figure}
\par\includegraphics[width=8.5cm]{14903fg3}
\end{figure} Figure 3:

Without the two body force, particle trajectories as witnessed from the comoving frame obey parabolas. Two trajectories are shown: one that passes through the origin (y0(x)) and one that just hits the target. The corresponding impact parameter b is denoted by arrows. For curved trajectories b is not conserved due to the changing slope of the curves, here indicated by the angle $\theta $. At S we have that $\Delta x_{\rm S} = b_{\rm S}/\!\sin \theta_{\rm S}$ and $\Delta y_{\rm S} = b_{\rm S}/\!\cos \theta_{\rm S}$. $\Delta y_{\rm S}$ is a conserved quantity.

Open with DEXTER
Since

\begin{displaymath}\frac{{\rm d}y}{{\rm d}x} = \frac{v_y}{v_x} = \frac{1}{2{\rm St}} + \frac{3(1+{\rm St}^2)}{4{\rm St}\zeta_{\rm w}}x
\end{displaymath} (20)

we immediately recognize that the particle's trajectory in the rotating frame obeys a parabola, y(x)=Ax2 +Bx +C with  $A =3(1+{\rm St}^2)/8{\rm St}\zeta_{\rm w}$, $B=1/2{\rm St}$, and C the integration constant, which is determined by the starting point S of the particle. We refer to the function that intersects the origin (C=0) as y0(x). The starting point $(x_{0S}, y_{\rm S})$ is defined to lie on this curve, $y_{\rm S} = y_0(x_{0S})$, see Fig. 3, where y0(x) is plotted by the upper parabola. Another parabola with the same A and B (i.e., for the same particle properties ${\rm St}$ and $\zeta _{\rm w}$) but with non-zero C is drawn in such a way that it just hits the ``target'' at the origin. Its ``launch point'' at $y=y_{\rm S}$ is shifted over a length $\Delta x_{\rm S}$. The vertical difference between the curves, $-C=\Delta y_{\rm S}$, is preserved.

If the particle trajectories were straight, impact parameters would be the same everywhere in the (x,y)-plane. However, due to the x2-term this statement no longer holds for the general case of nonzero A. See Fig. 3: the impact parameter near the origin or at the interaction region, $b_\sigma $, differs from that at S. Impact parameters are no longer conserved due to the change in ${\rm d}y/{\rm d}x$. The changing slope of the curve is indicated in Fig. 3 by the angle $\theta $, where $\theta $ is related to Eq. (20) as $\theta = \arctan ({\rm d}y/{\rm d}x)$.

Using the properties of the parabola solution we can relate the quantities at I to those at S. At the interaction point I the impact parameter is $b_\sigma $ (= $\alpha_{\rm p}$ in the geometrical case) and the associated vertical width is $\Delta y_I = b_\sigma/\cos \theta_I$ with $\theta_I$ the angle the parabola makes at this point. Similarly, $\Delta y_{\rm S} = b_{\rm S}/\cos \theta_{\rm S}$ and due to the invariance of $\Delta y$ we therefore have that $b_\sigma = b_{\rm S} \cos \theta_I/\cos \theta_{\rm S}$.

The associated change in x at the starting point is $\Delta x_{\rm S} = \Delta y_{\rm S}/\tan \theta_{\rm S}$. This can be expressed in terms of the impact parameter at the interaction point $b_\sigma $, i.e.,

\begin{displaymath}\Delta x_{\rm S} = \frac{\Delta y_{\rm S}}{\tan \theta_{\rm S...
...rm d}y/{\rm d}x \right)^2_I }}{({\rm d}y/{\rm d}x)^2_{\rm S}},
\end{displaymath} (21)

where we used that $1/\cos \theta = 1/ \cos [\arctan ({\rm d}y/{\rm d}x)] = \sqrt{1+({\rm d}y/{\rm d}x)^2}$. Of course, in the non-gravity limit we know already that the impact parameter is just the planet radius, $b_\sigma=\alpha_{\rm p}$. However, Eq. (21) is essential to interpret the numerical result of Sect. 4. That is, in our numerical integration (that includes the 2-body force) we will scan the x-axis for trajectories that lead to a collision and obtain a length scale $\Delta x_{\rm S}$ over which particles hit the target. Using the above equation we can then relate the obtained range of projected impact parameters at S, $\Delta x_{\rm S}$, to the impact parameter at the interaction point, $b_\sigma $.

3.3 Collision rates

The collision rate P in the 2D configuration is the product of the collision cross section, $2b_\sigma$, and the approach velocity, $v_{\rm a}=\sqrt{v_x^2+v_y^2}$. At S we have that $b_{\rm S} v_{\rm S} = \Delta y_{\rm S} v_x = \Delta x_{\rm S} v_y$. At I the collision rate equals $b_\sigma v_{\rm a} = \alpha_{\rm p} v_{\rm a}$. Now, using Eq. (21) we have that $b_\sigma v_I = b_\sigma v_x/\cos \theta_I = v_x \tan \theta_{\rm S} \Delta x_{\rm S} = v_x \Delta y_{\rm S}$ and we see that the rates at I and S are equal and independent of the choice for the starting point $y_{\rm S}$. This result just reflects mass conservation. Thus, we obtain the collision rate:

\begin{displaymath}P_{\rm geo} = 4\alpha_{\rm p} \frac{\zeta_{\rm w} {\rm St}}{(...
...\rm St}^2) +4\zeta_{\rm w})^2}{64{\rm St}^2 \zeta_{\rm w}^2}},
\end{displaymath} (22)

where $({\rm d}y/{\rm d}x)$ has been evaluated at $x=\alpha_{\rm p}$. The complexity of Eq. (22) may seem surprising for something as straightforward as a geometrical sweepup. However, this is entirely due to the fact that Eq. (22) covers several (velocity) regimes. In Appendix A we consider the asymptotic limits of Eq. (22) and show that these correspond to the expected sweepup rates (cross section $\times$ approach velocity) and to the findings of Kary et al. (1993).

4 Full 3-body integrations including gas drag and gravity

4.1 Description of the adopted algorithm

We perform a parameter study of Eq. (17), varying $\zeta _{\rm w}$ and ${\rm St}$. For the dimensionless headwind velocity runs were performed at $\zeta_{\rm w} = 0.01, 0.03, 0.1,\dots10^4$ and for the Stokes number values of ${\rm St}=10^{-4}, 3\times10^{-4}\dots10^4$ were sampled. Thus, we obtain a grid of $17\times13=221$ different combination of $\zeta _{\rm w}$ and ${\rm St}$. Not every combination is equally likely. Indeed, the parameter space samples areas where our key approximations (linear drag law, constant gas density) lose validity, but we intentionally sample a broad range of values to verify the validity of our analytical expressions (see Sect. 5).

\begin{figure}
\par\includegraphics[width=9cm]{14903fg4}
\end{figure} Figure 4:

The minimum distance in units of Hill radii to the origin (center of the protoplanet), $r_{\rm min}$, as function of the x-coordinate of the starting point, $x_{\rm S}$ ($y_{\rm S}$ is fixed at 40 Hill radii). Plotted is $r_{\rm min}$ for a gas-free system (left) and a system that is characterized by the parameters ${\rm St}=10$ and $\zeta _{\rm w}=1$ (right). Several bands are labeled. The inclusion of gas drag shifts the bands to larger $x_{\rm S}$ while merging several chaotic features within the chaotic c-band.

Open with DEXTER

For each combination of $\zeta _{\rm w}$ and ${\rm St}$ we numerically determine the function $r_{\rm min}(x_{\rm S})$. Particles are launched from a starting point $(x_{\rm S},y_{\rm S})$ where $y_{\rm S}$ is fixed and $x_{\rm S}$ is varied, see Fig. 3. The initial velocities are given by Eq. ([*]). Depending on the sign of $v_y(x_{\rm S})$ the initial y-position (+ or -) is determined, such that the initial motion in y is always directed towards the planet. Here, we fix $\vert y_{\rm S}\vert$ at 40 (Ida & Nakazawa 1989). Particles that leave the computational domain (when |y|>40 or x<-40) are no longer followed. For a certain $x_{\rm S}$ we numerically integrate Eq. (17) adopting a relative error of at most 10-8. As our integrator we use a fifth-order Runge-Kutta scheme with timestep control (Shampine et al. 1979; Fehlberg 1969). After the calculation has terminated we determine (and store) the minimum distance, $r_{\rm min}$. In this way $r_{\rm min}(x_{\rm S})$ is obtained, see Fig. 4. Projected impact parameters $\Delta x_{\rm S}$ are then obtained from the $r_{\rm min}(x_{\rm S})$ curve by summation over the intervals where $r_{\rm min} < \alpha_{\rm p}$, i.e.,

\begin{displaymath}\Delta x_{\rm S} = \int {\rm d}x_{\rm S} H( \alpha_{\rm p} - r_{\rm min}[x_{\rm S}]),
\end{displaymath} (23)

where H(t) is the Heaviside step function,

\begin{displaymath}H(t) = \left\{\begin{array}{cc}
1 & t \ge 0; \\
0 & t<0.
\end{array}\right.
\end{displaymath} (24)

(Simply put: we only include the orbits that hit the target.)

As $r_{\rm min}$ is occasionally found to vary steeply with $x_{\rm S}$, fine sampling of the x-axis is required. Therefore, we sample our parameters space ($x_{\rm S}$) adaptively. We start out with intervals of 1 Hill radii, e.g.,  $\{x_{\rm S}\} = 0, 1, 2, 3, \dots$. In the next level the interval spacing is reduced by a factor 10, $\delta = 0.1$. However, we only treat the points that fulfill the condition $r_{\rm min}(x_{\rm S}) < F \delta$, where F is empirically fixed at 103 (see below). For example, if $r_{\rm min}(0) = 300$ this point will be skipped in the next iteration of the algorithm. If $x_{\rm S}$ fulfills the condition, however, then both sides will be scanned; e.g., if $r_{\rm min}(2) = 15 < 0.1F$ then, we will additionally perform calculations for $x_{\rm S}=1.1,1.2,\dots 1.9$ and $x_{\rm S}=2.1,2.2,\dots2.9$. In this way we reduce the number of calculations but are still able to obtain a good assessment of $\Delta x_{\rm S}$ for low $\alpha_{\rm p}$.

Despite this optimization, we were forced to perform a relatively large number of integrations, i.e., a large F. The reason is the presence of very narrow, chaotic bands. In Fig. 4 band b near $x_{\rm S}=2.0$ is a regular band since $r_{\rm min}$ varies smoothly with $x_{\rm S}$. A low F value suffices to pick up this feature. However, bands a and d are very narrow and show (if one would zoom in) additional substructure. These chaotic bands are not resolved when choosing a low F. In fact, there is no guarantee that our algorithm will pick up every band since they can be very narrow. However, with F=103 we do obtain a good correspondence to previous works (Ida & Nakazawa 1989; Petit & Henon 1986), also matching the substructure within the narrow bands shown in Fig. 4.

Following the discussion in Sect. 3 we emphasize again that the starting points ($x_{\rm S}$ values) are not fundamental, but depend on the choice of the starting point $y_{\rm S}$. Taking a different value of $y_{\rm S}$, e.g., $y_{\rm S}=80$, will shift the features of Fig. 4b towards higher $x_{\rm S}$ values. In addition, the spacing (width of the features) will be different. The only conserved (physical) quantity is the mass flux, i.e., the integral of $\int v_y(x_{\rm S}){\rm d}x_{\rm S}$ over the width of the feature - independent of the choice of $y_{\rm S}$.

\begin{figure}
\par\includegraphics[width=17cm]{14903fg5}
\end{figure} Figure 5:

Examples of planet-particle interactions for different values of the dimensionless headwind velocity $\zeta _{\rm w}$ and coupling parameter ${\rm St}$. For typical nebula parameters particles of ${\rm St}=10$ correspond to loosely coupled m-size particles, whereas ${\rm St}=0.01$ are more strongly coupled cm-size particles, see Fig. 2. Likewise, $\zeta _{\rm w}=1$ corresponds to protoplanets of $R_{\rm p}\sim 10^3$ km in radius, while $\zeta _{\rm w}=100$ corresponds to $R_{\rm p}\sim 10$ km planetesimals. (A) Two particles of ${\rm St}=10$ experience a close encounter within the Hill sphere (dotted circle). The $x_{\rm S}=3.9$ particle is captures and settles to the planet, whereas the other particle is ejected from the Hill sphere (The Keplerian shear eventually causes it to resurface at the other side of the Hill sphere). (B) Strong gas coupling, ${\rm St}=0.01$. There is a competition between the gravitational pull of the planet and the drag force directed towards negative y. (C) Close encounters at large $\zeta _{\rm w}$ without settling (see inset). (D, E) Examples of particle trajectories originating from interior orbits. (F) Radially approaching orbits. (This figure is available in color in electronic form.)

Open with DEXTER

4.2 Orbits including gas drag

Figure 5 provides several examples of particle trajectories that experience gas drag. Figure 5a shows the trajectories for $\zeta _{\rm w}=1$ and ${\rm St}=10$ with different starting points $x_{\rm S}$ (Fig. 4 contains the same parameters). For a ``standard'' nebula setting, these parameters correspond to $\sim$m-size particles accreting onto a $\approx$103 km planet, see Fig. 2. Due to the large Stokes number, the influence of the gas is relatively weak and the orbits bear a close resemblance to the gas-free, three body regime. The $x_{\rm S}=3.8$ trajectory enters the Hill sphere, where it experiences a close encounter, at $r_{\rm min}=2.8\times10^{-2}$, before leaving the Hill sphere. It then re-emerges later at negative x due to the combined effects of inwards radial drift and Keplerian shear. However, the particle that started out at $x_{\rm S}=3.9$ is captured within the Hill sphere and experiences strong orbital decay due to gas drag.

Figure 5b shows orbits for smaller particles of Stokes number 0.01. Four orbits are shown of which two lead to accretion. Clearly, there is a contest between the gravitational pull of the planet and the aerodynamic pull of the gas flow. Once close enough, gravity always wins. All orbits with $0.38\le x_{\rm S} \le 0.74$ are accreted; there are no close encounters. Accretion is independent of the physical proportion or internal density of the protoplanet; once a particle's angular velocity about the planet is damped by drag, it settles radially at its terminal velocity. The only relevant physical quantity is the mass. This mode of accretion reflects the capture mechanism of Fig. 5a. We will refer to orbits like the $x_{\rm S}=3.9$ curve in Fig. 5a as gas drag induced orbital decay, whereas the accretion mode in Fig. 5b is referred to as settling and draw the dividing line at ${\rm St}=1$.

On the other hand Fig. 5c, which features a larger dimensionless headwind (meaning: a smaller protoplanet) of $\zeta _{\rm w}=100$, does not display the settling behavior. Here, particles can only be accreted due to the finite size of the target. The $x_{\rm S}=0.796$ trajectory has a minimum distance of $r_{\rm min}=5.0\times10^{-4}$; the $x_{\rm S}=0.8$ trajectory $r_{\rm min}=4.5\times10^{-3}$. Clearly, for a planet size $\alpha_{\rm p} \ll 1$ the impact parameter in Fig. 5c is much less than for the settling orbits of Fig. 5b. Since the Stokes numbers are the same, the reason must be due to the larger headwind velocity $\zeta _{\rm w}$. This is understandable since particles of ${\rm St}\ll1$ approach at the headwind velocity ( $v_{\rm a} \approx \zeta_{\rm w}$) and a large $v_{\rm a}$ is not conducive for accretion.

In the lower panels of Fig. 5 we vary either the Stokes number (particle size) or $\zeta _{\rm w}$ (protoplanet size) with respect to the panel above. For a Stokes number of 103, see Fig. 5d, the effects of gas-drag are even less pronounced and it becomes more difficult to capture these (big) particles within the Hill sphere. Moreover, if such a particle would be captured, it takes longer to finally accrete this particle due to orbital decay. Another difference with Fig. 5a is that the ${\rm St}=10^3$ particles can now also enter the Hill sphere from interior orbits (negative $x_{\rm S}$). In Fig. 5a the strong radial drift still prevents particles from entering the Hill sphere from the negative y-direction; however, for ${\rm St}=10^3$ the radial drift is sufficiently reduced to render the situation more akin to the symmetric gas-free limit.

The $\zeta_{\rm w}=0.1$ orbits in Fig. 5e also feature accretion from particles approaching the planet from interior orbits, which the $\zeta_{\rm w}=1.0$ orbits of Fig. 5b were not capable of. The dimensionless headwind parameter of $\zeta_{\rm w}=0.1$ corresponds to a very big planet (in the canonical model) for which, as we will discuss below, the constant gas density background is unrealistic. Alternatively, it can represent a smaller protoplanet in a nebula where the dimensional headwind is, for some reason, strongly reduced. In any case, we see that low $\xi_{\rm w}$ tends to make the interactions more symmetric. This can be seen from the % latex2html id marker 4011
$\mbox{Eq.~(\ref{eq:vy})}$: low $\zeta _{\rm w}$ or large ${\rm St}$ reduce the contribution from the non-symmetric headwind term, $\zeta_{\rm w}/(1+{\rm St}^2)$.

Figure 5f shows, however, that for ${\rm St}=10$ and $\zeta _{\rm w}=100$ the picture is anything but symmetric. The particles approach the planet from a very radial direction (x-direction) - at least, as seen from the perspective of the planet. The point is here that both ${\rm St}$ and $\zeta _{\rm w}$ are large. Thus, both planet and particle move at a Keplerian velocity (in the azimuthal direction) but, due to the large $\zeta _{\rm w}$, the particle still suffers a significant radial drift, which outweighs the effects of the Keplerian shear. As a result, the situation is similar to Fig. 5c: accretion does only proceed through close encounters.

4.3 Collision rates

We obtain the (dimensionless) collision rate from the encounters that hit the protoplanet, i.e.,

\begin{displaymath}P(\alpha_{\rm p}, {\rm St}, \zeta_{\rm w}) = 2\int {\rm d}x_{...
..._y(x_{\rm S})\vert H(\alpha_{\rm p} - r_{\rm min}[x_{\rm S}]),
\end{displaymath} (25)

vy(x) given by Eq. (19) and H(t) the Heaviside step function. Figure 6 plots contours of $P(\zeta_{\rm w}, {\rm St})$ for a planet size of $\alpha _{\rm p}=10^{-3}$, which corresponds to an heliocentric distance of $\approx$5 AU. The reader must realize that P is expressed in dimensionless units; the large rates that can be seen at large $\zeta _{\rm w}$ (and ${\rm St}<1$) are less impressive upon multiplication by $R_{\rm h} v_{\rm h} \propto \zeta_{\rm w}^{-2}$ (see Eq. (15)). In fact, these high P values are consistent with the geometrical sweepup rates of Eq. (22). However, the expression in terms of dimensionless units is useful since we can directly compare it to the gas-free limit for which $P_{\rm gf} \approx 11\alpha_{\rm p}^{1/2} = 0.35$, see Sect. 2.1. For large ${\rm St}$ and small $\zeta _{\rm w}$, P converges to $P_{\rm gf}$, the expected behavior. However, for the remainder P deviates significantly from $P_{\rm gf}$. We sum up the main features:

1.
Particles of ${\rm St}\sim1$ accrete very well when the headwind velocity is low. There is a distinct peak at $({\rm St},\zeta_{\rm w}) = (0,0)$; the accretion rate is here 20 times larger than Eq. (1). However, there is a very sharp transition between $1\la \zeta_{\rm w} \la 10$.
2.
For large $\zeta _{\rm w}$, $P=P_{\rm geo}$ is larger than $P_{\rm gf}$ although no gravitational focusing takes place. The sweepup is so effective due to the strong headwind.
3.
For ${\rm St}>1$, the band ${\rm St} \sim \zeta_{\rm w}$ features a maximum in P.
4.
For ${\rm St}\ll1$ and low $\zeta _{\rm w}$ (large planets), collision rates are lower than $P_{\rm gf}$. Tiny dust particles stay ``glued'' to the gas due to their strong coupling, preventing accretion.
\begin{figure}
\par\includegraphics[width=9cm]{14903fg6}
\end{figure} Figure 6:

Contour plot of collision rates obtained from the numerical integrations for $\alpha _{\rm p}=10^{-3}$. Contours of $\log_{10} P$ are shown as function of the Stokes number and the headwind velocity. Contour levels are indicated except for $\log_{10} P = 0$ (thick solid), -0.25 (dotted) and -0.75 (dashed). The accretion rate in the gas-free regime is $\log_{10} P_{\rm gf} \approx -0.46$. (This figure is available in color in electronic form.)

Open with DEXTER

5 A simple model for gravo-gas interactions

5.1 Model outline

We present a simple model for the impact parameter $b_\sigma $. The model is summarized in Fig. 7, where the three relevant regimes for the impact parameter are shown. In the hyperbolic regime encounters follow the two-body approximation. The usual gravitational focusing formula applies. Keplerian shear is unimportant. In the settling regime particles settle to the target and the impact parameter is independent of the planet size, $\alpha_{\rm p}$. For this reason, impact parameters can become rather large. Finally, in the three-body regime, the encounter proceeds along the lines of the drag-free three body encounters at low energy. However, the presence of the gas now causes some particles to be captured within the Hill sphere; these orbits decay and this enhances the accretion rate.

\begin{figure}
\par\includegraphics[width=9cm]{14903fg7}
\end{figure} Figure 7:

Illustration of the three accretion regimes. In the hyperbolic regime interactions are 2-body encounters and the standard formula for gravitational focusing applies. In the settling regime, accretion proceeds through settling, which enhances the impact parameter $b_\sigma $. Else, for ${\rm St}>{\rm max}(\zeta _{\rm w}, 1)$ the solar gravity affects the encounter and $b_\sigma $ is increased with respect to the gas-free limit ( $b_{\rm gf}$, Eq. (5)) due to particle capture within the Hill sphere. The gray band approximately indicates the zone where the settling solutions permeate into the hyperbolic regime. The letters a-f correspond to the parameters in the panels of Fig. 5.

Open with DEXTER

5.1.1 Importance of three body encounters

It is clear that for ${\rm St}\gg1$ the encounter cannot be described by a 2-body interaction, but should include the solar gravity. But what is the transition between the 2-body and the 3-body regime in the presence of gas drag? A passage through the Hill sphere typically takes a time of the order of the orbital period. Thus, at first sight, we can draw the boundary at ${\rm St}=1$ since particles of lower stopping time will be strongly affected by the gas. However, a large headwind velocity $\zeta _{\rm w}$ has the same effect. Particles that experience a drag force $\zeta_{\rm w}/{\rm St}>1$ are blown out of the Hill sphere. Thus, 3-body effects are reduced to the region of parameter space where ${\rm St}>1$ and ${\rm St}>\zeta_{\rm w}$, see Fig. 7.

5.1.2 Two body regime: settling- and hyperbolic interactions

We consider an interaction at impact parameter $b=b_\sigma$ at an approach velocity $v_{\rm a}$. The strength of the gravitational force is $f_{\rm g} = 3/b^2$ and the interaction timescale, $t_{\rm a} \simeq b/v_{\rm a}$. The latter quantity can be compared to the particle's response time ${\rm St}$. When $t_{\rm a}<{\rm St}$ gas drag can be neglected during the encounter; the induced velocity change is $\Delta v = f_{\rm g} t_{\rm a} = 3b/v_{\rm a}$. However, if $t_{\rm a}>{\rm St}$ the particle's velocity equilibrates towards $\Delta v = \Delta v_{\rm set} = f_{\rm g} {\rm St} = 3{\rm St}/b^2$.

For ${\rm St}<1$ the approach velocity can be approximately written as $v_{\rm a} = 3b/2 +\zeta_{\rm w}$. For low $\zeta _{\rm w}$ we therefore can expect settling since encounter timescales $t_{\rm a}$ are long. However, for large $\zeta _{\rm w}$ settling will be prohibited: either the induced change $\Delta v$ is too little (at large b) or the interaction timescale too short for the particle to obtain its settling velocity (at low b).

To see this quantitatively, the minimum impact parameter for settling is $b^\ast = \zeta_{\rm w} {\rm St}$ and the corresponding velocity change is $\Delta v^\ast = 3 /\zeta_{\rm w}^2 {\rm St}$. In order for the particle to settle to the central object, the direction of the particle has to change over a large angle, i.e., $\Delta v \sim v_{\rm a}$. In fact, we obtain a better correspondence with our numerical result if we set the required velocity change to $v_{\rm a}/4$. Analytically, it can be shown that this is the required change for ${\rm St}\ll1$, see Appendix B. Thus, accretion through settling takes place when $v_{\rm a}/4 \le \Delta v_{\rm set} \le \Delta v^\ast$ and disappears when $v_{\rm a}/4 >\Delta v^\ast$. At the boundary between the settling and hyperbolic regime it is allowed to take $v_{\rm a} = \zeta_{\rm w}$ (as can be verified a posteriori). We then have that for $\zeta_{\rm w}^3 {\rm St} > 12$ settling is no longer possible, corresponding to a critical Stokes number

\begin{displaymath}
{\rm St}^\ast = \frac{12}{\zeta_{\rm w}^3},
\end{displaymath} (26)

above which accretion through settling will no longer occur.

5.1.3 The settling regime

Assuming the settling regime, particles at impact parameter b experience a velocity impulse of $\Delta v_{\rm set} = 3{\rm St}/b^2$, which should equal $v_{\rm a}/4$ for accretion. Since $v_{\rm a} \approx 3b/2 + \zeta_{\rm w}$ in this regime (note that we do not neglect the shear term since it becomes important at low $\zeta _{\rm w}$) the condition $\Delta v_{\rm set} = v_{\rm a}/4$ requires us to solve the cubic equation

\begin{displaymath}b^3 + \frac{2\zeta_{\rm w}}{3}b^2 - 8{\rm St} = 0.
\end{displaymath} (27)

The (real, positive) solution to this equation is denoted $b_{\rm set}$.

5.1.4 Hyperbolic regime

At large $\zeta _{\rm w}$ the encounter is fast and the presence of gas drag can be ignored during the encounter. For accretion we now require (by conservation of angular momentum) that $\Delta v = 3b/v_{\rm a} \ge v_{\rm esc}$ which is much larger than in the settling regime. For the impact radius we can just take the standard expression of the gravitationally-enhanced cross section,

\begin{displaymath}b_{\rm hyp} = \alpha_{\rm p} \sqrt{ 1 + \left( \frac{v_{\rm e...
...alpha_{\rm p} \sqrt{1 + \frac{6}{\alpha_{\rm p} v_{\rm a}^2}},
\end{displaymath} (28)

where $v_{\rm esc} = \sqrt{6/\alpha_{\rm p}}$ in Hill units. For the approach velocity $v_{\rm a}$ we now also include the horizontal velocity components (i.e., vx) since for ${\rm St}>1$ vx becomes dominant. However, it is fine to neglect the shear term in Eq. (
19) since $b_{\rm hyp}\ll1$. Thus, for the approach velocity in the hyperbolic regime we can write

\begin{displaymath}v_{\rm a,hyp} = \zeta_{\rm w} \frac{\sqrt{1 + 4{\rm St}^2}}{1+{\rm St}^2}\cdot
\end{displaymath} (29)

5.1.5 The three body regime

Without gas drag the effective impact parameter for collisions is Eq. (5), $b_{\rm gf} = 1.7 \alpha_{\rm p}^{1/2}$, in dimensionless units (see Sect. 2.1). Gas drag adds another component to the impact parameter on top of Eq. (5). Figure 4 neatly illustrates this behavior. The chaotic band c in the gas-drag simulation has collapsed. Particles entering at the corresponding $x_{\rm S}$-values are captured and decay to the central object on a timescale $\sim$${\rm St}$. If the gas inside the Hill sphere is removed within this timescale, these particles will become satellites; however, here we will simply assume that all captured particles contribute to the collision rate.

Because the accretion in the dissipative 3-body regime is determined by the behavior of the chaotic zones, it is difficult to provide an analytic model for the enhanced $b_\sigma $. The chaotic zones are especially susceptible to collapse, because these particle trajectories are characterized by many revolutions, trough which a lot of energy can be dissipated. In the gas-free situation one requires a positive energy J to enter the Hill sphere,

\begin{displaymath}J = \frac{1}{2}v^2 -\frac{3}{r} -\frac{3}{2}x^2 +\frac{9}{2};
\end{displaymath} (30)

and once J becomes negative in the Hill sphere, e.g., by inelastic collisions, the body becomes trapped (Ohtsuki 1993). Here, we face a similar situation where the gas drag is responsible for the energy removal. Unfortunately, in our case an analysis in terms of the Jacobian is not so meaningful as the gas flow can also add energy; i.e., J is not conserved and bodies with J<0 can still be ``blown out'' of the Hill sphere. However, the picture - that gas drag can trap particles - is still the key.

Empirically, we find that the impact radius is increased by a term $1/{\rm St}$, which corresponds to the dissipated energy over a revolution. For these reasons, we add a term proportional to $1/{\rm St}$ to Eq. (5),

\begin{displaymath}b_{\rm 3b} = b_{\rm gf} + \frac{1.0}{{\rm St}} = 1.7\alpha_{\rm p}^{1/2} + \frac{1.0}{{\rm St}},
\end{displaymath} (31)

where the 1.0 constant is obtained empirically. As we will see in the next section Eq. (31) fits the general trend well, but it cannot reproduce the impact radius at every Stokes value.

5.2 Comparison to numerical results and fine tuning of the recipe

\begin{figure}
\par\includegraphics[width=9cm]{14903fg8}
\end{figure} Figure 8:

Impact radii from the numerical integrations (symbols) and analytic fits (curves) for a headwind velocity of $\zeta _{\rm w} = 10$. Vertical lines ${\rm St}^\ast $ and $\zeta _{\rm w}$distinguish the settling, hyperbolic, and three-body regimes. Analytic fits from our recipe are shown by the dashed-black, solid-gray, and solid black curves, respectively, and denoted $\tilde{b}_{\rm set}, b_{\rm hyp}$ and $b_{\rm 3b}$. The corresponding impact radii obtained from the numerical integrations are shown by crosses ( $\alpha _{\rm p}=10^{-3}$) and circles ( $\alpha _{\rm p}=10^{-5}$).

Open with DEXTER

Figure 8 compares the impact radii obtained from the numerical integrations (symbols) with the analytical prescriptions (curves) for a headwind velocity of $\zeta _{\rm w} = 10$ and for a planet size of $\alpha _{\rm p}=10^{-3}$ and 10-5. We have used Eq. (21) to convert the projected impact parameters $\Delta x_{\rm S}$ to true impact parameters at the interaction point ($b_\sigma $). To do so we used the parabola solution, Eq. (20), to evaluate the gradients $({\rm d}y/{\rm d}x)$ at the starting point ($x_{\rm S}$) and at the interaction point xI. We determine the maximum value of $x_{\rm S}$ that resulted in a collision with the planet at the specified $\alpha_{\rm p}$ and took this value to compute $({\rm d}y/{\rm d}x)_{\rm S}$. To compute $({\rm d}y/{\rm d}x)_I$ we evaluated Eq. (20) at the approach radius $b_{\rm app}$. Here, for $b_{\rm app}$ we took the impact radius obtained from our analytical model described above, except for interactions in the 3-body regime where we always use $b_{\rm app}=2.5$.

At low ${\rm St}$ the interactions take place in the settling regime. Impact radii are rather large, particularly near the ${\rm St}^\ast = 12/\zeta_{\rm w}^3$ transition line, and independent of $\alpha_{\rm p}$ (the cross and circle symbols overlap), implying that the physical impact parameter is larger at larger disk radii. For intermediate Stokes numbers the hyperbolic regime is valid and impact radii are much smaller. However, for ${\rm St}>\zeta_{\rm w}$ impact radii once again increase. The behavior is rather erratic, though, with peaks at ${\rm St}=10$ and 300 and a depression at 102, valid for both $\alpha_{\rm p}$. We found that this complex behavior can be attributed to the trajectories that originate from the third quadrant. Initially, for low Stokes numbers, these are absent due to the strong radial drift. However, at a critical Stokes number the contribution of particles approaching from interior orbits (negative $y_{\rm S}$) becomes important. We do not have a full understanding how these outliers can be modeled analytically.

The analytic fits to the various regimes are given by the dashed curve (for settling), solid gray curve (hyperbolic) and solid black curve (three-body). From Fig. 8 it is obvious that the transition between the settling and the hyperbolic regime is not so sharp. Even particles that have ${\rm St}>{\rm St}^\ast$ display settling behavior. For these reasons, we have extended the validity of the settling regime beyond ${\rm St}^\ast $ by adding an exponential term, i.e.,

\begin{displaymath}\tilde{b}_{\rm set} = b_{\rm set} \exp[-({\rm St}/{\rm St}^\ast)^\gamma],
\end{displaymath} (32)

where $b_{\rm set}$ is the solution of Eq. (27) and $\gamma$ a constant that we empirically fix at $\gamma=0.65$. The impact parameter in the hyperbolic regime, ${\rm St}^\ast<{\rm St}<\zeta_{\rm w}$, is then given by the maximum of $\tilde{b}_{\rm set}$ and $b_{\rm hyp}$. The little depression that can be seen at ${\rm St}=1$ is caused by the fact that the approach velocity Eq. (29) has a maximum here. For larger ${\rm St}$ the approach velocity strongly decreases and the gravitational focusing strongly increases. Nevertheless, impact parameters at ${\rm St}=\zeta_{\rm w}$ are even larger than $b_{\rm hyp}$ and are better fitted by $b_{\rm 3b}$. Initially, gas drag very effectively captures bodies within the Hill sphere and impact radii are rather large. However, the capture probability decreases as $1/{\rm St}$ and for large ${\rm St}$ we retrieve the gas-free limit, Eq. (5).

5.3 Collision rates

\begin{figure}
\par\includegraphics[width=9cm]{14903fg9}
\end{figure} Figure 9:

Contours of $\log P$ according to the analytic prescription. Curves are the same as in Fig. 6. (This figure is available in color in electronic form.)

Open with DEXTER

In Fig. 9 we plot contours of the collision rate, that is, we plot $P=2b_\sigma v_{\rm a}$ as function of ${\rm St}$ and $\zeta _{\rm w}$ using the prescription outlined in Table 2. This figure should be compared with Fig. 6.

The curves in Fig. 9 are much smoother due to the much finer grid that the analytic formulation permits. The transition lines, ${\rm St}={\rm St}^\ast$ and ${\rm St}=\zeta_{\rm w}$ are clearly identified. Our analytic formalism fails to reproduce the $\log P = -0.5$ band towards the upper-right of Fig. 6. However, the overall match is satisfactory; for 90% of the 221 grid points the analytic and numerical results lie within 30% of each other.

5.4 Summary of impact parameter recipe

Table 2:   Summary of the analytic recipe to obtain the impact radii $b_\sigma $ and approach velocities $v_{\rm a}$.

Table 2 provides an executive summary of how the collisional parameters can be obtained using the analytic prescription. First, one converts the physical parameters (headwind velocity, disk radius, friction time, etc.) into the dimensionless quantities $\zeta_{\rm w}, \alpha_{\rm p}$ and ${\rm St}$. The corresponding impact radii for the three regimes are calculated in the second step. Then, in step 3, the appropriate collision regime is determined by comparing the Stokes number with ${\rm St}^\ast $ (Eq. (26)) and $\zeta _{\rm w}$ (Eq. (15)), see also Fig. 7. Dependent on the applicable regime, the final impact parameter is obtained by taking the maximum of two impact radii (step 4). Other quantities ($v_{\rm a}$ and $b_{\rm app}$) also depend on the collision regime.

Then, these results can be converted back to physical units by multiplication of $R_{\rm h}$ and $v_{\rm h} = R_{\rm h} \Omega$ (Eq. (2)) for, respectively, lengths and velocities. The 2D-collision rate is then obtained from Eq. (3). The 3D-collision rate may be estimated by multiplication by a factor $\max(1, H_{\rm p}/b_\sigma)$ (see Sect. 6), where $H_{\rm p}$ is scaleheight of the particles.

6 Significance to the growth of pre-planetary bodies

In the previous sections we have outlined a general approach to analytically derive impact radii and collision rates in Hill coordinates. But what does all of this imply for the growth of preplanetary bodies? Perhaps the best way to illustrate this point is to calculate the accretion timescale

\begin{displaymath}T_{\rm ac}^{\rm 2D} = \frac{M}{{\rm d}M/{\rm d}t} = \frac{4\p...
...o_{\rm s} R_{\rm p} \alpha_{\rm p}^2}{3 \Sigma P} \Omega^{-1},
\end{displaymath} (33)

where $P_{\rm col}$ is the dimensional accretion rate and $P=P_{\rm col}/R_{\rm h} v_{\rm h}$ the dimensionless, the quantity plotted in Figs. 6 and 9. From Eq. (33) we see that the accretion timescale is inversely proportional to P but also scales with $R_{\rm p}$. We further write Eq. (33) in terms of semi-major axis a0 by substitution of Eq. (16) for $\alpha_{\rm p}$ and $\Omega(a)$ for a solar-mass star
                        $\displaystyle T_{\rm ac}^{2D}$ = $\displaystyle \frac{6.7}{P} \left( \frac{\rho_{\rm s}}{3\ {\rm g\ cm^{-3}}} \right)^{-1/3} \left( \frac{\Sigma}{1\ {\rm g\ cm^{-2}}} \right)^{-1}$  
    $\displaystyle \times\left( \frac{R_{\rm p}}{100\ {\rm km}} \right) \left( \frac{a}{{\rm AU}} \right)^{-1/2}\ {\rm yr}.$ (34)

The inverse dependence on disk radius may seem surprising but one has to realize that P via $\alpha_{\rm p}$ and $\zeta _{\rm w}$ implicitly depends on a. Nevertheless, Eq. (34) shows that 2D accretion may be especially advantageous in the outer disks.

The 2D regime, however, may not be applicable to small particles since any breath of turbulence will stir them up. The height of the particle layer may be obtained by equating particle diffusion and settling timescale; i.e.,

\begin{displaymath}\frac{H_{\rm p}}{H_{\rm g}} \approx \min\left(1, \sqrt{\frac{\alpha_{\rm t}}{{\rm St}}} \right),
\end{displaymath} (35)

(Dubrulle et al. 1995; Youdin & Lithwick 2007; Carballido et al. 2006) where $\alpha_{\rm t}$ is the Shakura & Sunyaev (1973) viscosity parameter for turbulent diffusion. In a 3D setting, the particle scaleheight can exceed the impact parameter $b_\sigma $; then, only a fraction, $b_\sigma/H_{\rm p}$, of the particles take part in the interaction and the accretion timescale is correspondingly longer,

\begin{displaymath}T_{\rm ac}^{\rm 3D} \approx T_{\rm ac}^{\rm 2D} \times \max\left(1, \frac{H_{\rm p}}{b_\sigma}\right)
\end{displaymath} (36)

(with $b_\sigma $ in physical units). The 3D correction factor significantly increases collision timescales for small particles ($H_{\rm p}$ is large) and the hyperbolic regime ($b_\sigma $ is small).

\begin{figure}
\par\mbox{\includegraphics[width=7.5cm]{14903fg10a.eps} \includegraphics[width=7.5cm]{14903fg10b.eps} }
\end{figure} Figure 10:

The 3D growth timescale $T_{\rm ac}^{\rm 3D}$ as function of Stokes number (or particle size) and size of the protoplanet. Contour lines of $\log_{10} T_{\rm ac}^{\rm 3D}/{\rm yr}$are shown. All the solid density is assumed to be contained in particles of the indicated size. (A) At 5 AU for a surface density of $\Sigma = 2~{\rm g~cm^{-2}}$. (B) At 35 AU for a surface density of $\Sigma=0.1~{\rm g~cm^{-2}}$. In (A) the wider spacing between the tickmarks of the upper axis indicates particles enter the Stokes drag regime. (This figure is available in color in electronic form.)

Open with DEXTER

In Fig. 10a we have plotted contours of the 3D growth timescales for a disk radius a0=5.2 AU, $\rho_{\rm s}=3\ {\rm g\ cm^{-3}}$ (making $\alpha _{\rm p}=10^{-3}$), $\Sigma = 2\ {\rm g\ cm^{-2}}$, $H_{\rm g} = 0.25\ {\rm AU}$, $\Delta v_{\rm hw} = 30\ {\rm m\ s^{-1}}$, and $\alpha_{\rm t} = 10^{-4}$. The calculated accretion timescales assume that all the solid density is in particles of a single size. Due the inclusion of the $H_{\rm p}/b_\sigma$ factor the structure is quite different from that of Fig. 9. However, the contrast between the hyperbolic and settling regimes is still clearly visible and has in fact even increased due to the correction factor for the vertical structure. Note that for the bigger bodies, which settle into a thin plane, Fig. 10 still assumes that their eccentricities and inclinations are absent (low velocity regime).

Figure 10 tells a few interesting points. First, it can be clearly seen that growth of $\sim$km-size planetesimals by accretion of small particles ( ${\rm St}<1$) takes a (perhaps prohibitively) long time. Two mechanisms conspire. First, the small particles couple effectively to the gas which dilutes their number densities near the midplane where the planetesimals are residing. Of course, this statement depends on the strength of the turbulence that prevents the particles from settling effectively; timescales will be shorter for lower turbulent strength parameter, $\alpha_{\rm t}$. However, even in a completely laminar disk we may expect shear turbulence to develop (Weidenschilling 1980), which strength may be equivalent to $\alpha_{\rm t}$-values of $\sim$10-6 (Johansen et al. 2006; Cuzzi & Weidenschilling 2006). Second, small particles, being strongly coupled, move with the gas, at a relative velocity of $\Delta v \approx v_{\rm hw} \gg v_{\rm esc}$, where $v_{\rm hw}$ is the velocity of the headwind and $v_{\rm esc}$ the escape velocity of the planetesimal. Therefore, small particles lack gravitational focusing and it is hard to avoid the conclusion that sweepup of small particles by $\sim$km size planetesimals is slow. In order to grow, planetesimals have to accrete among themselves.

However, the situation completely reverses when protoplanets sizes reach $\sim$103 km: for these bodies, accretion of cm to m-size particles becomes very rapid: in only $\sim$103 yr the protoplanet can double in size. This is entirely due to the increased cross section in the settling regime. For ``optimal'' parameters ( ${\rm St}\sim1$, $\zeta_{\rm w}\la1$) the combined effect of gravitational focusing and gas damping results in impact parameters of $b_\sigma \sim 0.5R_{\rm h}$ - larger than what hitherto has been thought possible ( $b_{\rm gf}$, see Eq. (5)). Accreting at impact parameters of the order of the Hill sphere is fast in any case but since $R_{\rm h}$ increases with disk radii it is especially impressive for the outer disk, see Fig. 10b.

For ${\rm St}\sim1$particles, accretion is fast - even though it is inefficient due to the strong radial drift. We denote the probability that radially-inward drifting particles become accreted by the protoplanet $P_{\rm eff}$. Since the drift flow is $2\pi a v_{\rm r} \Sigma$ ($\approx$ $0.15~M_\oplus~{\rm yr}^{-1}$ for ${\rm St}=1$ particles at 5 AU) $P_{\rm eff}$ is given as

\begin{displaymath}P_{\rm eff} = \frac{P_{\rm col}}{2\pi a v_{\rm r}} \times \max\left(1, \frac{H_{\rm p}}{b_\sigma}\right)\cdot
\end{displaymath} (37)

(where we have again included the scaleheight correction factor). For a protoplanet of $R_{\rm p}=10^3~{\rm km}$ we find that ${\rm St}\sim1$ particles are accreted at an efficiency of only $P_{\rm eff}\approx 10^{-4}$. The efficiency increases away from the ${\rm St}=1$ line ($v_{\rm r}$ is lower) and towards larger protoplanet size and larger disk radii (larger $R_{\rm h}$). Levison et al. (2010), using N-body techniques, also finds that accretion of small fragments was very inefficient (and concluded that it was therefore unlikely) except for a few specific particle sizes, that may have coincided with the peaks in Fig. 10. To increase the accretion efficiency, smaller particles have to be accreted. Dust fragmentation (Birnstiel et al. 2009) or radial diffusion (Ciesla 2009) may be promising mechanisms to retain mm-size particles in the (outer) disk, where they are observed on $\sim$Myr timescales (e.g., Ricci et al. 2010; Lommen et al. 2009).

It is instructive to compare the accretion timescales of Fig. 10a to detailed hydrodynamical simulations involving ${\rm St}\sim1$ particles (Johansen et al. 2007; Johansen & Lacerda 2010). In Johansen et al. (2007) a dense particle layer of ${\rm St}\sim1$ boulders collapses into a Ceres-mass planet ($R\sim500$ km), that rapidly accretes the remaining boulders on timescales of perhaps 10 yr. Although from Fig. 10a a Ceres-mass protoplanet in combinations with ${\rm St}\sim1$ particles form the optimal growth conditions, our accretion timescale of 103 yr is still two orders of magnitude higher than what can be inferred from Johansen et al. (2007). However, a direct comparison is perhaps not so meaningful since in the Johansen et al. (2007) simulations the ${\rm St}\sim1$ particles are highly clumped and exert a strong feedback effect on the gas (Johansen & Lacerda 2010 discuss some alternate settings). Feedback effects are not taken into account in this study.

7 Discussion

We summarize the key assumptions that have been employed in this study:

  • a drag law linear in velocity;
  • neglect of resonance trapping of particles;
  • a smooth, laminar disk (only drift motions) without local pressure fluctuations;
  • neglect of the gas flow around the protoplanet and of a possible atmosphere surrounding the proto(planet);
  • a dynamically cold protoplanet on a circular, non-migrating, orbit.
The assumption of a linear drag law implies that this study - and in particular the analytical prescriptions that have been derived - apply for particles smaller than $\sim$ $s_{\rm max}$ (see Eq. (13)) only. But this still covers an appreciable size range, especially for the outer disk.

The adopted flow pattern in our study is unrealistic since it does not take account of the presence of the planet. Of course, streamlines will have to bend around the object and this will affect the motion of the particle. Tiny dust grains can only collide with a dust aggregate when the aggregate size is less than the mean free path of the gas molecules (Wurm et al. 2001; Sekiya & Takeda 2005). However, this restriction probably applies only for small particles. If we assume that the flowlines change over the lengthscale of the protoplanet, it follows that particles of $t_{\rm s} < R_{\rm p}/v_{\rm hw}$ are too tightly coupled to the gas to become accreted. For example, for $R_{\rm p} = 100$ km only ${\rm St}<10^{-4}$ particles are affected and our model overestimates the (already low) accretion rates. More serious, perhaps, is our neglect of collective effects due to strong particle volume densities, as this will provide a feedback effect on the gas, affecting both the flow pattern as well as the drift rates. This will probably be important for ${\rm St}\sim1$ particles and could significantly enhance the accretion rate (see our discussion at the end of Sect. 6).

For simplicity, our analysis only included drift motions. We have, for example, neglected a systematic accretion (or decretion) flow of the gas. In the turbulent $\alpha$-model this gas moves in at a velocity $\sim$ $\alpha_{\rm t} c_{\rm g} H_{\rm g}/a$. Equating this expression with the radial drift velocity (Eq. (19)) we find that for particles ${\rm St}<\alpha_{\rm t}$ the systematic accretion flow will dominate. This will affect the expressions for the accretion rates. Likewise, turbulent motions can become more effective to move particles around than drift motions, which affects the input parameter $v_{\rm a}$ in our model described in Sect. 5. In the $\alpha_{\rm t}$-model large eddies transport particles at velocities $\sim$ $\alpha_{\rm t}^{1/2}c_{\rm g}$ (Cuzzi & Weidenschilling 2006) and turbulent velocities will dominate the drift motions for $\alpha_{\rm t}>\eta$ (see Eq. (10)). The presented model may still be valid though, if the turbulent motions are included in the definition of the approach velocity, $v_{\rm a}$.

Mean motion resonances may halt the particle long before it drifts to the Hill sphere (Weidenschilling & Davis 1985). In such a situation the inward-directed drag force is balanced by the outward resonant perturbations. Weidenschilling & Davis (1985) showed that the strength of the perturbations is proportional to the planet's mass and to the resonance number j. Thus, smaller particles (which experience a stronger drag force) move into a higher resonance. However, this trend will not pursue indefinitely as at some maximum j resonances will overlap and the effect is lost. Paardekooper (2007) simulated particle accretion onto gas giants and showed that $\sim$m-size particles ( ${\rm St}\sim1$) avoid resonance trapping for Jupiter-mass ($M_{\rm J}$) planet. For a planet of $0.1\ M_{\rm J}$ the critical Stokes number has risen to ${\rm St}=10$ and for lower mass planets it will even be larger. Therefore, our results are not so much affected by resonance trapping when the protoplanet and core-formation stages are considered.

Our assumptions of a completely ``inert'' protoplanet is also peculiar. In the oligarchic growth regime (where most of the mass resides in leftover planetesimals) dynamical friction will keep the eccentricities of the most massive bodies small (e.g., Kokubo & Ida 2000); however, if most of the mass is transferred to the protoplanets their motion will become eccentric (see Kary & Lissauer 1995 for accretion probabilities of protoplanets on eccentric orbits). Likewise, type-I radial migration (Tanaka et al. 2002) is not incorporated in our framework. These effects will again become important for already evolved protoplanets of masses > $0.1~M_\oplus$.

The feedback effect of the protoplanet on the structure of the gas disk is also neglected. The protoplanet's gravity influences the gas disk at larger distances, which could invalidate our approximation of a smooth (global) pressure structure. Indeed, particles have a tendency to drift to high pressures regions and the ${\rm St}\sim1$ particles may be most affected by this process, piling up at a pressure bump instead of proceeding to the protoplanet. Paardekooper (2007) finds that this effect (together with the resonant trapping of bigger bodies) shuts off all accretion of particles sizes above $\sim$ $10{-}100~\mu$m! However, this particle trapping is applicable for evolved planets only. Muto & Inutsuka (2009) found that the criterion for particle trapping

\begin{displaymath}M_{\rm p} \ga \eta \frac{H_{\rm g}}{a} M_\star \sim 10\ M_\oplus
\end{displaymath} (38)

for solar mass stars (Equation (38) is apparently independent of particle size or Stokes number).

Long before this size is reached, protoplanets bind the nebular gas and form atmospheres that will enhance the capture radius (Tanigawa & Ohtsuki 2010; Inaba & Ikoma 2003). According to the results of Inaba & Ikoma (2003) this will perhaps become important when oligarchs reach $0.1~M_\oplus$[*]. Our expression for the impact radius and collision rates, therefore, are lower limits when protoplanets are surrounded by a thick atmosphere.

In summary, most of the mentioned effects become relevant for evolved (gas) planets only. Most damaging to our analysis are the pressure fluctuations that could virtually shut off accretion, or make it very difficult to model it analytically. However, this effect may only become effective for large planet masses (Eq. (38)). For lower mass planets, the build-up of a dense atmospheres will enhance the accretion rates with respect to our prescription. We believe that for dynamically-cold protoplanets below $0.1~M_\oplus$ our prescription should be quantitatively correct if collective effects can be neglected. In future work we intend to test the validity of our analytic expressions in more convoluted environments that incorporate some of the above processes.

8 Summary

We have developed a framework for the calculations of particle-protoplanet interactions in a gaseous environment. This involves the integration of the equations of motions in the circularly restricted three body problem including drag forces. Using the above mentioned simplifications - most notably the assumption of a linear drag force, a smooth background density and headwind velocity, and a 2D setting - we were able to reduce the general problem to a state that includes only two dimensionless parameters: the dimensionless headwind velocity $\zeta _{\rm w}$ and the dimensionless stopping time (Stokes number, ${\rm St}$), see Sect. 2. A large parameter study of particle trajectories has been conducted, from which, as function of $\zeta _{\rm w}$ and ${\rm St}$, the impact radii are derived. We find that three accretion modes can be distinguished:

  • Settling encounters. Particles settle to the protoplanet and the impact radius is independent of the size of the latter;
  • Hyperbolic encounters. Accretion proceeds like the usual gravitational focusing with the approach velocity being influenced by gas drag.
  • (Drag enhanced) three body encounters. Interactions take place on the scale of the Hill radius and gas drag causes a fraction of the particles to become captured, which settle to the protoplanet.
We have developed analytic recipes for all three encounters and found them to match very well to the results of our numerical study (except perhaps for the three-body regime). The recipe is summarized in Table 2 and Sect. 5.4. In Sect. 6 we have extended our approach to the usual 3D-setting and calculated the accretion times of (proto)planets by sweepup of particles. We found that small particles are very unlikely candidates to grow to planetesimals of 10-102 km, since their trajectories are tightly coupled to the gas. However, if the protoplanet has reached a size of $\sim$103 km it can very quickly accrete ${\rm St}\sim1$ particles through the settling mechanism. Since these fragments do not suffer from the protoplanet gravitational scattering (they quickly circularize), accretion under such conditions represents an avenue for quick growth, especially in the outer disk - provided they are not lost by radial drift.

Acknowledgements
The authors appreciate the helpful comments of the referee, Stuart Weidenschilling. C.W.O. acknowledges a grant from the Alexander von Humboldt Foundation and the hospitality of the Max-Planck-Institute for Astronomy for hosting him. C.W.O. also appreciates stimulating discussions with Tilman Birnstiel, Kees Dullemond, Christoph Mordasini, and Marco Spaans.

Appendix A: Asymptotic limits of Eq. (22)

In this appendix we consider the asymptotic limits of Eq. (22) and show the correspondence to the findings of Kary et al. (1993).

Three limits of Eq. (22) can be identified

\begin{displaymath}P_{\rm geo} =\left\{ \begin{array}{ll}
2\alpha_{\rm p} \zeta...
..._{\rm p} {\rm St}/\zeta_{\rm w} \gg 1. \\
\end{array}\right.
\end{displaymath} (A.1)

These three regimes correspond to the cases where the square-root term of Eq. (22) evaluates to $1/2{\rm St}$, 1, and $3\alpha_{\rm p}{\rm St}^2/4\zeta_{\rm w}$, respectively. In all limits we have assumed that $\alpha_{\rm p}\ll\zeta_{\rm w}$.

Rewritten in physical units Eq. (A.1) reads (i.e., we multiply by $R_{\rm h}v_{\rm h}$)

\begin{displaymath}P_{\rm col}^{\rm geo} = \left\{ \begin{array}{ll}
2R_{\rm p}...
...\ R_{\rm p} {\rm St}/v_{\rm hw} \gg 1. \\
\end{array}\right.
\end{displaymath} (A.2)

The interpretation of the first two limits is straightforward. If ${\rm St}\ll1$ particles arrive with the headwind velocity, $v_{\rm hw}$. In the 2D-setting the cross section is $2R_{\rm p}$, so the collision rate is $2R_{\rm p} v_{\rm hw}$. Similarly, in the second limit the particles approach from the x-direction at a speed of $2v_{\rm hw}/{\rm St}$ (see Eq. (19)).

In the third limit the particles approach once again from the y-direction. However, for these very large particles, the approach velocity is given by the Keplerian shear ( $3\Omega x/2$) instead of the headwind, see Eq. (19). The approach velocity at the point of intersection, i.e., at $x=R_{\rm p}$, is then $v_{\rm a} = 3R_{\rm p}\Omega/2$. Multiplied by $2R_{\rm p}$ this reduces to the given expression. The dependence on $R_{\rm p}$ may seem counter intuitive but is natural in situations that involve shear.

The study of Kary et al. (1993) concerned massive particles (i.e., the third limit of Eq. (A.2)). Kary et al. (1993) gave an analytic expression for the impact probability or efficiency $P_{\rm eff}$ of a particles while crossing the semi-major axis of the protoplanet due to radial drift[*]:

\begin{displaymath}P_{\rm eff} = \frac{3R_{\rm p}^2 \Omega^2}{16\pi K a v_{\rm hw}},
\end{displaymath} (A.3)

where K is the drag constant for a drag law that quadratically depends on velocity.

In our case $P_{\rm eff}$ can be found by taking the ratio of the collision rate $P_{\rm col}$ to the mass inflow rate $\vert 2\pi a v_x\vert$ of the particles. For the third limit of Eq. (A.2) we obtain

\begin{displaymath}P_{\rm eff}^{\rm geo} = \frac{P_{\rm col}^{\rm geo}}{2\pi a v_x} = \frac{3 {\rm St} \Omega R_{\rm p}^2 }{4\pi a v_{\rm hw}},
\end{displaymath} (A.4)

where we used that $v_x = 2 v_{\rm hw}/{\rm St}$ in this regime. Equation (A.4) is different from Kary et al. (1993)'s result for three reasons:
1.
Kary et al. (1993) considers a drag force quadratic in velocity where particles move at a drift velocity $v_x = -2K v_{\rm hw}^2/\Omega$ (K has units ${\rm cm}^{-1}$). However, we can mimic the linear drift law by substitution of $K=\Omega/v_{\rm hw} {\rm St}$ (cf. Eq. (19)) into Eq. (A.3).
2.
In our approach we have not accounted for the variation of the approach velocity over the impact range. More correctly, the mean approach velocity is $\overline{v_{\rm a}}=3\Omega R_{\rm p}/4$.
3.
Kary et al. (1993) do not take account of the planetesimals coming from the negative y-direction (the third quadrant), as they (correctly) argue that any such body should already have impacted during its approach from the first quadrant. This is due to the highly symmetrical setting in this limit, which our naive reasoning above does not account for. Equation (A.4) is therefore too large by a factor of two.
With these corrections Eqs. (A.3) and (A.4) agree.

Appendix B: The settling path

In the limit of ${\rm St}\ll1$ we can use the approximation that the particle is always in the settling regime, i.e., its velocity is given by $\vec{v} = \vec{F}_{\rm g} t_{\rm s}$ (in Hill units):

\begin{displaymath}v_x = -\frac{3x}{r^3} {\rm St},
\end{displaymath} (B.1a)

\begin{displaymath}v_y = -\frac{3y}{r^3} {\rm St} -\zeta_{\rm w}.
\end{displaymath} (B.1b)

Thus, the particle path obeys the differential equation

\begin{displaymath}\frac{{\rm d}y}{{\rm d}x} = \frac{v_y}{v_x} = \frac{y}{x} + \frac{\zeta_{\rm w}}{3{\rm St}} \frac{(x^2+y^2)^{3/2}}{x},
\end{displaymath} (B.2)

which is slightly simplified if expressed in angular coordinates by the substitution $y=x \tan \theta$. Then, ${\rm d}y/{\rm d}x = \tan \theta + x(1 + \tan^2 \theta) {\rm d}\theta/{\rm d}x$ and Eq. (B.2) becomes in terms of x and $\theta $

\begin{displaymath}x (1+\tan^2 \theta) \frac{{\rm d}\theta}{{\rm d}x} = \frac{\zeta_{\rm w}}{3{\rm St}} (1+\tan^2 \theta)^{3/2},
\end{displaymath} (B.3)

which is equivalent to

\begin{displaymath}\cos \theta \frac{{\rm d}\theta}{{\rm d}x} = \frac{\zeta_{\rm w}}{3{\rm St}} x.
\end{displaymath} (B.4)

Straightforward integration gives the solution

\begin{displaymath}\sin \theta = \frac{\zeta_{\rm w}}{6{\rm St}} x^2 +C,
\end{displaymath} (B.5)

with C the integration constant, which we obtain by the requirement that at $\theta = \pi/2$ the particle starts out at x=x0. Thus, the full solution to the orbit under these conditions is

\begin{displaymath}1-\sin \theta = \frac{\zeta_{\rm w}}{6{\rm St}}(x_0^2 -x^2).
\end{displaymath} (B.6)

During the encounter $\theta $ decrease from $\pi/2$ to $-\pi/2$ and $1-\sin \theta$ increases from 0 to 2. The particle's x coordinate then decreases by an amount that depends on the numerical value of $\zeta_{\rm w}/6{\rm St}$ and x0. The particle settles to the origin if x=0 can be reached. The largest value of x0 for which this is possible is $b = x_0 = \sqrt{12 {\rm St}/\zeta_{\rm w}}$. This is the impact parameter $b_\sigma $.

For example, for ${\rm St}=10^{-2}$ and $\zeta _{\rm w}=1$, we obtain $b_\sigma=0.35$, a value that is reasonably close to the numerically derived $x_{\rm S} = 0.38$ (see Fig. 5). For lower Stokes number the agreement becomes better.

In Sect. 5.1.2 we have used the expression $\Delta v = 3{\rm St}/b^2$ for the impulse change in the settling regime. For $b=\sqrt{12 {\rm St}/\zeta_{\rm w}}$ this corresponds to a velocity change of $\Delta v = \zeta_{\rm w}/4$. We further argued that for settling encounter the approach velocity $v_{\rm a}$ should be changed by an amount $\Delta v \sim v_{\rm a}$. For small Stokes values $v_{\rm a} = \zeta_{\rm w}$. Therefore, for ${\rm St}\ll1$ the criterion for accretion by settling becomes $v_{\rm a} \ge \zeta_{\rm w}/4$. In Sect. 5.1.2 we have applied this criterion generally.

References

Footnotes

... ${\rm St}=t_{\rm s} \Omega_0$[*]
Note that the Stokes number in this study simply indicates the dimensionless friction time; it is not necessarily the same as the Stokes number used in turbulent studies, ${\rm St}=t_{\rm s}/t_{\rm L}$ where $t_{\rm L}$ is the turn-over timescale of the largest eddies. For $t_{\rm L} = \Omega^{-1}$ the definitions agree (Youdin & Lithwick 2007).
...[*]
 In this and the next two sections lengths ( tex2html_wrap_inline3599, etc.) and velocities (v) are expressed in dimensionless (Hill) units, unless otherwise specified.
...[*]
In this and the next two sections lengths ( $x,y,b_\sigma $, etc.) and velocities (v) are expressed in dimensionless (Hill) units, unless otherwise specified.
...[*]
In this and the next two sections lengths ( $x,y,b_\sigma $, etc.) and velocities (v) are expressed in dimensionless (Hill) units, unless otherwise specified.
...$0.1~M_\oplus$[*]
We remark that the atmosphere calculations of Inaba & Ikoma (2003) do not take into account headwind flow, perhaps important for low $M_{\rm p}$, which would destroy the spherical symmetry of the problem.
... drift[*]
Note that we give the dimensional form. Equation (7) of Kary et al. (1993) is expressed in dimensionless units (but not in Hill units).
Copyright ESO 2010

All Tables

Table 1:   Dimensional and dimensionless parameters.

Table 2:   Summary of the analytic recipe to obtain the impact radii $b_\sigma $ and approach velocities $v_{\rm a}$.

All Figures

  \begin{figure}
\par\includegraphics[width=8cm]{14903fg1}
\end{figure} Figure 1:

Sketch of particle trajectories in the comoving frame. We consider the motion of the third (test) particle m in the comoving frame of the second body ($M_{\rm p}$, the planet) while including the gravity of the central star. In the gas-free limit zero-eccentricity particles can enter the Hill sphere from both the first and the third quadrant (black curves) but only from specific impact parameters indicated by the hatched regions. Particles arriving at closer impact parameters move on horseshoe orbits. The magnitude and direction of the gas velocity $\vec{v}_{\rm gas}$ as seen from the comoving frame is indicated by the dashed arrows. Particle trajectories including gas drag (solid gray arrows) can be anything depending on the properties of the particle and the gas.

Open with DEXTER
In the text

  \begin{figure}
\par\includegraphics[width=8.5cm]{14903fg2}
\end{figure} Figure 2:

Relation between dimensionless and physical quantities. The Stokes number (${\rm St}$, solid lines) and the dimensionless headwind velocity $\zeta _{\rm w}$ (dashed lines) are plotted on the y-axis as function of the radius s of the particle and the radius $R_{\rm p}$ of the protoplanet. Note the different units of s and $R_{\rm p}$ on the x-axis. The black dot denotes $s_{\rm max}$ (Eq. (13)). Lines are shown for: (i) $c_{\rm g}=10^5~{\rm cm\ s^{-1}}$ and $\rho_{\rm g}=10^{-9}~{\rm g~cm^{-3}}$ at a disk radius of 1 AU (black lines) and (ii) $c_{\rm g} = 6\times10^4~{\rm cm~s^{-1}}$ and $\rho_{\rm g} = 10^{-11}~{\rm g~cm^{-3}}$ at a position of 10 AU (gray lines). The internal density of solids is fixed at $\rho_{\rm s} = 3~{\rm g~cm^{-3}}$, the nebula headwind is $v_{\rm hw}=30~{\rm m~s^{-1}}$, and the mass of the central object is solar.

Open with DEXTER
In the text

  \begin{figure}
\par\includegraphics[width=8.5cm]{14903fg3}
\end{figure} Figure 3:

Without the two body force, particle trajectories as witnessed from the comoving frame obey parabolas. Two trajectories are shown: one that passes through the origin (y0(x)) and one that just hits the target. The corresponding impact parameter b is denoted by arrows. For curved trajectories b is not conserved due to the changing slope of the curves, here indicated by the angle $\theta $. At S we have that $\Delta x_{\rm S} = b_{\rm S}/\!\sin \theta_{\rm S}$ and $\Delta y_{\rm S} = b_{\rm S}/\!\cos \theta_{\rm S}$. $\Delta y_{\rm S}$ is a conserved quantity.

Open with DEXTER
In the text

  \begin{figure}
\par\includegraphics[width=9cm]{14903fg4}
\end{figure} Figure 4:

The minimum distance in units of Hill radii to the origin (center of the protoplanet), $r_{\rm min}$, as function of the x-coordinate of the starting point, $x_{\rm S}$ ($y_{\rm S}$ is fixed at 40 Hill radii). Plotted is $r_{\rm min}$ for a gas-free system (left) and a system that is characterized by the parameters ${\rm St}=10$ and $\zeta _{\rm w}=1$ (right). Several bands are labeled. The inclusion of gas drag shifts the bands to larger $x_{\rm S}$ while merging several chaotic features within the chaotic c-band.

Open with DEXTER
In the text

  \begin{figure}
\par\includegraphics[width=17cm]{14903fg5}
\end{figure} Figure 5:

Examples of planet-particle interactions for different values of the dimensionless headwind velocity $\zeta _{\rm w}$ and coupling parameter ${\rm St}$. For typical nebula parameters particles of ${\rm St}=10$ correspond to loosely coupled m-size particles, whereas ${\rm St}=0.01$ are more strongly coupled cm-size particles, see Fig. 2. Likewise, $\zeta _{\rm w}=1$ corresponds to protoplanets of $R_{\rm p}\sim 10^3$ km in radius, while $\zeta _{\rm w}=100$ corresponds to $R_{\rm p}\sim 10$ km planetesimals. (A) Two particles of ${\rm St}=10$ experience a close encounter within the Hill sphere (dotted circle). The $x_{\rm S}=3.9$ particle is captures and settles to the planet, whereas the other particle is ejected from the Hill sphere (The Keplerian shear eventually causes it to resurface at the other side of the Hill sphere). (B) Strong gas coupling, ${\rm St}=0.01$. There is a competition between the gravitational pull of the planet and the drag force directed towards negative y. (C) Close encounters at large $\zeta _{\rm w}$ without settling (see inset). (D, E) Examples of particle trajectories originating from interior orbits. (F) Radially approaching orbits. (This figure is available in color in electronic form.)

Open with DEXTER
In the text

  \begin{figure}
\par\includegraphics[width=9cm]{14903fg6}
\end{figure} Figure 6:

Contour plot of collision rates obtained from the numerical integrations for $\alpha _{\rm p}=10^{-3}$. Contours of $\log_{10} P$ are shown as function of the Stokes number and the headwind velocity. Contour levels are indicated except for $\log_{10} P = 0$ (thick solid), -0.25 (dotted) and -0.75 (dashed). The accretion rate in the gas-free regime is $\log_{10} P_{\rm gf} \approx -0.46$. (This figure is available in color in electronic form.)

Open with DEXTER
In the text

  \begin{figure}
\par\includegraphics[width=9cm]{14903fg7}
\end{figure} Figure 7:

Illustration of the three accretion regimes. In the hyperbolic regime interactions are 2-body encounters and the standard formula for gravitational focusing applies. In the settling regime, accretion proceeds through settling, which enhances the impact parameter $b_\sigma $. Else, for ${\rm St}>{\rm max}(\zeta _{\rm w}, 1)$ the solar gravity affects the encounter and $b_\sigma $ is increased with respect to the gas-free limit ( $b_{\rm gf}$, Eq. (5)) due to particle capture within the Hill sphere. The gray band approximately indicates the zone where the settling solutions permeate into the hyperbolic regime. The letters a-f correspond to the parameters in the panels of Fig. 5.

Open with DEXTER
In the text

  \begin{figure}
\par\includegraphics[width=9cm]{14903fg8}
\end{figure} Figure 8:

Impact radii from the numerical integrations (symbols) and analytic fits (curves) for a headwind velocity of $\zeta _{\rm w} = 10$. Vertical lines ${\rm St}^\ast $ and $\zeta _{\rm w}$distinguish the settling, hyperbolic, and three-body regimes. Analytic fits from our recipe are shown by the dashed-black, solid-gray, and solid black curves, respectively, and denoted $\tilde{b}_{\rm set}, b_{\rm hyp}$ and $b_{\rm 3b}$. The corresponding impact radii obtained from the numerical integrations are shown by crosses ( $\alpha _{\rm p}=10^{-3}$) and circles ( $\alpha _{\rm p}=10^{-5}$).

Open with DEXTER
In the text

  \begin{figure}
\par\includegraphics[width=9cm]{14903fg9}
\end{figure} Figure 9:

Contours of $\log P$ according to the analytic prescription. Curves are the same as in Fig. 6. (This figure is available in color in electronic form.)

Open with DEXTER
In the text

  \begin{figure}
\par\mbox{\includegraphics[width=7.5cm]{14903fg10a.eps} \includegraphics[width=7.5cm]{14903fg10b.eps} }
\end{figure} Figure 10:

The 3D growth timescale $T_{\rm ac}^{\rm 3D}$ as function of Stokes number (or particle size) and size of the protoplanet. Contour lines of $\log_{10} T_{\rm ac}^{\rm 3D}/{\rm yr}$are shown. All the solid density is assumed to be contained in particles of the indicated size. (A) At 5 AU for a surface density of $\Sigma = 2~{\rm g~cm^{-2}}$. (B) At 35 AU for a surface density of $\Sigma=0.1~{\rm g~cm^{-2}}$. In (A) the wider spacing between the tickmarks of the upper axis indicates particles enter the Stokes drag regime. (This figure is available in color in electronic form.)

Open with DEXTER
In the text


Copyright ESO 2010

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.