CAN PLANETESIMALS FORM BY COLLISIONAL FUSION?

Planet formation is hypothesized to occur through the growth of protoplanetary bodies formed from gas, dust, and ice grains in an accretion disk around a central star (Armitage 2010). The complex scenario of the planet formation process involves the following four stages. First, the initial collapse of interstellar gas to create the central protostar (∼0.1 Myr); second, the slow accretion of mass onto the star and the formation of primary planetesimals within the evolving accretion disk (∼Myr); third, a phase (∼Myr) of reduced accretion rate allowing the photoevaporative wind to divide the disk into an inner and an outer region at a radius determined by the ratio of the stellar accretion rate to the mass loss rate due to photoevaporation; and finally, there is a clearing phase (∼0.1 Myr) during which the inner disk accretes onto the star while the lightest elements of the outer disk are removed due to direct exposure to photoevaporative UV flux. Recent cosmochemical evidence reveals that the long held view of a ∼Myr age difference between Ca-Al-rich inclusions (CAIs) within carbonaceous chondrite meteorites and chondrules within chondrites can be refuted (Connelly et al. 2012). To the extent that these data demonstrate commensurability over disk lifetime scales of CAI and chondrule formation, the detailed transient development of matter within circumstellar disks becomes all the more compelling for studies that can isolate essential physical processes. Here we focus on fundamental aspects of the second stage above. This stage is crucial for understanding how the material that forms the building blocks of planets can organize into bodies that thwart the radiative pressure effects in the subsequent stages that sweep the disk of small particles and gas.⁴

The accretion disk is treated as a two-phase system defined by a fluid phase ("gas") and solid particles ("dust") advected by the fluid. Ubiquitous attractive long range van der Waals and electrostatic interactions facilitate the agglomeration and growth of small (micron or smaller) dust grains that are brought into proximity by the turbulent flow of the gas. However, depending on the material and the mechanical and thermodynamic conditions of a particle–particle collision, sticking (through a number of mechanisms), fragmentation, or bouncing will determine the fate and the size distribution of accreting matter from the small scales upward (Blum & Wurm 2008; Wettlaufer 2010; Zsom et al. 2011).

Because the central star creates a radially decaying pressure gradient, the gas moves at a slightly sub-Keplerian speed. Thus, depending on the position in the disk, there are a range of particle sizes that experience a strong "headwind" and so lose angular momentum, thereby driving them into the central star on timescales as rapidly as a century (Armitage 2010; Youdin 2010). We are concerned with the long-standing problem of how, when objects grow and begin to experience the local headwind, they can accumulate matter sufficiently quickly to slow their drift inward. To focus the question, we examine in some detail how a meter-sized object grows by accretion of small particles mediated by turbulent flows of the gas.

The typical value of the "disk Mach number" $\mathcal {M}_d$ is based on the Keplerian velocity v_kepler, which in the thin disk approximation is

$$\begin{equation} \mathcal {M}_d = \frac{v_{\rm kepler}}{c_{\rm s}} \approx \frac{r}{h}, \end{equation} \tag{ 1 }$$

where r is the radial position in the disk and h is its vertical scale height. At 1 AU, h/r ≈ 0.02 and hence $\mathcal {M}_d \approx 50$ (see, e.g., Armitage 2010, p. 40). Now, as noted above, because the central star creates a radially decreasing pressure gradient, relative to a meter-sized object the gas moves at a sub-Keplerian speed v_wind = ηv_kepler where η can be as small as 10⁻³ depending on the position in the disk.

To understand the effects of the interaction between the dust and the gas, we begin by considering a solid body of spherical shape with radius R_SB, moving through a gas with kinematic viscosity ν and speed v_wind. We estimate its Reynolds number as

$$\begin{equation} \mbox{Re}_{ \rm SB}=\frac{v_{\rm wind}R_{\rm SB}}{\nu } = \frac{v_{\rm wind}}{c_{\rm s}} \, \frac{R_{\rm SB}}{\lambda } \, \frac{\lambda c_{\rm s}}{\nu } \approx \mathcal {M}\, \frac{R_{\rm SB}}{\lambda }, \end{equation} \tag{ 2 }$$

where $\mathcal {M}\equiv v_{\rm wind}/c_{\rm s}$ is the Mach number of the headwind and λ is the mean-free-path of the gas molecules. Importantly, for this estimate we have used the well known expression for the viscosity of gases ν ∼ c_sλ (see, e.g., Lifshitz & Pitaevskii 1981, Section 8). Now, because $\mathcal {M}= \eta \mathcal {M}_d$, we can have $\mathcal {M}\approx 0.05$, and hence, so long as $R_{\rm SB}< \mathcal {M}^{-1}\lambda \approx 20\, \lambda$, the local Reynolds number of the solid body is less than unity. For R_SB ∼ λ the size of the solid body is well below the smallest hydrodynamic length scale in the gas and its motion is then described by the simple drag law:

