Nothing Special   »   [go: up one dir, main page]

A publishing partnership

Articles

GLOBAL DRAG-INDUCED INSTABILITIES IN PROTOPLANETARY DISKS

Published 2013 July 8 © 2013. The American Astronomical Society. All rights reserved.
, , Citation Mir Abbas Jalali 2013 ApJ 772 75 DOI 10.1088/0004-637X/772/1/75

0004-637X/772/1/75

ABSTRACT

We use the Fokker–Planck equation and model the dispersive dynamics of solid particles in annular protoplanetary disks whose gas component is more massive than the particle phase. We model particle–gas interactions as hard sphere collisions, determine the functional form of diffusion coefficients, and show the existence of two global unstable modes in the particle phase. These modes have spiral patterns with the azimuthal wavenumber m = 1 and rotate slowly. We show that in ring-shaped disks, the phase-space density of solid particles increases linearly in time toward an accumulation point near the location of pressure maximum, while instabilities grow exponentially. Therefore, planetesimals and planetary cores can be efficiently produced near the peaks of unstable density waves. In this mechanism, particles migrating toward the accumulation point will not participate in the formation of planets, and should eventually form a debris ring like the main asteroid belt or classical Kuiper Belt objects. We present the implications of global instabilities to the formation of ice giants and terrestrial planets in the solar system.

Export citation and abstract BibTeX RIS

1. INTRODUCTION

Protoplanetary disks are multi-phase environments composed of solid particles and molecular gas. The motion of particles is mainly governed by the gravitational forces of the disk material and the central star. Gas molecules feel the pressure gradient as well: for a polytropic isothermal gas, whose density profile monotonically increases toward the central star, the pressure gradient ∇p is always negative. This yields a sub-Keplerian circular velocity and generates headwind on solid particles that move on Keplerian orbits. Solid particles are thus expected to inspiral toward the central star. Although adhesive and electrostatic forces can enhance the clustering of micron-sized particles (Blum & Wurm 2008), centimeter- to meter-sized particles seem to inspiral toward (and fall into) the central star sooner than the timescale needed for assembling planetary cores. Therefore, several collective processes like streaming instability (Youdin & Goodman 2005), turbulent vortices induced by Kelvin–Helmholtz instability (Johansen et al. 2006; Bai & Stone 2010a, 2010b), and magnetorotational instability (Johansen et al. 2007, 2011) have been proposed to be responsible for the formation of kilometer-scale planetesimals.

The inspiraling motion of particles, however, does not globally occur in disks with non-monotonic density profiles, which are likely to form through a combination of viscous accretion and photoevaporation by the central star (e.g., Matsuyama et al. 2003). A ring-like disk is the simplest system with non-monotonic density profile. An interesting property of such systems is that ∇p becomes positive in regions where the density profile is rising, and gas molecules move with super-Keplerian velocities. Consequently, solid particles that approximately move on Keplerian orbits are accelerated by gas and migrate outward. This means that solid particles do not necessarily fall into the central star and may instead migrate to regions where the pressure is maximum and both the head and tail winds are minimized (Haghighipour & Boss 2003). Migrations of individual particles have been well understood by solving their equations of motion in the presence of gravitational and drag forces, but we do not still know the collective effects of such migrating particles. Can they efficiently produce planetesimals and massive planetary cores?

In this paper, we generalize the analysis of Jalali & Tremaine (2012, hereafter JT12) to disks that include a locally isothermal gas component, and search for global instabilities in the particle phase. Our disks are self-gravitating, and the particle phase has non-zero radial velocity dispersion. The dynamics of particles is modeled by the Fokker–Planck equation, and the gas component is assumed to be in a steady-state rotation around the central star. For simplicity, we confine our study to disks with Σgp ≫ 1, where Σg and Σp are the surface densities of the gas and particle phases, respectively. By this assumption, the background gas component does not respond to the disturbances of the particle phase. We neglect collisions between solid particles, but those between gas molecules and solid particles are taken into account.

We present our simple model of protoplanetary disks in Section 2 and derive the circular velocities of solid particles and gas molecules. In Section 3, we model the dispersive dynamics of particles in the context of kinetic theory, and utilize a perturbation theory in Section 4 to solve the resulting Fokker–Planck equation. In Section 5, we apply our theory to planet formation in the solar system, and in Section 6 we estimate the physical ranges of parameters for which the perturbation solutions are valid. We conclude the paper in Section 7 by comparing our findings with previous works. Open problems for future studies are also discussed.

2. DISK MODEL

We assume a two-phase medium consisting of solid particles and gas molecules, and refer to them by subscripts p and g, respectively. Solid particles are assumed to be monodisperse hard spheres of mass mp and radius rp. The gas phase has the molecular mass mg, and the average radius of gas molecules is rg. We work with initially axisymmetric disks whose gas component has no streaming motion in the radial direction (no accretion), but particles can move on eccentric orbits. For both the particle and gas phases we use the annular ring model of JT12, whose dimensionless surface density is

Equation (1)

where Md = Mp + Mg is the total disk mass, M is the mass of the central star, Mp is the mass of the particles, Mg is the mass of the gas component, b is a length scale, and r is the radial distance to the central star. The top panel in Figure 1 shows the radial profile of Σ0/μ that peaks at $R=\sqrt{6}/3$. We suppose that the surface densities of the particle and gas phases are proportional to Σ0 so that $\Sigma _{\nu }(R)=\bar{M}_{\nu } \Sigma _{0}(R)$ with ν ≡ p, g and $\bar{M}_{\nu }=M_{\nu }/M_{\rm d}$. Defining G as the gravitation constant, the dimensionless total gravitational potential field arising from Σ0 and the central star becomes

Equation (2)

where Φ0 is the actual potential.

Figure 1. Refer to the following caption and surrounding text.

