GLOBAL DRAG-INDUCED INSTABILITIES IN PROTOPLANETARY DISKS

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 Σ_g/Σ_p ≫ 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.

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 m_p and radius r_p. The gas phase has the molecular mass m_g, and the average radius of gas molecules is r_g. 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

where M_d = M_p + M_g is the total disk mass, M_⋆ is the mass of the central star, M_p is the mass of the particles, M_g 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 Σ₀/μ that peaks at . We suppose that the surface densities of the particle and gas phases are proportional to Σ₀ so that with ν ≡ p, g and . Defining G as the gravitation constant, the dimensionless total gravitational potential field arising from Σ₀ and the central star becomes

where Φ₀ is the actual potential.

Zoom In Zoom Out Reset image size

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

where c_s is the normalized sound speed, k_B is the Boltzmann constant, and T_g 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 T_g is given by (Armitage 2010, Section 2.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 . We therefore find

With m_g being the molecular mass of H₂, 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.

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

and the radial momentum equation for the gas becomes

from which we obtain the circular velocity of gas:

The condition implies μ ⩽ 5^5/2, which is satisfied by protoplanetary disks. The bottom panel in Figure 1 illustrates the radial variation of Δv_c = (v_{p, c} − v_{g, c})/v_th for three values of q. The profile of Δv_c is almost flat for R ≲ 0.3. According to Equation (8), the gas streaming velocity exceeds the speed of particles for R² < 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 R² > 1/3 and inspiral toward the central star, they are accelerated when R² < 1/3 and migrate outward. Inward and outward migrating particles will then be accumulated near R² ≈ 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.

The dynamics of particles is described by the phase-space distribution function (DF) , where is the number density of particles. The vectors and (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),

to model the initial distribution of particles, before turning on the collisions between solid particles and gas molecules. Here, and 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 of the particle disk. The function g_K(E) takes physical (positive) values for K ⩾ 2, and 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)

where D[Δv_i] and D[Δv_iΔv_j] are diffusion coefficients and a_i 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 , where the perturbed potential V₁ is self-consistently calculated from the density disturbance Σ₁(R, ϕ, t) of particle distribution. Since we have assumed Σ_p ≪ Σ_g, the contribution of the gas component to V₁ 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 , where β = r_p + r_g. Let us define γ = m_g/(m_p + m_g) and denote the absolute velocity of gas molecules by . We obtain

where the potential functions and are given by the integrals

and 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 , where is the velocity component perpendicular to the disk plane. Defining

the elements of the stress tensor are determined as . The macroscopic quantities Σ_p, , and τ_ij are functions of 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 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, m_g/m_p ≪ 1, which implies . The Fokker–Planck equation can therefore be reduced to

with the dimensionless motion equations

Differentiating Equation (12) with respect to v_i gives

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