$$\begin{equation} \frac{d {\boldsymbol v}_{\rm SB}}{dt} = \frac{1}{\tau _{\rm SB}}\left({\boldsymbol v}_{\rm SB}- {\boldsymbol U}\right), \end{equation} \tag{ 3 }$$

where v_SB is the velocity of the particle, U is the local velocity of the gas, and τ_SB is the so-called stopping time describing the deceleration of particle motion relative to the gas. When a particle is smaller than the typical hydrodynamic length scale in the problem, τ_SB is given by the Epstein drag law,

$$\begin{equation} \tau _{\rm SB}^{\rm Ep} = \frac{\rho _{\rm SB}}{\rho _{\rm g}}\frac{R_{\rm SB}}{c_{\rm s}}, \end{equation} \tag{ 4 }$$

where ρ_SB is the material density of the solid particles and ρ_g is the gas density. When $\mathcal {M}^{-1}\lambda > R_{\rm SB}> \lambda$, the relevant drag law is that of Stokes and τ_SB is given by

$$\begin{equation} \tau _{\rm SB}^{\rm St} = \frac{2}{9} \frac{\rho _{\rm SB}}{\rho _{\rm g}}\frac{R_{\rm SB}^2}{\nu }. \end{equation} \tag{ 5 }$$

Despite the fact that when $R_{\rm SB}> \mathcal {M}^{-1}\lambda$, the simple drag law (Equation (3)) no longer describes the motion of the dust particles, most numerical approaches to these problems (see, e.g., Johansen et al. 2007; Armitage 2010; Nelson & Gressel 2010; Carballido et al. 2010, 2011) continue to use it because a more accurate description is computationally prohibitive. Here we will call bodies of approximately this size "boulders." The mean-free-path λ in an accretion disk varies with radius; e.g., according to the minimum mass solar nebula model λ ranges from ≈10 cm at approximately 1.5 AU to ≈10 m at 10 AU. Hence R_boulder ranges from ∼2 m in inner disk regions to ∼200 m at about 10 AU. A more accurate approximation of the motion of such particles is given by the Maxey–Riley equation (Maxey & Riley 1983), which assumes a spherical geometry. While the Maxey–Riley approach is appealing on fundamental grounds, it has yet to be used in simulations of fully developed turbulence.

A crucial and often used assumption is that all collisions have a sticking probability of unity. Indeed, under such an assumption, planetesimal growth under a wide range of disk conditions is sufficiently rapid that there is no loss to the central star. Clearly, however, the probability of sticking depends, among other things, on the collisional velocity, the material properties of the colliding bodies, the ambient temperature, and the relative particle size. It is a commonly accepted picture that for collisional velocities V_c above a certain threshold value, V_th ∼ 0.1–10 cm s⁻¹, particle agglomeration is not possible; and elastic rebound overcomes attractive surface and intermolecular forces (e.g., Chokshi et al. 1993). However, for bodies covered with ice, experimental (Blum & Wurm 2008) and theoretical (Wettlaufer 2010) studies of collisions between dust grains and meter-sized objects have elucidated the range of collisional velocities (which depends on the relative particle size) over which perfect sticking occurs. This latter work considers the basic role of the phase behavior of matter (phase diagrams, amorphs, and polymorphs) in leading to the so-called collisional fusion. In this fusion process, a physical basis for efficient sticking is provided through collisional melting/amphorphization/polymorphization and subsequent fusion/annealing to extend the collisional velocity range of sticking to ΔV_c ∼ 1–100 m s⁻¹ ≫ V_th, which encompasses both typical turbulent rms (root-mean-square) speeds and the velocity differences between boulders and small grains ∼1–50 m s⁻¹. Moreover, bodies of high melting temperature and multicomponent materials, such as silicon and olivine, can fuse in this manner depending on the details of their phase diagrams. Hence, in principle, the approach provides a framework for sticking from the inner to the outer nebula. Here, we explore the influence of such a range, ΔV_c, on the growth of a boulder in a simulated disk.