Figure 1. Top: density profile of ring-shaped disks. Bottom: relative circular velocity of particles with respect to gas, normalized to the mean thermal speed of gas molecules. The values of q have been extracted from Table 1.

Standard image High-resolution image

If the gas is locally isothermal, its dimensionless pressure will be determined from

Equation (3)

where cs is the normalized sound speed, kB is the Boltzmann constant, and Tg is the absolute gas temperature. In passive disks heated by the stellar radiation (not by accretion), and sufficiently far from the central star, the radial profile of Tg is given by (Armitage 2010, Section 2.4)

Equation (4)

where T is the effective temperature of the central star and r is its radius. The mean thermal speed of gas molecules is related to the sound speed as $v_{\rm th}^2=(8/\pi)c_s^2$. We therefore find

Equation (5)

With mg being the molecular mass of H2, we have computed the value of q, given in Table 1 for several choices of b in the solar system and around AB Aurigae (Hashimoto et al. 2011) and Fomalhaut. The reason for our special selection of b will be explained in Section 4.

Table 1. The Parameter q for Several Choices of b in Three Protoplanetary/Planetary Systems

  M/M r/r T/T b (AU) q
Solar system 1 1 1 3.8 0.023
Solar system 1 1 1 72 0.034
Fomalhaut 1.92 1.82 1.486 243 0.043
AB Aurigae 3.1 2.1 1.3−1.7 69 0.031

Download table as:  ASCIITypeset image

In the absence of collisions between solid particles and gas molecules, the velocity of solid particles on circular orbits is determined from

Equation (6)

and the radial momentum equation for the gas becomes

Equation (7)

from which we obtain the circular velocity of gas:

Equation (8)

The condition $v_{{\rm p},c}^2\ge 0$ implies μ ⩽ 55/2, which is satisfied by protoplanetary disks. The bottom panel in Figure 1 illustrates the radial variation of Δvc = (vp, cvg, c)/vth for three values of q. The profile of Δvc is almost flat for R ≲ 0.3. According to Equation (8), the gas streaming velocity exceeds the speed of particles for R2 < 1/3 and generates tailwind on them. This is a remarkable feature of ring-shaped disks, and has interesting implications for the dynamics of solid particles in protoplanetary disks: while particles experience a resistive headwind for R2 > 1/3 and inspiral toward the central star, they are accelerated when R2 < 1/3 and migrate outward. Inward and outward migrating particles will then be accumulated near R2 ≈ 1/3 where the gas pressure is maximum (see Haghighipour & Boss 2003). In the next sections, we show that such migrations are not toward the exact location of pressure maximum if solid particles move on eccentric orbits. Moreover, such migrations are shown to be accompanied by exponentially growing instabilities.

3. EVOLUTIONARY DYNAMICS OF SOLID PARTICLES

The dynamics of particles is described by the phase-space distribution function (DF) $f({\textbf {x}},{\textbf {v}},t)=m_{\rm p} \, {\cal N}_{\rm p}({\textbf {x}},{\textbf {v}},t)$, where ${\cal N}_{\rm p}$ is the number density of particles. The vectors ${\textbf {x}} =(x_1,x_2)$ and ${\textbf {v}} =(v_1,v_2)$ (both in Cartesian coordinates) are the position and velocity vectors of particles in the disk plane, and t is the time. We utilize the DF (Jalali & Tremaine 2012, Section 3),

Equation (9)

to model the initial distribution of particles, before turning on the collisions between solid particles and gas molecules. Here, $L=\vert {\textbf {x}} \times {\textbf {v}} \vert$ and $E=\frac{1}{2} {\textbf {v}} \cdot {\textbf {v}} +V_0$ are, respectively, the orbital angular momentum and energy of particles per unit mass. There is an invertible, one-to-one and onto map from the (E, L)-space to the space of orbital elements (a, e), where a and e are the orbital semimajor axis and eccentricity, respectively. K is a positive integer that controls the mean eccentricity $\bar{e}$ of the particle disk. The function gK(E) takes physical (positive) values for K ⩾ 2, and $\bar{e}$ decreases as K is increased. In the limit of K the disk becomes cold with all particles moving on circular orbits.

In this study we ignore particle–particle collisions, and assume that colliding solid particles and gas molecules are hard spheres. The evolution of f is therefore expressed by the Fokker–Planck equation (Binney & Tremaine 2008, Section 7.4)

Equation (10)

where Dvi] and DviΔvj] are diffusion coefficients and ai are the components of the acceleration vector. In Equation (10) and throughout the paper a repeated integer index stands for summation over that index from 1 to 2. The acceleration vector is computed from ${\textbf {a}} =-\nabla V=-\nabla [V_0+V_1(R,\phi,t)]$, where the perturbed potential V1 is self-consistently calculated from the density disturbance Σ1(R, ϕ, t) of particle distribution. Since we have assumed Σp ≪ Σg, the contribution of the gas component to V1 is neglected.

For local collisions, the diffusion coefficients are evaluated using the procedure of Rosenbluth et al. (1957), but for three-dimensional collisions of hard spheres with the cross section ${\cal S}_{\rm p,g}=\beta ^2/4$, where β = rp + rg. Let us define γ = mg/(mp + mg) and denote the absolute velocity of gas molecules by ${\textbf {v}}^{\prime }$. We obtain

Equation (11)

where the potential functions $h({\textbf {x}},{\textbf {v}},t)$ and $g({\textbf {x}},{\textbf {v}},t)$ are given by the integrals

Equation (12)

Equation (13)

and ${\cal N}_{\rm g}({\textbf {x}},{\textbf {v}}^{\prime },t)$ is the number density of gas molecules in the phase space. In deriving Equations (12) and (13), we have assumed three-dimensional scattering of gas molecules by solid objects whose motion is confined to the disk plane. Therefore, the velocity vector of gas molecules is ${\textbf {v}}^{\prime }=(v^{\prime }_1,v^{\prime }_2,v^{\prime }_3)$, where $v^{\prime }_3$ is the velocity component perpendicular to the disk plane. Defining