where is the orbital frequency of particles. We think of disks with v_th < v_{p, c}, v_{g, 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 . For , which corresponds to protoplanetary disks around solar-type stars, and for K ⩽ 29 experimented in JT12, one has

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 . The integral in Equation (17) is thus approximated by

where ρ_g = Σ_g(R)δ(z) and 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_⋆/b³ 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 , with ρ_s being the typical density of rocky material, one finds ξ₀∝1/r_p. With ξ₀ = (M_⋆/m_p)(r_p/b) the scattering of gas molecules takes place in the disk plane and the cross section becomes a line of the length r_p + r_g. We are not interested in this extreme unphysical case.

For near-circular orbits with e → 0, one obtains and D[Δv_i] transforms to the well-known form of Epstein drag:

During our numerical computations we use Equation (20) if e v_{p, c} > v_th 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 , and exactly performing the integral in Equation (17). Nonetheless, the missing factor is unimportant in our computations, for we will vary ξ₀ to explore the influence of drag force on the disk evolution, and one may suppose that any constant factors have already been included in ξ₀. To compute diffusion coefficients numerically, we soften Dirac's delta function using its normal distribution representation:

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

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 and their conjugate actions . We follow JT12 and set J₁ and J₂ to the radial action J_R and the orbital angular momentum J_ϕ = Rv_ϕ, respectively. In the -space, the motion equations (16) become

where is the Hamiltonian function, and the generalized forces and are determined using the virtual work of nonconservative forces:

with δ being the variational operator. In the Appendix, we explain the procedure of calculating and . They are real harmonic functions of w₁ and are smooth in the -space. We now write Equation (15) in the angle–action space:

where [ ·  , ·] denotes the Poisson bracket over the -space and the collision operator is defined by

We seek solutions of the form for Equation (27) so that |f₁| ≪ |f₀|. The Hamiltonian corresponding to f will become , where the perturbed potential V₁ (self-consistently arising from f₁) and are expressed in the angle–action space (Jalali 2010, Section 2). f₁ is obtained by solving the perturbed Fokker–Planck equation:

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 . For the small disturbances |f_d| ≪ |f₀|, we can ignore against and write

where the left-hand side is the linear approximation of the total derivative df_d/dt, and the potential V_d corresponds to f_d. From the definition of and Equations (A8)–(A10), we arrive at

whose terms lead to particular periodic solutions . However, the dominant collisional term

results in a secular drift f_d ∼ S₀t in the phase space: the DF of particles linearly increases in time if their actions satisfy , and it decreases for . Such drag-induced migrations can accumulate particles in regions where has a positive local maximum. Since S₀ (and therefore df_d/dt) depends on the initial DF, it is useful to normalize it to f₀ 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 S₀/(f₀ξ) in Figure 2 for mass parameters μ = 0.04 and μ_p = M_p/M_g = 1/8. The mean eccentricity of particle orbits of this model is . It is seen that S₀/(f₀ξ) 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 (a₀, e₀) ≈ (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 S₀/(f₀ξ) 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 . The mean semimajor axis of classical KBOs is . 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 , 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 a_AST ≈ 2.1 to 3.3 AU. The most populous group of these asteroids has a mean semimajor axis of . The second plausible length scale of our model is therefore 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 S₀/(f₀ξ), but proportionally decreases its maximum value at the accumulation point.

Zoom In Zoom Out Reset image size

4.2. Exponentially Growing Instabilities

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

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

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

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

and density Σ₁ 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 :

where is the precession frequency of particle orbits. The orbital frequencies (i = 1, 2) are determined in the unperturbed state. For nearly circular orbits one finds (JT12)

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 R_n = 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 A₁ (fundamental mode) and A₂ (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 J₁ = 0, and vanishes for radial orbits. Singular modes are associated with the inner Lindblad resonance 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 A₁ and A₂ versus ξ

Model	Model Parameters					Mode A1		Mode A2
	μ	μp	ϖmin × 103	ϖmax × 103	ξ	ωr/ϖmax	s × 105	ωr/ϖmax	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 A₁ and A₂ have been reported in Table 2 for 16 different models whose particle and gas phases are Toomre stable. It is seen that ω_r/ϖ_max is a function of μ_p, and s is an almost linear function of . 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 A₁ always rotates and grows faster than mode A₂. 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 M_p < M_d and the pattern speeds of the modes with the longest wavelengths have approached . 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).

Zoom In Zoom Out Reset image size

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 A₁ and A₂. 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%.

Zoom In Zoom Out Reset image size

For discrete slow modes, resonant cavities become smaller as μ_p falls off and mode A₂ disappears: mode A₂ has hardly managed to exist in model 13, which is stable. Nonexistence of mode A₂ 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 A₂ (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 A₂ 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 A₁ can marginally tolerate ω_r/ϖ_max ≲ 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 . The number of slow modes also depends on the mean eccentricity of particle orbits. Our numerical experiments with μ_p = 1/8 show that mode A₂ completely disappears as its pattern speed drops below ϖ_max by increasing .

Except for ξ = 0, 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 S₀/f₀ decides which process wins. According to Figures 2 and 3, near the major peaks of modes A₁ and A₂ we have S₀/f₀ ≈ 2 × 10⁻⁴ξ, 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 e^st can therefore overwhelm secular migrations, consume most solid particle reserves of the disk within R < a₀, 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.

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 (top panel in Figure 5). Interestingly, the semimajor axis of Uranus matches the location of the major density bump of mode A₁ and Neptune is clearly associated with the outer bump. If we assume that both planets were formed exactly at the peaks of mode A₁, 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.

Zoom In Zoom Out Reset image size

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 . Doing so, present orbital semimajor axes of Venus, the Earth, and Mars will lie in the region affected by the density bumps of modes A₁ and A₂ (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:

where is the epicyclic frequency of near-circular orbits. The function Q(R) has a global minimum at R_cr = 1.856, and its minimum value is . Near the orbit of Jupiter we have q = 0.024. The gas phase thus becomes unstable there if Q_min < 1, which implies . 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 , short-wavelength instabilities in the gas phase can be triggered around r_cr = R_crb = 7.07 AU, which is halfway between the orbits of Jupiter and Saturn. For , equation Q(R) − 1 = 0 has two roots at r₁ = 4.95 and r₂ = 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.

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⁻³ 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 . Moreover, the radial velocity dispersion of particles in our disks is determined from

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

where h_p is the scale height of the particle disk. In low-mass disks particles can be scattered up to a vertical distance of , and Equation (41) yields

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

with . Since particle–particle collisions can be ignored only for , 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 r_p ≫ 1.2 m. This size threshold reduces to r_p ≫ 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)

Note that ρ_g is the dimensionless spatial density of the gas disk (see Section 3). For molecular hydrogen, we have m_g = 3.32 × 10⁻²⁷ kg and the collision cross section is . At R = 0.3 and in the disk mid-plane we obtain

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 r_p 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:

where the drag coefficient C_D 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 , the drag coefficient approximately becomes C_D ≈ 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

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 A₁ 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, C_D = 0.92, and b = 3.81 AU, we require r_p ≳ 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 C_D = 2 (collisional/Epstein drag regime) and b = 71.43 AU, we obtain r_p ≳ 170 m, which is smaller than the mean free path of gas molecules in the outer solar system.

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

satisfies τ_s ≫ 1 in regions affected by unstable modes. The stopping time t_stop 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.

The infall timescale of 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

where α is the dimensionless turbulent diffusivity. For α ≲ 10⁻³ 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 f₀ and a small mass ratio μ_p, we investigated global modes only for K = 29 that gives . 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 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.

We define (v_R, v_ϕ) as the velocity components of particles in the polar coordinates (R, ϕ), and the streaming velocity of the gas component will become , where U_ϕ(R) has the functional form of v_{g, c}. Using the action variables and their conjugate angles , the radial distance and azimuthal position of a test particle are calculated from

where the Fourier coefficients are given by

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

The virtual work of the drag force reads

where C_D(R) = ξ₀ ev_{p, c} ρ_g if the orbital eccentricity e satisfies the inequality e > v_th/v_{p, c} and C_D(R) = ξ₀v_thρ_g otherwise. We now utilize the following Fourier expansions:

and obtain