The fate of the boulder is studied from a reference frame fixed to it, while the gas flows around it. We provide an accurate description of the flow by performing a direct numerical simulation (DNS) with no-slip and non-penetrating boundary conditions on the boulder surface using a numerical technique called the Immersed Boundary Method (Peskin 2002). Hence, there is no ad hoc approximation involved in describing the mutual interaction between the boulder and the gas flow. However, at present, it is computationally prohibitive to solve for more than one boulder using this DNS scheme. Consequently, we focus our study on the flow mediated collisions between one boulder and many "effectively" point-sized dust grains whose sizes are much smaller than R_boulder. Our principal approximations in treating the motion of the dust grains are (1) to use Equation (3) and (2) to ignore the back-reaction of the dust grains onto the flow. Advected by the turbulent disk flow, the dust grains collide with the boulder and we compute the probability density function (PDF) of the normal component of the collisional velocity.

Our computational domain is a rectangular box divided into two equal parts (Figure 1). In the right half, fluid turbulence is generated by external forcing that is non-zero between the two dashed lines shown in Figure 1. The turbulence thus generated is moved toward the "boulder" by the action of a body force g in the direction of the arrow shown in the figure. This body force is responsible for generating a mean flow, which models the head-wind faced by a boulder—the circular object at the left half of the domain. The boundary layer around the boulder is fully resolved by imposing non-penetrating and no-slip boundary conditions using the immersed boundary method. After the flow has reached a stationary state, we introduce N_p = 2 × 10⁴ particles into the right half of the domain as depicted in Figure 1. The motion of these particles obeys the simple drag law,

$$\begin{equation} \frac{d {\boldsymbol v}_{\rm p}}{dt} = \frac{1}{\tau _{\rm p}}({\boldsymbol v}_{\rm p} - {\boldsymbol U}), \end{equation} \tag{ 6 }$$

with the characteristic drag time of the "dust particles" τ_p. As noted before and as is clear from context, no such assumption need be made for the boulder. The back-reaction from the dust grains to the gas is ignored. When a dust grain collides with the boulder it is removed from the simulation and a new dust grain is introduced in the right half of the domain. We use the Pencil Code⁵ in which the immersed boundary method was implemented by Haugen & Kragset (2010).

Zoom In Zoom Out Reset image size

The characteristic large-scale velocity is the root-mean-square velocity in the streamwise direction, $v_{\rm wind}\equiv \langle v^2_y\rangle ^{1/2}$. We always use the Reynolds number corresponding to the central solid body, defined by

$$\begin{equation} \mbox{Re}_{ \rm SB}\equiv v_{\rm wind}R_{\rm SB}/\nu. \end{equation} \tag{ 7 }$$

The Stokes number of the "dust particles" is defined by

$$\begin{equation} \mbox{St}\equiv \tau _{\rm p}/\tau _{\rm L}, \end{equation} \tag{ 8 }$$

where τ_L = L_y/v_wind, with L_y being the length of our domain along the streamwise direction; from right to left in Figure 1. By virtue of limiting our simulations to two dimensions, we can access a larger range of particle Reynolds numbers Re_SB, from 30 to 1000 with resolutions ranging from 128 × 512 to 512 × 2048 grid points. The surface of the boulder is resolved with 100–400 grid points.

The criterion for a collision is that the distance between a dust grain and the boulder becomes less than a grid point. After this collision, we remove the dust grain from the simulation. For Re_SB ≈ 1000 in Figure 4, we plot the PDF, P(v_n), of the component of the velocity of the dust grain normal to the surface of the boulder, v_n. At small v_n, $P(v_{\rm n}) \sim v_{\rm n}^2$ (Figure 4(a)) and at large v_n the falloff is ∼exp [ − (v_n/v₀)²] (Figure 4(b)). However, as shown in Figure 4(c), over the whole range it is difficult to fit the PDF with a Maxwellian distribution.

Zoom In Zoom Out Reset image size