Equation (14)

the elements of the stress tensor are determined as $\tau _{ij}= \Sigma _{\rm p} (\overline{v_i v_j} - \bar{v}_i \bar{v}_j )$. The macroscopic quantities Σp, $\bar{v}_i$, and τij are functions of ${\textbf {x}}$ and t, and the second term on the right-hand side of Equation (10) integrates to zero up to the first-order moment equations. Collisional terms involving $g({\textbf {x}},{\textbf {v}},t)$ correspond to random motions, and contribute to the evolutionary equations of τij (second-order moments of the Fokker–Planck equation). We neglect them in the present study because the mass ratio of gas molecules to solid particles is small, mg/mp ≪ 1, which implies $g/h \sim {\cal O}(m_{\rm g}/m_{\rm p})$. The Fokker–Planck equation can therefore be reduced to

Equation (15)

with the dimensionless motion equations

Equation (16)

Differentiating Equation (12) with respect to vi gives

Equation (17)

Evaluating this integral requires the explicit form of ${\cal N}_{\rm g}$. Nonetheless, such details are not necessary if we make some further simplifying assumptions: let $({\textbf {e}} _R,{\textbf {e}} _{\phi })$ be unit base vectors in the polar coordinates (R, ϕ). For e ≪ 1, the velocities of particles and gas molecules will be approximated as

Equation (18)

where $\Omega =R^{-3/2}+{\cal O}(\mu)$ is the orbital frequency of particles. We think of disks with vth < vp, c, vg, c over 0.03 ≲ R ≲ 10 (cf. Figure 1) to guarantee the existence of bound orbits, and avoid gas dispersal through photoevaporation. The mean eccentricity corresponding to Equation (9) is almost constant over the entire disk space (see JT12), and it is given by $\bar{e}=[\pi /(4K+10)]^{1/2}+{\cal O}(\mu)$. For $q \sim {\cal O}(10^{-2})$, which corresponds to protoplanetary disks around solar-type stars, and for K ⩽ 29 experimented in JT12, one has

Equation (19)

From Equations (8), (18), and (19) we conclude that the bulk of particles move with supersonic speeds with respect to the gas stream, and they satisfy $|{\textbf {v}}-{\textbf {v}}^{\prime }|\approx e R \Omega$. The integral in Equation (17) is thus approximated by

Equation (20)

Equation (21)

where ρg = Σg(R)δ(z) and $U_i({\textbf {x}},t)=\bar{v}^{\prime }_i$ are the spatial density and streaming velocity components of the gas phase, respectively. Here, δ(z) is Dirac's delta function and z measures the height above the disk mid-plane. The gas density and velocity components have been normalized to M/b3 and [GM/b]1/2, respectively. Equation (20) is equivalent to the drag force expression of Kwok (1975) in supersonic regimes. Assuming that the mass of each particle is computed from $m_{\rm p}=(4/3)\pi \rho _{\rm s} r_{\rm p}^3$, with ρs being the typical density of rocky material, one finds ξ0∝1/rp. With ξ0 = (M/mp)(rp/b) the scattering of gas molecules takes place in the disk plane and the cross section ${\cal S}_{\rm p,g}$ becomes a line of the length rp + rg. We are not interested in this extreme unphysical case.

For near-circular orbits with e → 0, one obtains $|{\textbf {v}}-{\textbf {v}}^{\prime }|\approx v_{\rm th}$ and Dvi] transforms to the well-known form of Epstein drag:

Equation (22)

During our numerical computations we use Equation (20) if evp, c > vth and apply Equation (22) otherwise. A factor 4/3 is missing on the right-hand side of Equation (22). It can be recovered through assuming a Maxwell–Boltzmann distribution in the velocity space for ${\cal N}_{\rm g}$, and exactly performing the integral in Equation (17). Nonetheless, the missing factor is unimportant in our computations, for we will vary ξ0 to explore the influence of drag force on the disk evolution, and one may suppose that any constant factors have already been included in ξ0. To compute diffusion coefficients numerically, we soften Dirac's delta function using its normal distribution representation:

Equation (23)

where h ≪ 1 is the dimensionless scale height of the disk and can be a function of R. The three-dimensional structure and evolution of circumstellar disks have not been modeled in this study; we thus work with a constant h and set

Equation (24)

in the disk mid-plane where the Fokker–Planck equation governs the evolution of the particle phase.

Solving Equation (15) in a four-dimensional phase space, with particle motions confined to the disk plane, is facilitated by utilizing the angle variables ${\textbf {w}} =(w_1,w_2)$ and their conjugate actions ${\textbf {J}} =(J_1,J_2)$. We follow JT12 and set J1 and J2 to the radial action JR and the orbital angular momentum Jϕ = Rvϕ, respectively. In the $({\textbf {w}},{\textbf {J}})$-space, the motion equations (16) become

Equation (25)

where ${\cal H}=({1}/{2}) {\textbf {v}} \cdot {\textbf {v}} +V({\textbf {x}},t)$ is the Hamiltonian function, and the generalized forces $F_{J_i}$ and $F_{w_i}$ are determined using the virtual work of nonconservative forces:

Equation (26)

with δ being the variational operator. In the Appendix, we explain the procedure of calculating $F_{w_i}({\textbf {J}},w_1)$ and $F_{J_i}({\textbf {J}},w_1)$. They are real harmonic functions of w1 and are smooth in the ${\textbf {J}}$-space. We now write Equation (15) in the angle–action space:

Equation (27)

where [ ·  , ·] denotes the Poisson bracket over the $({\textbf {w}},{\textbf {J}})$-space and the collision operator ${\cal D}_F$ is defined by

Equation (28)

4. UNSTABLE MODES

We seek solutions of the form $f=f_0({\textbf {J}})+f_1({\textbf {w}},{\textbf {J}},t)$ for Equation (27) so that |f1| ≪ |f0|. The Hamiltonian corresponding to f will become ${\cal H}={\cal H}_0+V_1$, where the perturbed potential V1 (self-consistently arising from f1) and ${\cal H}_0=({1}/{2}) {\textbf {v}} \cdot {\textbf {v}} +V_0(R)$ are expressed in the angle–action space (Jalali 2010, Section 2). f1 is obtained by solving the perturbed Fokker–Planck equation:

Equation (29)

which is a non-homogenous partial differential equation.

4.1. Secular Migrations

The particular solution of Equation (29) is a radial drift of the form $f_d({\textbf {J}},w_1,t)$. For the small disturbances |fd| ≪ |f0|, we can ignore ${\cal D}_F f_d$ against ${\cal D}_F f_0$ and write

Equation (30)

where the left-hand side is the linear approximation of the total derivative dfd/dt, and the potential Vd corresponds to fd. From the definition of ${\cal D}_F$ and Equations (A8)–(A10), we arrive at

Equation (31)

whose $l\not= 0$ terms lead to particular periodic solutions $f_d({\textbf {J}},w_1)=f_d({\textbf {J}},w_1+2\pi)$. However, the dominant collisional term

Equation (32)

results in a secular drift fdS0t in the phase space: the DF of particles linearly increases in time if their actions satisfy $S_0({\textbf {J}}) > 0$, and it decreases for $S_0({\textbf {J}})<0$. Such drag-induced migrations can accumulate particles in regions where $S_0({\textbf {J}})$ has a positive local maximum. Since S0 (and therefore dfd/dt) depends on the initial DF, it is useful to normalize it to f0 and investigate the relative variation of particle distribution in the (a, e)-space.

For a model with K = 29 and q = 0.034, we have plotted the contours of S0/(f0ξ) in Figure 2 for mass parameters μ = 0.04 and μp = Mp/Mg = 1/8. The mean eccentricity of particle orbits of this model is $\bar{e}=0.156$. It is seen that S0/(f0ξ) is negative for eccentric orbits with e ≳ 0.22 and the population of those orbits is falling in time. Meanwhile, the number of particles increases toward an accumulation point at (a0, e0) ≈ (0.63, 0.128). The eccentricity of this accumulation point is less than the mean eccentricity of the initial model, and its semimajor axis is very close to R = 0.577, where the gas pressure is maximum. The local maximum of S0/(f0ξ) in the (a, e)-space is reminiscent of the distribution of asteroids between Mars and Jupiter with the mean orbital eccentricity ≈0.14, and classical Kuiper Belt objects (KBOs) with $\bar{e}_{\rm KBO} \approx 0.1$. The mean semimajor axis of classical KBOs is $\bar{a}_{\rm KBO} \approx 45 \, {\rm AU}$. If we assume that they are remnants of planet formation that reside at the accumulation point, the length scale b of our annular disk model reads $b= \bar{a}_{\rm KBO}/a_0 \approx 71.43 \, {\rm AU}$, which gives q ≈ 0.034 used in computations of Figure 2 (see also Table 1). Moreover, for the asteroids between Mars and Jupiter the semimajor axis ranges from aAST ≈ 2.1 to 3.3 AU. The most populous group of these asteroids has a mean semimajor axis of $\bar{a}_{\rm AST} \approx 2.4 \, {\rm AU}$. The second plausible length scale of our model is therefore $b=\bar{a}_{\rm AST}/a_0 \approx 3.81 \, {\rm AU}$ that corresponds to q ≈ 0.023. Similarly, the values of b and q given in Table 1 for the Fomalhaut and AB Aurigae systems have been calculated based on the observed semimajor axes of their dust rings. Decreasing μ and μp does not considerably change the pattern of S0/(f0ξ), but proportionally decreases its maximum value at the accumulation point.

Figure 2. Refer to the following caption and surrounding text.

Figure 2. Isocontours of the collisional term S0/(f0ξ) in the space of orbital elements. Horizontal dashed line corresponds to the mean eccentricity $\bar{e}=0.156$ of particle orbits in the initial model. Mass parameters have been set to μ = 0.04 and μp = 1/8. Particles migrate linear in time toward the accumulation point (a0, e0) ≈ (0.63, 0.128) shown by a cross. Contour lines mark the levels 0.0 to 4 × 10−4 with steps of 10−4. Note the logarithmic scale of the a-axis.

Standard image High-resolution image

4.2. Exponentially Growing Instabilities

We now search for non-axisymmetric homogeneous solutions of Equation (29) that depend on (w1, w2, t) and satisfy

Equation (33)

We consider unsteady DFs of the form (JT12, Section 6)

Equation (34)

which corresponds to a slowly rotating density wave with the azimuthal wavenumber m:

Equation (35)

Here, ωr/m and s are the pattern speed and growth/decay rate of density perturbations, respectively. The radial profiles of the wave amplitude A(R) and phase angle ϑ(R) are time-invariant in the linear regime, and ϑ vanishes for stable waves with s = 0. The perturbed potential

Equation (36)

and density Σ1 are related through Poisson's integral. In this paper we work with lopsided m = 1 waves, which accelerate the central star. Therefore, the reference frame attached to the central star is not inertial, and we include the indirect gravitational potential (see JT12) in our formulation.

Substituting Equations (34) and (36) into Equation (33) yields the following linear eigenvalue problem for ω and its associated eigenfunction $\tilde{f}_1$:

Equation (37)

where $\varpi ({\textbf {J}})=\Omega _2({\textbf {J}})-\Omega _1({\textbf {J}})$ is the precession frequency of particle orbits. The orbital frequencies $\Omega _i({\textbf {J}})=\partial {\cal H}_0/\partial J_i$ (i = 1, 2) are determined in the unperturbed state. For nearly circular orbits one finds (JT12)

Equation (38)