Now we consider how the PDF changes as the Stokes and Reynolds numbers of the flow change. We vary the Reynolds number by changing the viscosity of the flow. Hence, a change in Reynolds number also changes v_wind, and this leads to a change in the Stokes number.⁶ Therefore, in our approach to the numerical treatment of the flow, it is not possible to perform a systematic study of the Reynolds number dependence of the PDF at the fixed Stokes number. However, in order to produce an effective treatment of such a circumstance, we present in Figure 5 the PDFs for different Reynolds numbers wherein the Stokes numbers are not too different from each other. We see that for small Re_SB the peak of the PDF lies very close to v_wind, but as Re_SB increases the peak moves to smaller velocities by only a very small amount. Although Re_SB changes by almost a factor of 20, the position of the peak (normalized by v_wind) only changes from 0.6 to 0.3. A more dramatic change is observed in the PDF as the Stokes number is changed from 0.5 to 0.1 when Re ≈ 29 is held fixed, as shown in Figure 6. In particular, the tail of the PDF at high v_n is severely cutoff as the Stokes number is decreased by a factor of five. One can understand this as follows: when the Stokes number decreases, the dust grains begin to follow streamlines and hence never collide with the boulder. The implication of this is clearly that a smaller Stokes number implies a smaller number of high-impact collisions. Nevertheless, the most striking result for the problem at hand is the insensitivity of the PDF to Re_SB and St.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Our simulations take place in the reference frame of the boulder. Although the boulder is also comoving with the local gas with velocity v_kepler, it is the head wind, which corresponds to v_wind in our simulations, that sets the velocity scale. The radius of the boulder is taken to be ≈10 m. The magnitude of the headwind in the disk is estimated to be v_wind ≈ 10⁻³v_kepler ≈ 3 × 10³ cm s⁻¹ (see, e.g., Armitage 2010, p. 130). In the astrophysical literature it is common to non-dimensionalize τ_p with Ω_kepler, the Keplerian frequency, to define the orbital Stokes number, St_kepler. Here, we use the largest eddy timescale τ_L = L_y/v_wind, to obtain St. These two Stokes numbers are related by

$$\begin{equation} \mbox{St}_{\rm kepler}= \mbox{St}\frac{\tau _L}{\tau _{\rm orb}}, \end{equation} \tag{ 9 }$$

where τ_orb is the characteristic timescale of the Keplerian orbit defined by

$$\begin{equation} \tau _{\rm orb}= \frac{R_{\rm orb}}{v_{\rm kepler}}, \end{equation} \tag{ 10 }$$

where R_orb is the orbital radius. Using the definition of the two timescales τ_L and τ_orb, we obtain the ratio of the two Stokes numbers to be

$$\begin{equation} \frac{\mbox{St}_{\rm kepler}}{\mbox{St}} = \frac{L_y}{R_{\rm orb}} \frac{v_{\rm kepler}}{v_{\rm wind}} \approx 10^{-5}, \end{equation} \tag{ 11 }$$

where we have used R_orb = 1 AU, L_y = 50R_SB ≈ 500 meter, and v_wind = ηv_kepler with η = 10⁻³. We have used the Stokes number ranging from St = 0.1 to 2 which in turn gives St_kepler ≈ 10⁻⁶ to 2 × 10⁻⁵. We use the same conventions used in the supplementary information of Johansen et al. (2007) to convert the value of St_kepler to a radius of the dust grain; this implies that our "dust particles" are of the size of tenths of millimeters or smaller. Clearly, the "dust particles" are smaller than hydrodynamic scales, and hence it is justified to consider them as point objects whose motions are described by the Epstein drag law.

The PDFs of collisional velocities show that, irrespective of the Reynolds number and the Stokes number within the range considered by us, most collisions occur at velocities rather near to v_wind. To illustrate this, in Figure 5 we have drawn two vertical dashed lines at (v_n/v_wind)² = 0.2 and v_n/v_wind = 1.2. The area under the PDF between the two lines includes approximately 95% of the total number of collisions. Translated to parameters in the disk, this implies that, if there is a mechanism by which dust grains with velocities ranging from 0.4 v_wind to 1.1 v_wind would stick to a boulder, then we could consider 95% of collisions to have a perfect sticking probability.

Roughly speaking, this implies a range of velocities 6–36 m s⁻¹. These collisional velocities are far too high for the bodies to fuse by attractive intermolecular forces. An alternative scenario by which the colliding bodies can fuse at high speed has been suggested by Wettlaufer (2010). As discussed in Section 1.3, the very high local pressures that occur during a collision can lead to phase change. If, when the pressure begins to relax during rebound the momentarily liquified (or disordered) interfacial material re-freezes (or anneals) before particle separation, then fusion can occur. The idea was demonstrated when the colliding bodies are covered by ice, but the theory is generally applicable to all materials whose phase diagram is known in detail. An example of the process in a high melting temperature material (silicon) was noted in Wettlaufer (2010). Hence, whether the range of collisional velocities over which such process can occur in a material such as olivine matches with the range we find here is a topic of ongoing research. Note that here the particle Reynolds number Re_SB varies linearly with the particle radius but the range over which most of the collisions occur does not depend sensitively on Re_SB, and hence not on the particle radius. Thus, runaway growth of the boulder through the accretion of dust grains is a viable mechanism in areas of the disk where collisional fusion can operate in the range we obtain.