which has the maximum value ϖmax = 0.05861μ at R = 0.2859. We utilize the finite element method of Jalali (2010) and JT12, and compute the eigenfrequency spectrum of Equation (37) for a model with K = 29 and q = 0.034. We vary μ, μp, and ξ to investigate the effects of particle size and mass fraction on the development of density waves. Our finite element model has N = 200 ring elements whose radial nodes are located at Rn = 10(4n − 2N − 4)/(N + 1) for n = 1, 2, ..., N. Using this mesh, eigenfrequencies are calculated with a relative accuracy ⩽0.5%.

In the absence of gas drag, ξ = 0, and for 1/32 ⩽ μp ⩽ 1/8, the spectrum contains two stable slow modes, which are labeled by A1 (fundamental mode) and A2 (secondary mode). The pattern speeds of these modes satisfy the inequality ωr ≳ ϖmax. All other modes are singular and form a continuum over the range ϖmin < ωr < ϖmax. The minimum and maximum precession frequencies, ϖmin and ϖmax, occur at the circular orbit boundary of the action space with J1 = 0, and $\varpi ({\textbf {J}})$ vanishes for radial orbits. Singular modes are associated with the inner Lindblad resonance $\varpi ({\textbf {J}})-\omega _r =0$ and can engage both circular and non-circular orbits. Table 2 shows the variation of ϖmin and ϖmax in terms of μ and μp.

Table 2. Pattern Speeds and Growth Rates of Slow Modes A1 and A2 versus ξ

Model Model Parameters Mode A1 Mode A2
μ μp ϖmin × 103 ϖmax × 103 ξ ωrmax s × 105 ωrmax s × 105
1 0.04 1/8 −8.239 2.344 0 1.105 0 1.053 0
2 0.04 1/8 −8.239 2.344 0.10 1.105 1.719 1.053 1.908
3 0.04 1/8 −8.239 2.344 0.20 1.106 3.428 1.052 3.821
4 0.04 1/8 −8.239 2.344 0.40 1.110 6.777 1.051 7.638
5 0.04 1/16 −8.239 2.344 0 1.034 0 0.9995 0
6 0.04 1/16 −8.239 2.344 0.10 1.034 1.908 0.9997 1.848
7 0.04 1/16 −8.239 2.344 0.20 1.036 3.787 1.000 3.516
8 0.04 1/16 −8.239 2.344 0.40 1.040 7.386 ... ...
9 0.02 1/16 −4.137 1.172 0 1.034 0 0.9994 0
10 0.02 1/16 −4.137 1.172 0.10 1.034 0.954 0.9995 0.924
11 0.02 1/16 −4.137 1.172 0.20 1.035 1.893 1.000 1.758
12 0.02 1/16 −4.137 1.172 0.40 1.040 3.693 ... ...
13 0.01 1/32 −2.073 0.586 0 0.998 0 0.977 0
14 0.01 1/32 −2.073 0.586 0.10 0.999 0.494 ... ...
15 0.01 1/32 −2.073 0.586 0.20 1.001 0.947 ... ...
16 0.01 1/32 −2.073 0.586 0.40 1.005 1.756 ... ...

Note. In all models we have set K = 29 and q = 0.034.

Download table as:  ASCIITypeset image

By setting ξ > 0 and turning on the gas drag, all singular and long-wavelength modes become unstable. The pattern speeds and growth rates of modes A1 and A2 have been reported in Table 2 for 16 different models whose particle and gas phases are Toomre stable. It is seen that ωrmax is a function of μp, and s is an almost linear function of $\xi \bar{M}_{\rm g}=\xi \mu /(1+\mu _{\rm p})$. Gas drag has no notable contribution to the pattern speeds of modes and only controls the growth rate. Figure 3 displays the perturbed density patterns of stable model 1 and unstable model 4. Except in the model with μp = 1/8, mode A1 always rotates and grows faster than mode A2. Comparing the modal content of our ξ = 0 disks with the results of JT12 (see their Figure 4) shows that short-wavelength slow modes have disappeared by setting Mp < Md and the pattern speeds of the modes with the longest wavelengths have approached $\varpi _{\rm max}^{+}$. This is very similar to the behavior of disk galaxies: increasing the mass of the dark matter halo stabilizes tightly wound spiral modes, and only modes with the longest wavelengths, especially the bar mode, survive (Jalali 2007).

Figure 3. Refer to the following caption and surrounding text.

Figure 3. Mode shapes of global stable (top) and unstable (bottom) density waves of the particle phase in a protoplanetary disk with μ = 0.04, Mp/Mg = 1/8, and $\bar{e}=0.156$. Contour plots show the positive parts of the perturbed density Σ1 at t = 0. In all cases, density waves rotate counterclockwise, and their corresponding pattern speeds are given in Table 2. The outer density bump of unstable mode A1 is a trailing spiral. The inner and outer wave packets of unstable mode A2 lead and trail their stable counterparts, respectively. Solid circle shows the location of the central star.

Standard image High-resolution image

Figure 4 demonstrates the amplitude function A(R) for several models. It is seen that by increasing ξ the local minimum between two density maxima increases. This is because of the enhanced spirality that smoothly connects density bumps. The local minimum is exactly equal to zero in models with ξ = 0 and corresponds to a node of stable oscillatory waves. Varying μ does not change the mode shape (this had already been pointed out by JT12), but decreasing μp shortens the wavelength of both modes A1 and A2. Variation of q has a negligible effect on the eigenmodes. Our experiments show that reducing q from 0.034 to 0.023 changes the eigenfrequencies by ≈0.2%.

Figure 4. Refer to the following caption and surrounding text.

Figure 4. Variation of the amplitude profile A(R) of modes A1 and A2 in models with different mass parameters and drag coefficients. See Table 2 for the specifications of the models. The difference between models 6 and 10 is indistinguishable in the plots. The R-axis is in logarithmic scale.

Standard image High-resolution image

For discrete slow modes, resonant cavities become smaller as μp falls off and mode A2 disappears: mode A2 has hardly managed to exist in model 13, which is stable. Nonexistence of mode A2 in models with μp = 1/16 and ξ ≳ 0.2 shows that the development of unstable wave packets is allowed only for ξ < ξcr, where the critical drag parameter ξcr depends on both μ and μp. Also note the nonexistence of unstable secondary modes in models with μp = 1/32. We have a simple physical explanation for the existence of a critical drag parameter: slow modes are supported by the precession of orbits, and the leading or trailing nature of developing unstable spiral patterns can be estimated using dϖc/dR, which switches sign at R = 0.2859. Spiral patterns will be trailing if dϖc/dR < 0 and leading otherwise (cf. Binney & Tremaine 2008, Section 6.1.3). In the region occupied by the central wave packet of mode A2 (see Figure 4), and when the mode is unstable, the quantity dϖc/dR can take both negative and positive values and that wave packet will be sheared. The larger the value of ξ, the higher the imposed shearing. The wave packet can thus resist disruption only for small deformations corresponding to ξ < ξcr. It is evident from Figure 3 that the inner wave packet of mode A2 leads its stable counterpart because it lies in the region with dϖc/dR < 0. The opposite phenomenon is happening for the outer wave packet.

Only mode A1 can marginally tolerate ωrmax ≲ 1 and ξ > 0 because the local minimum of A(R) is located near the maximum of ϖc, and its two main wave packets lie in regions where dϖc/dR has a definite sign. Moreover, its resonant zone is large enough to trap non-circular orbits with $\varpi ({\textbf {J}}) < \varpi _{\rm max}$. The number of slow modes also depends on the mean eccentricity of particle orbits. Our numerical experiments with μp = 1/8 show that mode A2 completely disappears as its pattern speed drops below ϖmax by increasing $\bar{e}$.

Except for ξ = 0, $f_0({\textbf {J}})$ is not an equilibrium DF and it is regarded as the initial condition for the perturbed Fokker–Planck equation. It is hard to imagine an equilibrium state at early epochs of protoplanetary disks when most ingredients of planet formation are transported due to dissipative forces. There is indeed a competition between the homogenous and particular solutions of Equation (29), and the relative magnitude of s with respect to S0/f0 decides which process wins. According to Figures 2 and 3, near the major peaks of modes A1 and A2 we have S0/f0 ≈ 2 × 10−4ξ, which has exactly the same order of magnitude of s for unstable μp = 1/8 models of Table 2. Drag-induced instabilities that grow proportional to est can therefore overwhelm secular migrations, consume most solid particle reserves of the disk within R < a0, and rapidly form bigger objects. Near the accumulation point, the amplitude of unstable density waves diminishes significantly, but since secular migration is a very slow process (linear in time), particles are expected to form only a debris ring as the gas is depleted at the later stages of disk evolution. A small fraction of solid particles will eventually live in the vicinity of the accumulation point, and the majority of them are transported through spiral arms to unstable regions.

5. APPLICATIONS TO THE SOLAR SYSTEM

To this end, we discuss the implications of global drag-induced instabilities to planet formation in the solar system. We use the amplitude functions of modes in model 1 because their profiles do not change significantly by varying the mass parameters and ξ (see Figure 4).

According to simulations of circumstellar disks, viscous accretion, photoevaporation, and stellar winds create a gap structure near the gravitational radius (Matsuyama et al. 2003), and the disk is split into two annuli. For the solar nebula, the gravitational radius is between the orbits of Saturn and Uranus (Shu et al. 1993), and therefore the outer ring of the solar system would contain ice giants and KBOs. If we suppose that classical KBOs are the debris material near the accumulation point, our ring-like disks and their modal content can be fit to the structure of the outer solar system by assuming $b=\bar{a}_{\rm KBO}/a_0=71.43 \, {\rm AU}$ (top panel in Figure 5). Interestingly, the semimajor axis of Uranus matches the location of the major density bump of mode A1 and Neptune is clearly associated with the outer bump. If we assume that both planets were formed exactly at the peaks of mode A1, Uranus and Neptune should have migrated outward for about 1.5 and 6 AU, respectively. This is consistent with Malhotra's (1995) resonant capture theory that explains the orbital dynamics of Pluto and Plutinos.

Figure 5. Refer to the following caption and surrounding text.

Figure 5. Normalized amplitude profiles A(r)/Σ0(r) of unstable modes A1 and A2 for model 1 in Table 2. Top: classical KBOs reside at the accumulation point, and we have set $b=\bar{a}_{\rm KBO}/a_0=71.43 \, {\rm AU}$. The current semimajor axes of the orbits of Saturn and Neptune are indicated by vertical dotted lines. Bottom: unstable modes are fit to the inner solar system with $b=\bar{a}_{\rm AST}/a_0=3.81 \, {\rm AU}$. The mean semimajor axis of the most populous asteroid belt and the present semimajor axes of the terrestrial planets have been indicated by vertical dotted lines.

Standard image High-resolution image

If we now assume that the inner annulus (emerged from the gap formation) contained terrestrial planets and the main asteroid belt, and that the asteroids between 2.1 and 3.2 AU are the remnants of secular migrations, one can fit our annular disk to the structure of the inner solar system by setting $b=\bar{a}_{\rm AST}/a_0=3.81 \, {\rm AU}$. Doing so, present orbital semimajor axes of Venus, the Earth, and Mars will lie in the region affected by the density bumps of modes A1 and A2 (bottom panel in Figure 5). We note that the mean radial distance between the density bumps of unstable modes and the accumulation point is approximately equal to the distance from the position of maximum orbital precession (with dϖc/dR = 0) to the position of maximum pressure. This characteristic length depends on the radial variations of surface density, sound speed, and gravitational potential, but the relative positions of unstable modes and the accumulation point seem to be model independent because the maximum precession rate occurs where the surface density is rising, and the region with maximum gas pressure is close to the region with maximum surface density.

Our theory cannot be directly applied to the formation of gas giants, Jupiter and Saturn, because the gas disk was not responsive to perturbations in the particle phase. Nonetheless, we can make useful predictions about the possible origins of gas giants and see whether they could have interfered with the formation of rocky planets. First of all, Jupiter and Saturn lie well within the inner ring characterized by b = 3.81 AU, and since we have not detected any instability of the particle phase at their current orbital distances, they have probably formed from an instability in the gas phase. We compute Toomre's Q over the inner ring and in terms of dimensionless variables:

Equation (39)

where $\kappa =\Omega +{\cal O}(\mu)$ is the epicyclic frequency of near-circular orbits. The function Q(R) has a global minimum at Rcr = 1.856, and its minimum value is $Q_{\rm min} = 3.168 \, q/\bar{M}_{\rm g}$. Near the orbit of Jupiter we have q = 0.024. The gas phase thus becomes unstable there if Qmin < 1, which implies $\bar{M}_{\rm g} > 0.076$. This corresponds to a more massive disk than the minimum solar nebula and models investigated in Table 2. Nevertheless, the predicted mass threshold falls well in the mass range of circumstellar disks observed around nearby stars. The most striking point of this calculation is that when one sets $b=\bar{a}_{\rm AST}/a_0=3.81\, {\rm AU}$, short-wavelength instabilities in the gas phase can be triggered around rcr = Rcrb = 7.07 AU, which is halfway between the orbits of Jupiter and Saturn. For $\bar{M}_{\rm g} = 0.084$, equation Q(R) − 1 = 0 has two roots at r1 = 4.95 and r2 = 10.6 AU, and the entire region between the orbits of Jupiter and Saturn is unstable in Toomre's sense. Such short-wavelength instabilities will not affect inner regions where the particle phase is unstable and rocky planets are being assembled.

6. PHYSICAL RANGES OF PARAMETERS

In this section, we determine over which physical ranges of parameters the perturbation solutions of the Fokker–Planck equation are acceptable. Throughout the calculations of this section, we set M = M and use a density of ρs = 3 g cm−3 for solid particles.

We have ignored particle–particle collisions in writing Equation (29), and this simplification is legitimate if particles collide after several orbital periods. In a monodisperse system of spherical particles, the collisional cross section is ${\cal S}_{\rm p}=\pi (2r_{\rm p})^2$. Moreover, the radial velocity dispersion of particles in our disks is determined from

Equation (40)

with κ being the epicyclic frequency of near-circular orbits. In the absence of gas, the collision time tcollide (in the disk mid-plane) normalized to the orbital period torbit reads

Equation (41)

where hp is the scale height of the particle disk. In low-mass disks particles can be scattered up to a vertical distance of $h_{\rm p}=\sigma _R/\Omega +{\cal O}(\mu)$, and Equation (41) yields

Equation (42)

At R = 0.3, which is the position of the node of mode A1 and the mean orbital distance of its two density bumps, we obtain

Equation (43)

with $\bar{M}_{\rm p}=\mu _{\rm p}\mu /(1+\mu _{\rm p})$. Since particle–particle collisions can be ignored only for $\bar{t}_{\rm collide} \gg 1$, our governing equations are valid for all models of Table 2 and for b = 3.81 AU (this is the length scale of the inner solar system) if rp ≫ 1.2 m. This size threshold reduces to rp ≫ 3.5 mm by adopting b = 71.43 AU for the outer solar system. On the other hand, the mean free path of gas molecules is defined as (e.g., Blundell & Blundell 2006, Section 8.3)

Equation (44)

Note that ρg is the dimensionless spatial density of the gas disk (see Section 3). For molecular hydrogen, we have mg = 3.32 × 10−27 kg and the collision cross section is ${\cal S}_{\rm g} \approx 2\times 10^{-19}\, {\rm m}^2$. At R = 0.3 and in the disk mid-plane we obtain

Equation (45)

Assuming a scale height h = 0.01 (e.g., Goldreich & Ward 1973) in the models of Table 2, we find λ ≈ 4–16 cm for b = 3.81 AU and λ ≈ 300–1000 m for b = 71.43 AU. From the acceptable values of rp for having a collisionless particle phase and the physical range of λ, we conclude that particles interact with gas molecules through Stokes drag (skin friction) if we apply our model to the inner solar system. In such a circumstance, one must use the following drag parameter:

Equation (46)

where the drag coefficient CD depends on the Reynolds number. In the outer solar system, the drag force is computed from Equations (20) and (22) up to kilometer-size objects; it then switches to Stokes drag. When particles move with supersonic speeds with the Mach number ${\cal M}\sim eR\Omega /c_{\rm s} \gtrsim 3$, the drag coefficient approximately becomes CD ≈ 0.92 (Liu 1958).

Perturbation theory fails for large values of drag force, and one needs to constrain particle sizes (to which our results are applied) by the value of ξ. From Equation (46) one can write

Equation (47)

Consequently, the maximum value of ξ reachable by perturbation theory puts a minimum threshold on the allowed particle sizes. Our numerical experiments show that by increasing the drag parameter to ξmax ∼ 1 the accuracy of mode A1 drops significantly and the amplitude function A(R) of that mode loses its smoothness, especially for a smaller fraction of solid particles. Using this empirical upper limit of ξ, and with h = 0.01, CD = 0.92, and b = 3.81 AU, we require rp ≳ 28 km to ensure the validity of perturbation theory. A nonlinear Fokker–Planck equation solver is thus needed to understand the physics of global instabilities for subkilometer- and kilometer-sized particles in the inner solar system. Choosing CD = 2 (collisional/Epstein drag regime) and b = 71.43 AU, we obtain rp ≳ 170 m, which is smaller than the mean free path of gas molecules in the outer solar system.

For the valid ranges of rp discussed above, one can readily verify that the dimensionless stopping time parameter

Equation (48)

satisfies τs ≫ 1 in regions affected by unstable modes. The stopping time tstop in Equation (48) has been defined based on Epstein drag law. One still finds τs ≫ 1 if Stokes drag force applies. Therefore, solid particles are not dynamically coupled to the gas flow, and modeling the collective dynamics of particles in the context of kinetic theory is justified.

7. DISCUSSIONS

The infall timescale of ${\cal O}(10^2) \, {\rm yr}$ for meter-sized solid bodies puts a strong constraint on the formation of planetesimals from dust grains and pebbles (Weidenschilling 1977). Recent simulations have shown that the back-reaction of particles on gas can trigger out-of-plane Kelvin–Helmholtz instability, which boosts local particle density and helps self-gravity assemble kilometer-sized bodies. Turbulence, on the other hand, imposes stochastic forces on planetesimals, increases their collision frequency, and disrupts those with ⩽10 km radius (Nelson & Gressel 2010). Although the random forcing of planetesimals can be suppressed in the presence of a dead zone (Gressel et al. 2011), alternative and simpler processes may also be involved in the formation of super kilometer-scale planetesimals.

Slow density waves exist in all near-Keplerian, self-gravitating rings and can be excited by encounters (JT12). Addition of a gas component, however, destabilizes the particle phase without any external disturbance. We showed that the drag-induced infall of solid bodies into the central star is not a universal phenomenon and the direction of particle migration highly depends on the disk structure. There will be no infall if at some region the disk density profile, including its solid particle and gas components, rises outward. Therefore, with a preserved source of solid particles in the disk, global instabilities explored in this study will have time to boost the particle density to arbitrarily large levels and enhance the formation of bigger objects through gravitational collapse. One of the fundamental achievements of this study was how secular migrations and global instabilities can be used to identify possible planet-forming regions in observed protoplanetary systems and debris disks.

Although the gas flow in our disks was in the laminar regime, turbulence does not seem to considerably change our fundamental results. Youdin (2011) estimates turbulent eddy length as

Equation (49)

where α is the dimensionless turbulent diffusivity. For α ≲ 10−3 used by Youdin (2011), we see that ℓeddy is smaller, at least by three orders of magnitude, than the scale of the wave packets (and therefore the wavelength) of unstable modes. Therefore, turbulent diffusion is unimportant in the development of global density waves and their growth.

Michikoshi et al. (2010) have also studied the formation of planetesimals through gravitational instability. They assumed a non-responsive gas component, as we did, and introduced fluid dynamical equations to model the dynamics of dust particles in a local simulation box that rotates with Keplerian angular velocity. They then derived a dispersion relation by linearizing the continuity, momentum, and Poisson equations, and showed the existence of secularly unstable long-wavelength modes for a dissipative dust layer. This is somehow consistent with our findings that unstable modes have a long-wavelength nature. However, their perturbation theory and N-body simulations that utilize a rotating simulation box with periodic boundary conditions are not able to provide a global picture of particle migrations and information about possible planet-forming regions. Moreover, the analytical results of Michikoshi et al. (2010) are valid only for particles dynamically coupled to gas; otherwise, fluid dynamical equations with an isotropic pressure tensor could not be applied to the particle phase. Our findings apply to dispersive particle disks where the evolution of orbital eccentricity does matter.

Due to computational difficulties of working with large values of K in f0 and a small mass ratio μp, we investigated global modes only for K = 29 that gives $\bar{e}=0.156$. In accordance with WKB theory (JT12), decreasing the mean eccentricity of the particle disk is expected to preserve global stable modes, which then bifurcate to unstable density waves in the presence of gas drag. For smaller mean eccentricities, the number of global modes may even increase as the resonant cavities of slow modes become thinner. A useful future experiment would be to decrease $\bar{e}$ using a Schwarzschild DF and determine the number and shape of exponentially growing modes for μp ≲ 0.01.

Our results have been obtained for self-gravitating disks whose particle phase is constituted from monodisperse hard spheres. Including the size distribution of particles and the effect of particle–particle collisions, especially catastrophic disruptions, is an interesting open problem. Moreover, by assuming Σp ≪ Σg we neglected the back-reaction of particles on gas flow. A more accurate procedure is to simultaneously perturb the hydrodynamic and Fokker–Planck equations for the gas and particle phases, respectively. We anticipate angular momentum exchange between the particle and gas phases, and any particle migration should induce radial mass transfer in the gas phase.

I thank Scott Tremaine for his stimulating discussions during the course of this project. I also thank the referee for useful comments that inspired me to carry out new computations and improve the presentation of the paper.

APPENDIX: GENERALIZED FORCES IN THE ANGLE–ACTION SPACE

We define (vR, vϕ) as the velocity components of particles in the polar coordinates (R, ϕ), and the streaming velocity of the gas component will become $U_{\phi }(R) {\textbf {e}} _{\phi }$, where Uϕ(R) has the functional form of vg, c. Using the action variables ${\textbf {J}} =(J_1,J_2)$ and their conjugate angles ${\textbf {w}} =(w_1,w_2)$, the radial distance and azimuthal position of a test particle are calculated from

Equation (A1)

where the Fourier coefficients are given by

Equation (A2)

These integrals are taken over a full cycle of rosette orbits. From Equation (A1) one can compute the variations δR and δϕ as

Equation (A3)

Equation (A4)

The virtual work of the drag force reads

Equation (A5)

where CD(R) = ξ0evp, c ρg if the orbital eccentricity e satisfies the inequality e > vth/vp, c and CD(R) = ξ0vthρg otherwise. We now utilize the following Fourier expansions:

Equation (A6)

Equation (A7)

and obtain

Equation (A8)

Equation (A9)

Equation (A10)
Please wait… references are loading.
10.1088/0004-637X/772/1/75