A Thermodynamic View of Dusty Protoplanetary Disks

We consider the limit of short cooling times, ${t}_{\mathrm{cool}}\to 0$, appropriate for the outer parts of an irradiated protoplanetary disk (Chiang & Goldreich 1997; Lin & Youdin 2015). Then the disk temperature $T={T}_{\mathrm{ref}}$ at all times, and so we may adopt a locally isothermal equation of state

$$\begin{eqnarray}&&P={c}_{{\rm{s}}}^{2}(r,z){\rho }_{{\rm{g}}}={c}_{{\rm{s}}}^{2}(1-{f}_{{\rm{d}}})\rho ,\end{eqnarray} \tag{ 13 }$$

where ${c}_{{\rm{s}}}(r,z)=\sqrt{{ \mathcal R }{T}_{\mathrm{ref}}/\mu }$ is a prescribed sound-speed profile fixed in time. In most applications we consider vertically isothermal disks with

$$\begin{eqnarray}&&{c}_{{\rm{s}}}^{2}(r)\propto {r}^{q},\end{eqnarray} \tag{ 14 }$$

where q is the power-law index for the disk temperature. For q = 0 the disk is strictly isothermal.

Notice that Equation (13) resembles an ideal-gas equation of state but with a reduced temperature $\widetilde{T}={T}_{\mathrm{ref}}(1-{f}_{{\rm{d}}})$. This reduced temperature decreases with dust-loading. Since ${f}_{{\rm{d}}}$ typically decreases away from the midplane, we expect vertically isothermal dusty disks to behave as if the temperature increased away from z = 0.

Although we have deleted the true energy equation by fixing a locally isothermal equation of state, we show that the mixture nevertheless obeys an effective energy evolution equation. This is because advection of the dust fraction, described by Equation (6), can be transformed into an energy-like equation. The equation of state, Equation (13), implies

$$\begin{eqnarray*}&&{f}_{{\rm{d}}}=1-\displaystyle \frac{P}{{c}_{{\rm{s}}}^{2}(r,z)\rho }.\end{eqnarray*}$$

Then Equation (6) becomes

$$\begin{eqnarray}&&\displaystyle \frac{\partial P}{\partial t}+{\boldsymbol{v}}\cdot {\rm{\nabla }}P=-P{\rm{\nabla }}\cdot {\boldsymbol{v}}+P{\boldsymbol{v}}\cdot {\rm{\nabla }}\mathrm{ln}{c}_{{\rm{s}}}^{2}+{ \mathcal C },\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}&&{ \mathcal C }\equiv {c}_{{\rm{s}}}^{2}{\rm{\nabla }}\cdot \left[{t}_{{\rm{s}}}\left(1-\displaystyle \frac{P}{{c}_{{\rm{s}}}^{2}\rho }\right){\rm{\nabla }}P\right].\end{eqnarray} \tag{ 16 }$$

For comparison with the standard energy equation in hydrodynamics, rewriting Equation (8) for a pure ideal gas with $P\propto {\rho }_{{\rm{g}}}T$, and reverting back to gas velocities, gives

$$\begin{eqnarray}&&\displaystyle \frac{\partial P}{\partial t}+{{\boldsymbol{v}}}_{{\rm{g}}}\cdot {\rm{\nabla }}P=-\gamma P{\rm{\nabla }}\cdot {{\boldsymbol{v}}}_{{\rm{g}}}+\displaystyle \frac{P}{T}({ \mathcal H }-{\rm{\Lambda }}).\end{eqnarray} \tag{ 17 }$$

We can thus interpret Equation (15) as the energy equation for an ideal gas with adiabatic index $\gamma =1$, but now with the imposed temperature gradient ${\rm{\nabla }}{c}_{{\rm{s}}}^{2}$ and dust–gas drag ${ \mathcal C }$ acting as source terms.

If we denote the pressure forces by

$$\begin{eqnarray}&&{\boldsymbol{F}}\equiv -\displaystyle \frac{{\rm{\nabla }}P}{\rho },\end{eqnarray} \tag{ 18 }$$

then

$$\begin{eqnarray*}&&{ \mathcal C }=-{c}_{{\rm{s}}}^{2}{\rm{\nabla }}({f}_{{\rm{d}}}{t}_{{\rm{s}}}\rho {\boldsymbol{F}}),\end{eqnarray*}$$

which is in the same form as cooling by radiative diffusion. In protoplanetary disks the corresponding "heat flux," proportional to ${\boldsymbol{F}}$, is directly radially outward and vertically upward. This is simply a reflection of particle drift toward increasing pressure (inward and downward). Particle flux into a region contributes to "cooling" of that region because the reduced temperature is lowered.

The specific entropy of an ideal gas is $S={C}_{P}\mathrm{ln}({P}^{1/\gamma }/{\rho }_{{\rm{g}}})$, where C_P is the heat capacity at constant pressure.

Since we have shown that a locally isothermal dusty gas effectively has $\gamma =1$ we can define an effective entropy for the mixture as

$$\begin{eqnarray}&&{S}_{\mathrm{eff}}\equiv \mathrm{ln}\displaystyle \frac{P}{\rho }=\mathrm{ln}[{c}_{{\rm{s}}}^{2}(1-{f}_{{\rm{d}}})],\end{eqnarray} \tag{ 19 }$$

where the constant heat capacity has been absorbed into ${S}_{\mathrm{eff}}$. Then combining Equations (15) and (5) gives

$$\begin{eqnarray*}&&\displaystyle \frac{{{DS}}_{\mathrm{eff}}}{{Dt}}={\boldsymbol{v}}\cdot {\rm{\nabla }}\mathrm{ln}{c}_{{\rm{s}}}^{2}+\displaystyle \frac{{c}_{{\rm{s}}}^{2}}{P}{\rm{\nabla }}\cdot ({t}_{{\rm{s}}}{f}_{{\rm{d}}}{\rm{\nabla }}P),\end{eqnarray*}$$

which is equivalent to entropy evolution in standard hydrodynamics, albeit with source terms. With an effective entropy defined this way, many of the results concerning the stability of (locally) isothermal dusty gas will have identical form and interpretations to those for pure gas.

The physical reason for this analogy is that with strong drag (${t}_{{\rm{s}}}\to 0$), dust is almost perfectly entrained in the gas, but there is some gain/loss of dust particles between different gas parcels. This is analogous to the entropy of an ideal pure gas subject to heating/cooling: entropy is conserved following a gas parcel, except if there is heat exchange between a fluid parcel and its surroundings. For strictly isothermal gas perfectly coupled to dust (constant c_s² and ${t}_{{\rm{s}}}=0$) we have ${{DS}}_{\mathrm{eff}}/{Dt}={{Df}}_{{\rm{d}}}/{Dt}=0$. In this case the effective entropy is exactly conserved because the dust fraction is "frozen in" the flow.

We can now define the vertical buoyancy frequency N_z of the mixture as

$$\begin{eqnarray}&&{N}_{z}^{2}\equiv -\displaystyle \frac{1}{\rho }\displaystyle \frac{\partial P}{\partial z}\displaystyle \frac{\partial {S}_{\mathrm{eff}}}{\partial z}={c}_{{\rm{s}}}^{2}(r)\displaystyle \frac{\partial \mathrm{ln}{\rho }_{{\rm{g}}}}{\partial z}\displaystyle \frac{\partial {f}_{{\rm{d}}}}{\partial z},\end{eqnarray} \tag{ 20 }$$

where the second equality assumes a vertically isothermal disk. Protoplanetary disks have ${\partial }_{z}{\rho }_{{\rm{g}}}\lt 0$, thus stability against vertical convection (${N}_{z}^{2}\gt 0$) requires ${\partial }_{z}{f}_{{\rm{d}}}\lt 0$, i.e., the dust density should drop faster than the gas density away from the midplane. This is equivalent to entropy increasing away from the midplane. A similar expression exists for the radial buoyancy frequency N_r, but $| {N}_{r}| \ll | {N}_{z}|$ in thin, smooth disks.

A vertical buoyancy force exists even in vertically isothermal dusty disks, because coupling gas to dust particles increases the fluid's inertia, but pressure (i.e., restoring) forces are unchanged. The increased weight of the fluid resists vertical oscillations and hence there is an associated buoyancy force. However, if ${t}_{{\rm{s}}}\ne 0$ so the gas–dust coupling is imperfect, then gas is no longer "weighed down" by the dust and can slip past it. Thus finite drag diminishes the dust-induced buoyancy. This is similar to reducing gas buoyancy through cooling (Lin & Youdin 2015).

The above correspondence between isothermal dusty gas and adiabatic pure gas is derived under the terminal velocity approximation, which applies to short stopping times, ${t}_{{\rm{s}}}{{\rm{\Omega }}}_{{\rm{K}}}\ll 1$, and the locally isothermal approximation, which applies to short cooling times, ${t}_{\mathrm{cool}}{{\rm{\Omega }}}_{{\rm{K}}}\ll$ 1. In typical protoplanetary disk models such as the Minimum Mass Solar Nebula, ${t}_{{\rm{s}}}$ decreases for smaller particles and/or toward smaller radii (Youdin 2011). However, t_cool is generally small only in the outer disk (Lin & Youdin 2015; Malygin et al. 2017). Combining these results, we estimate that the correspondence will be applicable to particles of less than millimeter size at a few to tens of au in protoplanetary disks. However, note that the isothermal approximation may be relaxed to other fixed equations of state (see Section 7.2).

For a given distribution of the dust fraction ${f}_{{\rm{d}}}$ (or dust-to-gas ratio ), the mass and momentum Equations (5)–(7) admit solutions with $\rho (r,z)$ and ${\boldsymbol{v}}=r{\rm{\Omega }}(r,z)\hat{{\boldsymbol{\phi }}}$ where ${\rm{\Omega }}={v}_{\phi }/r$, which satisfy

$$\begin{eqnarray}&&r{{\rm{\Omega }}}^{2}=\displaystyle \frac{\partial {\rm{\Phi }}}{\partial r}+\displaystyle \frac{1}{\rho }\displaystyle \frac{\partial P}{\partial r},\end{eqnarray} \tag{ 21 }$$

$$\begin{eqnarray}&&0=\displaystyle \frac{\partial {\rm{\Phi }}}{\partial z}+\displaystyle \frac{1}{\rho }\displaystyle \frac{\partial P}{\partial z},\end{eqnarray} \tag{ 22 }$$

with $P=P({f}_{{\rm{d}}},\rho )$ given by the equation of state (Equation (13)). An explicit solution is presented in Section 5.1, when we analyze the stability of protoplanetary disks.

The mixture possesses vertical shear. To see this, we eliminate Φ between Equations (21) and (22) to obtain

$$\begin{eqnarray}&&r\displaystyle \frac{\partial {{\rm{\Omega }}}^{2}}{\partial z}=\displaystyle \frac{1}{\rho }\left(\displaystyle \frac{\partial P}{\partial r}\displaystyle \frac{\partial {S}_{\mathrm{eff}}}{\partial z}-\displaystyle \frac{\partial P}{\partial z}\displaystyle \frac{\partial {S}_{\mathrm{eff}}}{\partial r}\right).\end{eqnarray} \tag{ 23 }$$

and recall that ${S}_{\mathrm{eff}}=\mathrm{ln}[{c}_{{\rm{s}}}^{2}(1-{f}_{{\rm{d}}})]$. It is well appreciated that vertical stratification of the dust layer contributes to vertical shear (Chiang 2008). Equation (23) shows that a radial dust stratification (${\partial }_{r}{f}_{{\rm{d}}}$) also contributes to vertical shear.

Our equilibrium solutions satisfy the hydrostatic constraints of Equations (21)–(22) but are not in general steady-state solutions to the energy Equation (15), because the source term ${ \mathcal C }\ne 0$ for realistic disks. Limiting cases with ${ \mathcal C }\equiv 0$ include (1) unstratified or 2D, razor-thin disk models with finite dust–gas drag; (2) perfectly coupled dust with ${t}_{{\rm{s}}}=0$. The following analyses thus apply strictly to these limiting cases of ${ \mathcal C }=0$ in equilibrium. Our analysis is also a good approximation when ${t}_{{\rm{s}}}$ is sufficiently small such that the evolution of the background disk (e.g., dust settling) occurs on much longer timescales than the instability growth timescales.

We consider axisymmetric Eulerian perturbations to a variable X of the form

$$\begin{eqnarray}&&\mathrm{Re}[\delta X(r,z)\exp (-i\sigma t)],\end{eqnarray} \tag{ 24 }$$

where σ is the complex mode frequency. We write

$$\begin{eqnarray}&&\sigma ={is}-\omega ,\end{eqnarray} \tag{ 25 }$$

where s and ω are the growth rate and real frequencies, respectively. Then perturbations have time dependence ${e}^{{st}+i\omega t}$. Thus for $\omega \gt 0$, perturbations rotate anticlockwise in the complex plane.

In Appendix A we linearize Equations (5), (7), and (15) to derive the following integral relation:

$$\begin{eqnarray}{\sigma }^{2}{{ \mathcal I }}^{2} & = & \displaystyle \int [\rho | \delta {v}_{r}{| }^{2}A+\rho (\delta {v}_{z}\delta {v}_{r}^{* }+\delta {v}_{z}^{* }\delta {v}_{r})B+\rho | \delta {v}_{z}{| }^{2}D\\ & & \left.+\displaystyle \frac{1}{P}{| {\rm{\nabla }}\cdot (P\delta {\boldsymbol{v}})| }^{2}\right]{dV}-\displaystyle \int ({\rm{\nabla }}\cdot \delta {{\boldsymbol{v}}}^{* })\delta { \mathcal C }{dV}\\ & & -\displaystyle \int P({\rm{\nabla }}\cdot \delta {{\boldsymbol{v}}}^{* })(\delta {\boldsymbol{v}}\cdot {\rm{\nabla }}\mathrm{ln}{c}_{{\rm{s}}}^{2}){dV},\end{eqnarray} \tag{ 26 }$$

where * denotes the complex conjugate,

$$\begin{eqnarray}&&{{ \mathcal I }}^{2}\equiv \int \rho (| \delta {v}_{r}{| }^{2}+| \delta {v}_{z}{| }^{2}){dV}\end{eqnarray} \tag{ 27 }$$

is the meridional kinetic energy, and coefficients $A,B,D$ can be read off Equation (77). We now consider various limits of Equation (26).

When c_s² is a constant and ${t}_{{\rm{s}}}=0$ (so that ${ \mathcal C }=0$), the dusty gas equations are exactly equivalent to those for adiabatic hydrodynamics with unit adiabatic index. Although the gas is strictly isothermal, the mixture behaves adiabatically because the dust fraction ${f}_{{\rm{d}}}$ is advected with the gas. This is similar to entropy conservation following an adiabatic gas without heating or cooling.

In this case, the last two integrals in Equation (26) vanish. Then the condition for axisymmetric stability, ${\sigma }^{2}\gt 0$, is met if the integrand coefficients satisfy $A+D\gt 0$ and ${AD}-{B}^{2}\gt 0$ (Ogilvie 2016, Section 11.6), which yields

$$\begin{eqnarray}&&{\kappa }^{2}+\displaystyle \frac{1}{{\rho }_{{\rm{g}}}}{\rm{\nabla }}P\cdot {\rm{\nabla }}{f}_{{\rm{d}}}\gt 0,\end{eqnarray} \tag{ 28 }$$

$$\begin{eqnarray}&&-\displaystyle \frac{1}{{\rho }_{{\rm{g}}}}\displaystyle \frac{\partial P}{\partial z}\left(-{\kappa }^{2}\displaystyle \frac{\partial {f}_{{\rm{d}}}}{\partial z}+r\displaystyle \frac{\partial {{\rm{\Omega }}}^{2}}{\partial z}\displaystyle \frac{\partial {f}_{{\rm{d}}}}{\partial r}\right)\gt 0,\end{eqnarray} \tag{ 29 }$$

where ${\kappa }^{2}\equiv {r}^{-3}{\partial }_{r}({r}^{4}{{\rm{\Omega }}}^{2})\gt 0$ is the square of the epicyclic frequency. Note that ${\rm{\nabla }}P/{\rho }_{{\rm{g}}}={c}_{{\rm{s}}}^{2}{\rm{\nabla }}\mathrm{ln}{\rho }_{{\rm{g}}}$ for strictly isothermal gas considered here.

Equations (28)–(29) can in fact be obtained by inserting ${C}_{P}\mathrm{ln}(1-{f}_{{\rm{d}}})$ for the entropy into the standard expression for the Solberg–Hoiland criteria (Tassoul 1978) for axisymmetric stability of adiabatic gas. This substitution is consistent with our definition of the effective entropy ${S}_{\mathrm{eff}}$ since we are considering constant c_s.

We expect Equation (28) to be satisfied in typical protoplanetary disks where the dust-to-gas ratio increases in the same direction as the local pressure gradient, which is equivalent to the gas density gradient for isothermal gas.

On the other hand, Equation (29) can be violated in disks if the dust is vertically well-mixed but radially stratified, such that ${\partial }_{z}{f}_{{\rm{d}}}=0$ but ${\partial }_{r}{f}_{{\rm{d}}}\ne 0$. In this case the left of Equation (29) becomes

$$\begin{eqnarray}&&-{\left(\displaystyle \frac{1}{{\rho }_{{\rm{g}}}}\displaystyle \frac{\partial P}{\partial z}\right)}^{2}{\left(\displaystyle \frac{\partial {f}_{{\rm{d}}}}{\partial r}\right)}^{2}\lt 0.\end{eqnarray} \tag{ 30 }$$

Such an isothermal dusty disk is unstable because there is no effective vertical buoyancy to stabilize vertical motions (N_z = 0), which can tap into the free energy associated with vertical shear due to the radial gradient in the dust fraction. This situation is identical to the adiabatic simulations of Nelson et al. (2013) for pure gas, where the gas entropy is vertically uniform but has a radial gradient. The authors indeed find instability. We give a numerical example in the dusty context in Section 6.6.

The physical interpretation of Section 2 is here derived in detail. Consider constant c_s but ${t}_{{\rm{s}}}\ne 0$ (so $\delta { \mathcal C }\ne 0$). Equation (26) indicates that ${\sigma }^{2}$ is generally complex, unless the second integral on the right is real. Thus we may have growing oscillations, or overstability, due to dust–gas friction. This is seen by taking the imaginary part of Equation (26):

$$\begin{eqnarray}&&s=\displaystyle \frac{\mathrm{Im}\int ({\rm{\nabla }}\cdot \delta {{\boldsymbol{v}}}^{* })\delta { \mathcal C }{dV}}{2\omega {{ \mathcal I }}^{2}},\end{eqnarray} \tag{ 31 }$$

assuming $\omega \ne 0$. Since drag ($\delta { \mathcal C }$) appears as a source term in our effective energy equation, the quantity $\mathrm{Im}[({\rm{\nabla }}\cdot \delta {{\boldsymbol{v}}}^{* })\delta { \mathcal C }]$ represents correlations between compression/expansion and heating/cooling.

It is well known that such correlations may lead to pulsational instabilities in stars (Cox 1967). We thus interpret dust-drag overstabilities in a similar way, adapting from the treatment of stellar pulsations by Cox (1967) and lecture notes by Samadi et al. (2015) and J. Christensen-Dalsgaard.⁴

4.4.1. Work Done by Dusty Gas

The physical interpretation of Equation (31) is that work done by pressure forces in the dusty gas leads to growth ($s\gt 0$) or decay ($s\lt 0$) in oscillation amplitudes. To demonstrate this, we calculate the average work done assuming periodic oscillations, and show that if the average work done is positive, then the oscillation amplitude would actually grow.

Consider oscillations in the dusty gas with period T_p. The average rate of work done is

$$\begin{eqnarray}&&{ \mathcal W }=\displaystyle \frac{1}{{T}_{p}}{\int }_{t}^{t+{T}_{p}}{{dt}}^{{\prime} }{\int }_{M}P\displaystyle \frac{D\upsilon }{{{Dt}}^{{\prime} }}{dm}\end{eqnarray} \tag{ 32 }$$

(Cox 1967, see their Equation (4.10) and related discussions), where $\upsilon =1/\rho$ is the specific volume of the mixture, D/Dt is the Lagrangian derivative, the spatial integral is taken over the total mass M of the mixture and dm is a mass element.

Here we consider Lagrangian perturbations such that

$$\begin{eqnarray}&&P\to P+\mathrm{Re}({\rm{\Delta }}{{Pe}}^{-i\sigma t})\end{eqnarray} \tag{ 33 }$$

is the pressure following a fluid element of the mixture (and similarly for ρ). The Lagrangian perturbation Δ of a variable X is ${\rm{\Delta }}X=\delta X+{\boldsymbol{\xi }}\cdot {\rm{\nabla }}X$ and ${\boldsymbol{\xi }}$ is the Lagrangian displacement, so ${\xi }_{x,z}=i\delta {v}_{x,z}/\sigma$. Inserting the above pressure and density fields into Equation (32), and noting that only products of perturbations contribute to ${ \mathcal W }$ after time-averaging, we find for periodic oscillations (time dependence ${e}^{i\omega t}$ and real ω) that

$$\begin{eqnarray}&&{ \mathcal W }=-\displaystyle \frac{\omega }{2}\int \mathrm{Im}\left({\rm{\Delta }}P\displaystyle \frac{{\rm{\Delta }}{\rho }^{* }}{\rho }\right){dV},\end{eqnarray} \tag{ 34 }$$

Equation (34) shows that a phase difference between gas pressure and the total density leads to work done (${ \mathcal W }\ne 0$).

Now, Equations (67) and (68) imply that the integrand of the numerator in Equation (31) is

$$\begin{eqnarray}&&\mathrm{Im}(\delta { \mathcal C }{\rm{\nabla }}\cdot \delta {{\boldsymbol{v}}}^{* })=-| \sigma {| }^{2}\mathrm{Im}\left({\rm{\Delta }}P\displaystyle \frac{{\rm{\Delta }}{\rho }^{* }}{\rho }\right),\end{eqnarray} \tag{ 35 }$$

for constant c_s. Then combining Equations (34), (35), and (31) gives

$$\begin{eqnarray}&&s=\displaystyle \frac{{ \mathcal W }}{{{ \mathcal I }}^{2}},\end{eqnarray} \tag{ 36 }$$

where we have set $\sigma =-\omega$ in Equation (35). Equation (36) states that if the average work done is positive during an oscillation, ${ \mathcal W }\gt 0$, then its amplitude will actually grow ($s\gt 0$).

Positive work is done by a fluid parcel if $-\omega \mathrm{Im}({\rm{\Delta }}P{\rm{\Delta }}{\rho }^{* })\gt 0$. Without loss of generality, take $\omega \gt 0$ and consider a mass element with ${\rm{\Delta }}\rho =1$. Then positive work requires $\mathrm{Im}({\rm{\Delta }}P)\lt 0$. This corresponds to Lagrangian pressure perturbations lagging behind those in density; see Figure 3 for the case of dusty gas.

Zoom In Zoom Out Reset image size

4.4.2. Physical Properties of Overstabilities in Dust–Gas Drag

The above discussion applies to any single fluid with a pressure and density. In the case of interest—dusty gas—the work done is attributed to finite dust–gas drag. The relative drift between gas and dust, which only exists if ${t}_{{\rm{s}}}{\rm{\nabla }}P\ne 0$ (Equation (4)), causes a phase difference between the two components, and hence between pressure and total density.

A parcel of the strictly isothermal dusty gas mixture does positive work if

$$\begin{eqnarray*}&&-\mathrm{sgn}(\omega )\mathrm{Im}({\rm{\Delta }}{\rho }_{{\rm{g}}}{\rm{\Delta }}{\rho }_{{\rm{d}}}^{* })\gt 0,\end{eqnarray*}$$

meaning that gas follows dust (Figure 3). Overstabilities are thus not possible if the gas does not respond to dust (i.e., no back-reaction). The pressure–density lag shown in Figure 3 is achieved if, just after the total density of a parcel maximizes, its gas content is increasing, which requires a sufficiently large particle flux out of the parcel. See A to B in Figure 2.

This thermodynamic interpretation does not explain why drag forces cause gas pressure to lag behind the dust density, but it shows that this must be the case for any growing oscillations associated with the dust–gas drag. To rigorously understand how dust drag causes this lag requires an explicit solution to the linearized equations with detailed treatment of the function ${ \mathcal C }$. However, given the complexity of ${ \mathcal C }$ (see Appendix C), we might generally expect dusty disks to support a range of stable and overstable modes, with the latter being associated with pressure–density lag. In Section 5.3.2 we check that the SI fits into this thermodynamic interpretation in the limit of strong drag.

If ${c}_{{\rm{s}}}(r,z)$ is nonuniform but ${t}_{{\rm{s}}}=0$, Equation (26) gives

$$\begin{eqnarray}&&s=\displaystyle \frac{\mathrm{Im}\int P({\rm{\nabla }}\cdot \delta {{\boldsymbol{v}}}^{* })(\delta {\boldsymbol{v}}\cdot {\rm{\nabla }}\mathrm{ln}{c}_{{\rm{s}}}^{2}){dV}}{2\omega {{ \mathcal I }}^{2}},\end{eqnarray} \tag{ 37 }$$

again assuming $\omega \ne 0$. This instability represents VSI caused by vertical shear arising from a radial temperature gradient (Nelson et al. 2013; Barker & Latter 2015; Lin & Youdin 2015). We present numerical solutions of the VSI with perfectly coupled dust in Section 6.

We assume a Gaussian profile in the dust-to-gas ratio,

$$\begin{eqnarray}&&\epsilon (r,z)={\epsilon }_{0}(r)\exp \left[-\displaystyle \frac{{z}^{2}}{2{H}_{\epsilon }^{2}(r)}\right].\end{eqnarray} \tag{ 38 }$$

Inserting Equation (38) into vertical hydrostatic equilibrium, Equation (22), and integrating with the approximate gravitational potential (Equation (10)), we obtain the gas density as

$$\begin{eqnarray}&&{\rho }_{{\rm{g}}}(r,z)\\ &&\quad {=\rho }_{{\rm{g}}0}(r)\exp \left\{-\displaystyle \frac{{z}^{2}}{2{H}_{{\rm{g}}}^{2}}-{\epsilon }_{0}\displaystyle \frac{{H}_{\epsilon }^{2}}{{H}_{{\rm{g}}}^{2}}\left[1-\exp \left(-\displaystyle \frac{{z}^{2}}{2{H}_{\epsilon }^{2}}\right)\right]\right\},\end{eqnarray} \tag{ 39 }$$

where

$$\begin{eqnarray}&&{H}_{{\rm{g}}}=\displaystyle \frac{{c}_{{\rm{s}}}}{{{\rm{\Omega }}}_{{\rm{K}}}},\quad {{\rm{\Omega }}}_{{\rm{K}}}\equiv \sqrt{\displaystyle \frac{{{GM}}_{* }}{{r}^{3}}}.\end{eqnarray} \tag{ 40 }$$

Here H_g is the gas scale-height in the dust-free limit and ${{\rm{\Omega }}}_{{\rm{K}}}$ is the Keplerian frequency.

In gas-dominated disks with ${\epsilon }_{0}\lt 1$ the gas distribution ${\rho }_{{\rm{g}}}(r,z)$ is close to Gaussian, as in the dust-free case, and the dust density is approximately

$$\begin{eqnarray}&&{\rho }_{{\rm{d}}}\simeq {\epsilon }_{0}{\rho }_{{\rm{g}}0}(r)\exp \left(-\displaystyle \frac{{z}^{2}}{2{H}_{{\rm{d}}}^{2}}\right),\end{eqnarray} \tag{ 41 }$$

with

$$\begin{eqnarray}&&\displaystyle \frac{1}{{H}_{{\rm{d}}}^{2}}=\displaystyle \frac{1}{{H}_{\epsilon }^{2}}+\displaystyle \frac{1}{{H}_{{\rm{g}}}^{2}},\end{eqnarray} \tag{ 42 }$$

and H_d is the dust scale-height.

Finally, we define

$$\begin{eqnarray}&&Z\equiv {\epsilon }_{0}\displaystyle \frac{{H}_{{\rm{d}}}}{{H}_{{\rm{g}}}}\simeq \displaystyle \frac{{{\rm{\Sigma }}}_{{\rm{d}}}}{{{\rm{\Sigma }}}_{{\rm{g}}}}\end{eqnarray} \tag{ 43 }$$

as a measure of the local metallicity, where ${{\rm{\Sigma }}}_{{\rm{d}}}$ and ${{\rm{\Sigma }}}_{{\rm{g}}}$ are the dust and gas surface densities, respectively. The second equality holds for ${\epsilon }_{0}\ll 1$.

5.1.1. Orbital Frequency

From Equation (21) the disk orbital frequency is

$$\begin{eqnarray}&&{\rm{\Omega }}(r,z)={{\rm{\Omega }}}_{{\rm{K}}}(r){\left[1-\displaystyle \frac{3}{2}\displaystyle \frac{{z}^{2}}{{r}^{2}}+\displaystyle \frac{{h}_{{\rm{g}}}^{2}}{(1+\epsilon )}\displaystyle \frac{\partial }{\partial \mathrm{ln}r}\mathrm{ln}({c}_{{\rm{s}}}^{2}{\rho }_{{\rm{g}}})\right]}^{1/2},\end{eqnarray} \tag{ 44 }$$

where

$$\begin{eqnarray}&&{h}_{{\rm{g}}}\equiv \displaystyle \frac{{H}_{{\rm{g}}}}{r}\end{eqnarray} \tag{ 45 }$$

is the characteristic disk aspect ratio, with ${h}_{{\rm{g}}}\ll 1$ for protoplanetary disks (PPDs).

5.1.2. Vertical Shear

Writing Equation (23) in terms of the gas density and dust-to-gas ratio with a power-law temperature profile (Equation (14)) gives the disk's explicit vertical shear profile:

$$\begin{eqnarray}r\displaystyle \frac{\partial {{\rm{\Omega }}}^{2}}{\partial z} & = & \displaystyle \frac{{c}_{{\rm{s}}}^{2}(r)}{{(1+\epsilon )}^{2}}\left\{\displaystyle \frac{\partial \epsilon }{\partial r}\displaystyle \frac{\partial \mathrm{ln}{\rho }_{{\rm{g}}}}{\partial z}-\displaystyle \frac{\partial \epsilon }{\partial z}\displaystyle \frac{\partial \mathrm{ln}P}{\partial r}\right.\\ & & \left.-\displaystyle \frac{q}{r}(1+\epsilon )\displaystyle \frac{\partial \mathrm{ln}{\rho }_{{\rm{g}}}}{\partial z}\right\}.\end{eqnarray} \tag{ 46 }$$

The first two terms correspond to vertical shear caused by spatial variations in the dust-to-gas ratio. The third term proportional to q corresponds to vertical shear due to the radial temperature gradient, which survives in the dust-free limit.

We can compare these sources by evaluating them using the equilibrium solutions in Section 5.1. We assume the disk is radially smooth so that ${\partial }_{r}\sim 1/r$. This gives

$$\begin{eqnarray}&&\displaystyle \frac{{| r{\partial }_{z}{\rm{\Omega }}| }_{\mathrm{dust}/\mathrm{gas}\mathrm{gradient}}}{{| r{\partial }_{z}{\rm{\Omega }}| }_{\mathrm{temp}.\mathrm{gradient}}}\sim \epsilon \displaystyle \frac{\max ({\delta }^{2},1)}{| q| (1+\epsilon ){\delta }^{2}},\end{eqnarray} \tag{ 47 }$$

where $\delta \equiv {H}_{\epsilon }/{H}_{{\rm{g}}}$. Since $| q| =O(1)$ in PPDs, Equation (47) indicates that vertical shear due to variations in the dust-to-gas ratio dominates over that due to the radial temperature gradient for thin dust layers such that ${\delta }^{2}\ll \epsilon$.

5.1.3. Dusty Vertical Buoyancy

For the above equilibrium the vertical buoyancy frequency is given explicitly as

$$\begin{eqnarray*}&&{N}_{z}^{2}=\displaystyle \frac{\epsilon }{{(1+\epsilon )}^{2}}{\left(\displaystyle \frac{z}{{H}_{\epsilon }}\right)}^{2}{{\rm{\Omega }}}_{{\rm{K}}}^{2},\end{eqnarray*}$$

where we have used ${c}_{{\rm{s}}}={c}_{{\rm{s}}}(r)$. Then

$$\begin{eqnarray*}&&{N}_{z}\lesssim O(\sqrt{{\epsilon }_{0}}{{\rm{\Omega }}}_{{\rm{K}}}).\end{eqnarray*}$$

However, for well-mixed dust layers such that ${H}_{\epsilon }\gg {H}_{{\rm{g}}}$, $\max ({N}_{z})$ may occur outside a finite vertical domain.

Here we apply the dusty analog of the Solberg–Hoiland criteria derived in Section 4.3 assuming strictly isothermal gas. As discussed there, the first criterion is generally satisfied. Thus we only consider the second condition, Equation (29).

We assume that the disk is approximately Keplerian so that $\kappa \simeq {{\rm{\Omega }}}_{{\rm{K}}}$, and that the vertical gas distribution is Gaussian. In terms of the dust-to-gas ratio , Equation (29) becomes

$$\begin{eqnarray}&&1-\displaystyle \frac{\epsilon {h}_{{\rm{g}}}^{2}}{{(1+\epsilon )}^{2}}\displaystyle \frac{\partial \mathrm{ln}\epsilon }{\partial \mathrm{ln}r}\left(-\displaystyle \frac{\partial \mathrm{ln}P}{\partial \mathrm{ln}r}+\displaystyle \frac{{H}_{\epsilon }^{2}}{{H}_{{\rm{g}}}^{2}}\displaystyle \frac{\partial \mathrm{ln}\epsilon }{\partial \mathrm{ln}r}\right)\gt 0,\end{eqnarray} \tag{ 48 }$$

for stability. The first term in brackets is usually stabilizing in PPDs for reasons given in Section 4.3. The second term in brackets is always destabilizing, but is small in smooth, thin disks with radial gradients $O(1/r)$, ${h}_{{\rm{g}}}\ll 1$, provided that ${H}_{\epsilon }/{H}_{{\rm{g}}}$ is not large (e.g., some dust settling has occurred). This means that typical dusty PPDs are stable to axisymmetric perturbations, even for arbitrarily thin dust layers.

To violate Equation (48) and obtain instability, notice that the left side is a quadratic in ${\partial }_{r}\epsilon$. Thus instability is possible for sufficiently large (in magnitude) radial gradients in the dust-to-gas ratio,

$$\begin{eqnarray}&&\displaystyle \frac{\partial \mathrm{ln}\epsilon }{\partial \mathrm{ln}r}\gt {\chi }_{+}\quad \mathrm{or}\quad \displaystyle \frac{\partial \mathrm{ln}\epsilon }{\partial \mathrm{ln}r}\lt {\chi }_{-},\end{eqnarray} \tag{ 49 }$$

for instability, where

$$\begin{eqnarray}&&{\chi }_{\pm }=\displaystyle \frac{1}{2}\displaystyle \frac{{H}_{{\rm{g}}}^{2}}{{H}_{\epsilon }^{2}}\left[\displaystyle \frac{\partial \mathrm{ln}P}{\partial \mathrm{ln}r}\pm \sqrt{{\left(\displaystyle \frac{\partial \mathrm{ln}P}{\partial \mathrm{ln}r}\right)}^{2}+4\displaystyle \frac{{H}_{\epsilon }^{2}}{{H}_{{\rm{g}}}^{2}}\displaystyle \frac{{(1+\epsilon )}^{2}}{\epsilon {h}_{{\rm{g}}}^{2}}}\,\right].\end{eqnarray} \tag{ 50 }$$

In typical accretion disks where ${\partial }_{r}P\lt 0$, instability is easier for increasing dust-to-gas ratios (${\partial }_{r}\epsilon \gt 0$).

We can neglect pressure gradients in Equation (50) if $r{\partial }_{r}\mathrm{ln}P\sim O(1)$ and ${H}_{\epsilon }/{H}_{{\rm{g}}}\gg \sqrt{\epsilon }{h}_{{\rm{g}}}/2(1+\epsilon )$. For example, if $\epsilon \simeq 0.01$ and ${h}_{{\rm{g}}}\simeq 0.05$, then we require ${H}_{\epsilon }/{H}_{{\rm{g}}}\gg 2\times {10}^{-3}$. This can be met if the dust is not well settled (e.g., due to a small amount of turbulence). Then instability requires

$$\begin{eqnarray}&&\left|\displaystyle \frac{\partial \mathrm{ln}\epsilon }{\partial r}\right|\gtrsim \displaystyle \frac{1}{{H}_{\epsilon }}\displaystyle \frac{(1+\epsilon )}{\sqrt{\epsilon }}.\end{eqnarray} \tag{ 51 }$$

That is, if the radial lengthscale of the dust-to-gas ratio is much less its vertical lengthscale, ${L}_{\epsilon }\ll O({H}_{\epsilon })$, then the system is unstable. Taking $\epsilon \sim 0.01$, we find that for thin dust layers with ${H}_{\epsilon }\simeq {H}_{{\rm{d}}}\ll {H}_{{\rm{g}}}$, instability requires ${L}_{\epsilon }\ll O(0.1{H}_{{\rm{g}}})$, i.e., the dust-to-gas ratio must vary on an extremely short radial lengthscale. This might be achieved, for example, at sharp edges associated with gaps opened by giant planets.

Technically, the above discussion is applicable only when c_s is constant and ${t}_{{\rm{s}}}=0$ (see Section 4.3). However, since dusty edges translate to sharp entropy gradients (Section 3.3), we may generally expect sharp features in the dust distribution of protoplanetary disks to be unstable.

We now specialize further and compute explicit solutions to the linear problem. We consider radially localized axisymmetric disturbances of the form

$$\begin{eqnarray}&&\delta X(r,z)=\delta {X}_{1}(r,z)\exp ({{ik}}_{x}r),\end{eqnarray} \tag{ 52 }$$

where k_x is a real wavenumber such that $| {k}_{x}r| \gg 1$, and the amplitude $\delta {X}_{1}(r,z)$ is a slowly varying function of r. Then ${\partial }_{r}\to {{ik}}_{x}$ when acting on the above primitive perturbations, and we may neglect curvature terms. We take ${k}_{x}\gt 0$ without loss of generality. Hereafter, we drop the subscript 1 on the amplitudes.

Introducing

$$\begin{eqnarray}&&W\equiv \displaystyle \frac{\delta \rho }{\rho },\quad Q\equiv \displaystyle \frac{\delta P}{\rho },\end{eqnarray} \tag{ 53 }$$

the linearized equations for vertically isothermal dusty gas with the pressure equation in place of the dust fraction (Equations (5), (7), (11), (15)) are then

$$\begin{eqnarray}&&i\sigma W={{ik}}_{x}\delta {v}_{r}+\delta {v}_{z}^{{\prime} }+\delta {v}_{r}{\partial }_{r}\mathrm{ln}\rho +\delta {v}_{z}{\partial }_{z}\mathrm{ln}\rho ,\end{eqnarray} \tag{ 54 }$$

$$\begin{eqnarray}&&-i\sigma \delta {v}_{r}=2{\rm{\Omega }}\delta {v}_{\phi }-{{WF}}_{r}-{{ik}}_{x}Q-{{ik}}_{x}\delta \psi ,\end{eqnarray} \tag{ 55 }$$

$$\begin{eqnarray}&&i\sigma \delta {v}_{\phi }=\displaystyle \frac{{\kappa }^{2}}{2{\rm{\Omega }}}\delta {v}_{r}+\displaystyle \frac{\partial {v}_{\phi }}{\partial z}\delta {v}_{z},\end{eqnarray} \tag{ 56 }$$

$$\begin{eqnarray}&&-i\sigma \delta {v}_{z}=-{{WF}}_{z}-[{Q}^{{\prime} }+Q{(\mathrm{ln}\rho )}^{{\prime} }]-\delta {\psi }^{{\prime} },\end{eqnarray} \tag{ 57 }$$

$$\begin{eqnarray}i\sigma Q & = & \displaystyle \frac{P}{\rho }({{ik}}_{x}\delta {v}_{r}+\delta {v}_{z}^{{\prime} })+\displaystyle \frac{1}{\rho }(\delta {v}_{r}{\partial }_{r}P+\delta {v}_{z}{\partial }_{z}P)\\ & & -\displaystyle \frac{P}{\rho }\delta {v}_{r}{\partial }_{r}\mathrm{ln}{c}_{{\rm{s}}}^{2}-\displaystyle \frac{\delta { \mathcal C }}{\rho },\end{eqnarray} \tag{ 58 }$$

$$\begin{eqnarray}&&\delta {\psi }^{{\prime\prime} }=4\pi G\rho W+{k}_{x}^{2}\delta \psi ,\end{eqnarray} \tag{ 59 }$$

where ${}^{{\prime} }\equiv {\partial }_{z}$ and recall ${\boldsymbol{F}}\equiv -{\rm{\nabla }}P/\rho$. We have temporarily restored self-gravity to discuss SGI in the next section. The linearized dust diffusion function $\delta { \mathcal C }$ is given in Appendix C.

Equations (54)–(59) are a set of ordinary differential equations in z. All coefficients and amplitudes are evaluated at a fiducial radius $r={r}_{0}$, but their full z-dependence is retained. We next discuss solutions to these equations. We first show that the above equations yield the SGI and SI in the strong drag regime in Sections 5.3.1 and 5.3.2, respectively. We then consider 3D, stratified disks in Section 6 to study how the addition of dust affects the VSI.

5.3.1. Secular Gravitational Instability

Consider a razor-thin, self-gravitating disk so that $\rho ={\rm{\Sigma }}\delta (z)$, where δ is the Dirac delta function and Σ is the total surface density. Here $\epsilon ={{\rm{\Sigma }}}_{{\rm{d}}}/{{\rm{\Sigma }}}_{{\rm{g}}}$ and ${f}_{{\rm{d}}}={{\rm{\Sigma }}}_{{\rm{d}}}/{\rm{\Sigma }}$, where ${{\rm{\Sigma }}}_{{\rm{d}},{\rm{g}}}$ are the dust and gas surface densities. The background disk is uniform and we neglect the vertical dimension in Equations (54)–(58). The linearized dust–gas drag term is then $-\delta { \mathcal C }/\rho ={t}_{{\rm{s}}}{f}_{{\rm{d}}}{c}_{{\rm{s}}}^{2}{k}_{x}^{2}Q$. The thin-disk solution to Equation (59) is $\delta \psi (z=0)=-2\pi G{\rm{\Sigma }}W/| {k}_{x}|$.

These simplifications yield the dispersion relation

$$\begin{eqnarray*}&&(i\sigma -{t}_{{\rm{s}}}{f}_{{\rm{d}}}{c}_{{\rm{s}}}^{2}{k}_{x}^{2})(2\pi G{\rm{\Sigma }}| {k}_{x}| -{\kappa }^{2}+{\sigma }^{2})=i\sigma {c}_{{\rm{s}}}^{2}{k}_{x}^{2}(1-{f}_{{\rm{d}}}).\end{eqnarray*}$$

Searching for slowly and purely growing modes, $\sigma ={is}$ with $| s| \ll \kappa$, we find

$$\begin{eqnarray}&&s=\displaystyle \frac{{t}_{{\rm{s}}}{f}_{{\rm{d}}}{c}_{{\rm{s}}}^{2}{k}_{x}^{2}(2\pi {\rm{\Sigma }}G| {k}_{x}| -{\kappa }^{2})}{{\kappa }^{2}-2\pi {\rm{\Sigma }}G| {k}_{x}| +{c}_{{\rm{s}}}^{2}{k}_{x}^{2}(1-{f}_{{\rm{d}}})}.\end{eqnarray} \tag{ 60 }$$

This is SGI mediated by strong dust–gas drag with negligible turbulent dust diffusion (Takahashi & Inutsuka 2014, their Equation (13) becomes our Equation (60) in this limit after a change of variables). A similar effect occurs in viscous self-gravitating gas disks (Gammie 1996; Lin & Kratter 2016). In fact, if we identify $\nu \equiv {t}_{{\rm{s}}}{f}_{{\rm{d}}}{c}_{{\rm{s}}}^{2}$ as a kinematic viscosity, then Equation (60) is identical to Gammie's Equation (18).

This exercise shows that the one-fluid framework, further simplified by the terminal velocity approximation, is sufficient to capture SGI in the limit of strong drag.

5.3.2. Streaming Instability

We now consider 3D disks without self-gravity. We neglect the vertical component of the stellar gravity, appropriate for studying regions near the disk midplane. This allows us to Fourier analyze in z to obtain an algebraic dispersion relation of the form ${\sum }_{j=0}^{5}{c}_{j}({k}_{x},{k}_{z}){\sigma }^{j}=0$, where k_z is a real vertical wavenumber. The coefficients c_j can be read off Equation (96) in Appendix D. There we also show that this dispersion relation reduces to that for the SI in the limit of incompressible gas and small ${t}_{{\rm{s}}}$ (Youdin & Goodman 2005; Jacquet et al. 2011).

In Table 1 we solve the full dispersion relation (Equation (96)) numerically for selected cases where analytic SI growth rates have been verified with particle–gas numerical simulations (from Youdin & Johansen 2007; Bai & Stone 2010b). Following previous works on SI, we use normalized wavenumbers ${K}_{x,z}=\eta {{rk}}_{x,z}$ where

$$\begin{eqnarray}&&\eta \equiv -\displaystyle \frac{1}{2{\rho }_{{\rm{g}}}r{{\rm{\Omega }}}_{{\rm{K}}}^{2}}\displaystyle \frac{\partial P}{\partial r}=\displaystyle \frac{1}{2(1-{f}_{{\rm{d}}})}\displaystyle \frac{{F}_{r}}{r{{\rm{\Omega }}}_{{\rm{K}}}^{2}}\end{eqnarray} \tag{ 61 }$$

measures the pressure offset of Keplerian rotation. We fix $\eta =0.05{c}_{{\rm{s}}}/r{{\rm{\Omega }}}_{{\rm{K}}}$. In this section we also quote the particle stopping time ${\tau }_{{\rm{s}}}={t}_{{\rm{s}}}/(1-{f}_{{\rm{d}}})$.

Table 1.  Selected Modes of the Streaming Instability

Mode (${\tau }_{{\rm{s}}}{{\rm{\Omega }}}_{{\rm{K}}}$, ${K}_{x,z},{\rho }_{{\rm{g}}}/{\rho }_{{\rm{d}}}$)	Complex Frequency, $\sigma /{{\rm{\Omega }}}_{{\rm{K}}}$		Work Done, ${ \mathcal W }/| {\rm{\Delta }}P{\rm{\Delta }}{\rho }^{* }/\rho |$		Pressure–density lag, φ
	Two-fluid	One-fluid	Two-fluid	One-fluid	Two-fluid	One-fluid
linA (1) ($0.1,30,3$)	0.3480 + 0.4190i	0.3640 + 0.4249i	0.078	0.090	$27^\circ$	$30^\circ$
linB (1) ($0.1,6,0.2$)	−0.4999 + 0.0155i	−0.4981 + 0.0054i	0.025	0.0054	$5\buildrel{\circ}\over{.} 8$	$1\buildrel{\circ}\over{.} 2$
linC (2) (${10}^{-2},1500,2$)	0.1049 + 0.5981i	0.1338 + 0.6650i	0.0076	0.013	$8\buildrel{\circ}\over{.} 3$	$11^\circ$
linD (2) (${10}^{-3},2000,2$)	0.3225 + 0.3154i	0.3219 + 0.3154i	0.061	0.061	$22^\circ$	$22^\circ$

References. (1) Youdin & Johansen (2007); (2) Bai & Stone (2010b).

Download table as:  ASCIITypeset image

The eigenfrequencies obtained from the one-fluid dispersion relation are compared with those from a full, two-fluid analysis (similar to Youdin & Goodman 2005; Kowalik et al. 2013). As expected, eigenfrequencies agree better with decreasing ${\tau }_{{\rm{s}}}$ since in that limit the mixture behaves more like a single fluid. Most importantly, we find the work done ${ \mathcal W }\gt 0$ in all cases, and hence find growing oscillations.

In Table 1 we also calculate the phase difference between the Lagrangian pressure perturbation and density perturbations as

$$\begin{eqnarray*}&&\varphi \equiv \mathrm{sgn}[\mathrm{Re}(\sigma )]\arg ({\rm{\Delta }}P{\rm{\Delta }}{\rho }^{* }).\end{eqnarray*}$$

(Note that it is important to include the global radial pressure gradient in ${\rm{\Delta }}P=\delta P+i\delta {v}_{x}{\partial }_{r}P/\sigma$.) Then $\varphi \gt 0$ indicates gas pressure lagging behind total density, which is true for all the cases. Thus SI is indeed associated with such a phase lag.

In Figure 4 we calculate the most unstable SI mode as a function of ${\tau }_{{\rm{s}}}$ at fixed K_z = 30 and $\epsilon =3$. Growth rates are maximized over K_x. We compare results between the full, two-fluid linear analysis and the one-fluid framework. We also include an analytic model, developed in Appendix D.2, based on the one-fluid dispersion relation with additional approximations (orange diamonds).

Zoom In Zoom Out Reset image size

Both results based on the one-fluid framework compare well with the full, two-fluid analysis, only breaking down at a relatively large ${\tau }_{{\rm{s}}}{{\rm{\Omega }}}_{{\rm{K}}}\gtrsim 0.1$. For larger ${\tau }_{{\rm{s}}}$ the two-fluid phase lag drops, along with the growth rate. This suggest that a non-vanishing phase lag is indeed necessary for instability. However, the magnitude of the phase lag does not correlate with growth rates. In fact, φ remains finite as $| \sigma | ,{\tau }_{{\rm{s}}}\to 0$ (but non-zero). This arises because the optimum ${K}_{x}\propto {\tau }_{{\rm{s}}}^{-1/2}$ diverges (see Appendix D.2).

In a future work we will perform a more detailed parameter survey to compare the simplified one-fluid and full two-fluid frameworks in calculating SI (see comparisons using other problems, Laibe & Price 2014; Price & Laibe 2015).

LY15 found that the appropriate way to compare vertical shear (destabilizing) and vertical buoyancy (stabilizing) is $r{\partial }_{z}{\rm{\Omega }}/{{\rm{\Omega }}}_{{\rm{K}}}$ against ${N}_{z}^{2}/{{\rm{\Omega }}}_{{\rm{K}}}^{2}$. From Sections 5.1.2 and 5.1.3 we find $| r{\partial }_{z}{\rm{\Omega }}| \sim {{qh}}_{{\rm{g}}}{{\rm{\Omega }}}_{{\rm{K}}}$ for a thin, gas-dominated disk, while ${N}_{z}^{2}\sim \epsilon {{\rm{\Omega }}}_{{\rm{K}}}^{2}$. Thus we expect dust-induced buoyancy forces to stabilize the disk against the VSI where $\epsilon \gtrsim {h}_{{\rm{g}}}$.

We first vary the midplane dust-to-gas ratio ${\epsilon }_{0}\in [{10}^{-3},1]$, fixing the dust thickness to ${H}_{{\rm{d}}}=0.99{H}_{{\rm{g}}}$. Then is roughly constant with height. We set the vertical domain to ${z}_{\max }=5{H}_{{\rm{g}}}$.

Figure 5 compares the vertical shear rate in the basic state, which is destabilizing, and the vertical buoyancy frequency, which is stabilizing. For the nearly dust-free disk with ${\epsilon }_{0}={10}^{-3}$ the vertical shear dominates buoyancy for all $| z| \gt 0$. However, a heavy dust-load with ${\epsilon }_{0}=1$ causes the buoyancy to dominate over vertical shear in the disk atmosphere $| z| \gtrsim 2.5{H}_{{\rm{g}}}$. We thus expect instability to occur at all heights for ${\epsilon }_{0}={10}^{-3}$, but to be restricted to the midplane for ${\epsilon }_{0}=1$.

Zoom In Zoom Out Reset image size

Figure 6 shows unstable modes for different values of ${\epsilon }_{0}$ with fixed perturbation wavenumber ${k}_{x}{H}_{{\rm{g}}}=30$. The eigenvalue distributions for ${\epsilon }_{0}\leqslant {10}^{-2}$ are similar to the dust-free fiducial case considered by LY15. This is expected since ${\epsilon }_{0}\lt {h}_{{\rm{g}}}$ (Section 6.1). Eigenvalues consists of the roughly horizontal "body modes," and the nearly vertical "surface modes" (which are associated with the imposed vertical boundaries, Barker & Latter 2015).

Zoom In Zoom Out Reset image size

We find that increasing the dust-to-gas ratio reduces VSI growth rates. Notably, surface modes, which are typically fastest growing in the dust-free case, are suppressed in dusty disks for ${\epsilon }_{0}\geqslant 0.1$ (i.e., ${\epsilon }_{0}\gt {h}_{{\rm{g}}}$). The body modes' growth rates remain $O({h}_{{\rm{g}}}{{\rm{\Omega }}}_{{\rm{K}}})$ but their oscillation frequency increases with dust-loading, i.e., it increases with the vertical buoyancy. The total number of modes does not change significantly. This is in contrast with the effect of increasing cooling times in an adiabatic gas disk, for which LY15 find fewer unstable modes.

The lowest frequency "fundamental" body mode is energetically dominant because the entire disk column is perturbed (unlike surface modes, which only disturb the disk boundaries, Umurhan et al. 2016a). In Figure 7 we compare the fundamental mode between the nearly dust-free case ${\epsilon }_{0}={10}^{-3}$ and a dusty disk with ${\epsilon }_{0}=1$. Dust-loading preferentially stabilizes the disk atmosphere against the VSI, restricting meridional motions to $| z| \lesssim 2{H}_{{\rm{g}}}$. This is consistent with Figure 5 comparing the vertical shear in the basic state and the buoyancy.

Zoom In Zoom Out Reset image size

In Figure 8 we plot the growth rates as a function of ${\epsilon }_{0}$ for different perturbation wavenumbers k_x. Dust-loading stabilizes the VSI more effectively for shorter wavelength perturbations. This is because for high wavenumbers the dominant modes are surface modes, which are effectively stabilized by dust-loading because buoyancy forces are largest near the vertical boundaries. The figure suggest that VSI becomes much less efficient for ${\epsilon }_{0}\gtrsim 0.1$ and ${k}_{x}{H}_{{\rm{g}}}\gtrsim 50$.

Zoom In Zoom Out Reset image size

We now vary ${H}_{{\rm{d}}}$ but fix the metallicity $Z\equiv {\epsilon }_{0}{H}_{{\rm{d}}}/{H}_{{\rm{g}}}\,=0.03$ to obtain ${\epsilon }_{0}$. Since we will consider thin dust layers, here we use a smaller domain with ${z}_{\max }=2{H}_{{\rm{g}}}$ so that does not become too small.

We analyze two disks with ${H}_{{\rm{d}}}=0.1{H}_{{\rm{g}}}$ and ${H}_{{\rm{d}}}=0.99{H}_{{\rm{g}}}$. Figure 9 compares the vertical shear rate and buoyancy frequency. For $| z| \gtrsim 0.4{H}_{{\rm{g}}}$ the two disks have the same profile, with vertical shear dominating over buoyancy. We thus expect perturbations away from the disk midplane in both cases. For $| z| \lesssim 0.4{H}_{{\rm{g}}}$, a thin dust layer with ${H}_{{\rm{d}}}=0.1{H}_{{\rm{g}}}$ boosts the vertical shear rate, but the associated buoyancy is larger still, implying that the midplane should be stable.

Zoom In Zoom Out Reset image size

Figure 10 compares the fastest growing VSI body modes with ${k}_{x}{H}_{{\rm{g}}}=30$ for the two cases above. (The thinner domain adopted here eliminates surface modes, LY15.) We find very similar mode frequencies

$$\begin{eqnarray*}\sigma =\left\{\begin{array}{ll}(0.3053{\text{}}i-0.8142){h}_{{\rm{g}}}{{\rm{\Omega }}}_{{\rm{K}}} & {H}_{{\rm{d}}}=0.99{H}_{{\rm{g}}},\\ (0.3178{\text{}}i-1.2237){h}_{{\rm{g}}}{{\rm{\Omega }}}_{{\rm{K}}} & {H}_{{\rm{d}}}=0.1{H}_{{\rm{g}}},\end{array}\right.\end{eqnarray*}$$

since the vertical shear profile is similar throughout most of the disk. However, meridional motions are suppressed near the midplane of the ${H}_{{\rm{d}}}=0.1{H}_{{\rm{g}}}$ disk, as expected from the larger buoyancy frequency relative to vertical shear there. This leads to a structure analogous to PPD dead zones: a quiescent midplane between active surface layers (Gammie 1996).

Zoom In Zoom Out Reset image size

Figure 11 shows the maximum VSI growth rates as a function of ${H}_{{\rm{d}}}$. As before, we find that growth rates are most affected by the vertical structure of the dust layer when the perturbation wavenumber is large. Notice that VSI growth rates converge as ${H}_{{\rm{d}}}\to 0$. Thus a thin dust layer, however large its associated vertical shear, does not affect VSI growth rates. The non-monotonic behavior for ${H}_{{\rm{d}}}\gtrsim 0.5{H}_{{\rm{g}}}$ arises because the vertical buoyancy frequency

$$\begin{eqnarray*}{N}_{z}^{2}({H}_{{\rm{d}}};z,Z) & \simeq & {{ZH}}_{{\rm{g}}}{z}^{2}\exp \left(-\displaystyle \frac{{z}^{2}}{2{H}_{{\rm{g}}}^{2}}\right){{\rm{\Omega }}}_{{\rm{K}}}^{2}\\ & & \times \displaystyle \frac{1}{{H}_{{\rm{d}}}}\left(\displaystyle \frac{1}{{H}_{{\rm{d}}}^{2}}-\displaystyle \frac{1}{{H}_{{\rm{g}}}^{2}}\right)\exp \left(-\displaystyle \frac{{z}^{2}}{2{H}_{{\rm{d}}}^{2}}\right)\end{eqnarray*}$$

is a non-monotonic function of ${H}_{{\rm{d}}}$ at fixed z. At $z={H}_{{\rm{g}}}$ and $z=2{H}_{{\rm{g}}}$ the buoyancy frequency is maximized for ${H}_{{\rm{d}}}\simeq 0.5{H}_{{\rm{g}}}$ and ${H}_{{\rm{d}}}\simeq 0.8{H}_{{\rm{g}}}$, respectively. This is consistent with the minima in growth rates in Figure 11.

Zoom In Zoom Out Reset image size

The discussion in Section 4.3 implies that strictly isothermal disks with a radially uniform dust-to-gas ratio are stable against axisymmetric perturbations, no matter how thin the dust layer is. We now demonstrate this numerically.

To connect with similar studies, here we also use the Richardson number $\mathrm{Ri}\equiv {N}_{z}^{2}/{(r{\partial }_{z}{\rm{\Omega }})}^{2}$ to label calculations (Youdin & Shu 2002). Numerical simulations show that non-axisymmetric instabilities develop when $\mathrm{Ri}\lesssim 0.1$, brought about by very thin dust layers (Chiang 2008; Lee et al. 2010). We show that axisymmetric instabilities never develop in radially uniform disks, however small $\mathrm{Ri}$.

We shall consider ultrathin dust layers with ${H}_{{\rm{d}}}\leqslant 0.01{H}_{{\rm{g}}}$ and thus restrict the vertical domain to ${z}_{\max }=0.02{H}_{{\rm{g}}}$. We fix the metallicity at Z = 0.01 so the midplane dust-to-gas ratio ${\epsilon }_{0}$ is O(1). We consider radial wavenumbers with ${k}_{x}{H}_{{\rm{d}}}=1$.

Figure 12 shows the maximum growth rate as a function of the radial temperature gradient, q. For all cases the vertical shear is dominated by that due to the dust layer (see Section 5.1.2). However, we see that $s\propto | q|$, i.e., growth rates vanish in the strictly isothermal limit. In particular, this holds for $\mathrm{Ri}\lt 0.1$, the critical value for non-axisymmetric instabilities.

Zoom In Zoom Out Reset image size

Axisymmetric instability here is associated with the thermal contribution to vertical shear: as $q\to 0$, ${\partial }_{z}{\rm{\Omega }}$ becomes entirely due to ${\partial }_{z}\epsilon$ and there is no instability (s→0). This result is independent of $\mathrm{Ri}$, so the Richardson number does not characterize the axisymmetric stability of dust layers.

We now consider a radially varying dust-to-gas ratio. Specifically we let ${\epsilon }_{0}\propto {r}^{-1}$ and ${H}_{\epsilon }\propto {H}_{{\rm{g}}}$ (compare with constant values in the previous calculations). Then

$$\begin{eqnarray*}&&\displaystyle \frac{\partial \epsilon }{\partial r}=\left(\displaystyle \frac{{z}^{2}}{{H}_{\epsilon }^{2}}\displaystyle \frac{d\mathrm{ln}{H}_{{\rm{g}}}}{{dr}}-\displaystyle \frac{1}{r}\right)\epsilon .\end{eqnarray*}$$

Here we fix Z = 0.01 and ${H}_{{\rm{d}}}=0.8{H}_{{\rm{g}}}$.

Figure 13 shows that a radially varying dust-to-gas ratio increases the magnitude of the vertical shear rate away from the midplane (see also Equation (46)). Thus we typically find higher VSI growth rates, as shown in Figure 14. Surfaces modes, appearing at high k_x, are more effectively enhanced by the additional vertical shear induced by ${\partial }_{r}\epsilon$. Low-frequency body modes with small k_x are more affected by radial variations in the dust-to-gas ratio than high-frequency body modes with high k_x.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

However, all growth rates remain $O({h}_{{\rm{g}}}{{\rm{\Omega }}}_{{\rm{K}}})$. Importantly, the increase in the growth rate of the fundamental mode is small. We do not expect the radial dependence in to significantly affect the VSI in a dusty disk.

In Section 4.3 we found that for strictly isothermal disks, a vertically uniform dust-to-gas ratio can be unstable if $d\epsilon /{dr}\ne 0$. To demonstrate this numerically we set q = 0 and H_d such that ${H}_{\epsilon }={10}^{3}{H}_{{\rm{g}}}$. We let ${\rho }_{{\rm{d}}}/{\rho }_{{\rm{g}}}\propto {r}^{-d}$ with different power indices d. The metallicity is fixed at Z = 0.03. In these cases vertical shear ${\partial }_{z}{\rm{\Omega }}$ is attributed to the radially varying dust-to-gas ratio (see Equation (46)).

We show unstable modes in Figure 15. As expected from the discussion in Section 4.3, the disk admits purely growing modes with $\omega =0$. This is distinct from the growing oscillations associated with classic VSI discussed above.

Zoom In Zoom Out Reset image size

We find that in smooth disks large k_x is needed for appreciable growth rates, e.g., ${k}_{x}{H}_{g}=1800$ with ${\rho }_{{\rm{d}}}/{\rho }_{{\rm{g}}}\propto {r}^{-1}$. Local shearing-box simulations may be required to study the nonlinear evolution of such short wavelengths (e.g., Bai & Stone 2010b; Yang & Johansen 2016). Alternatively, as shown in Figure 15, a rapidly varying dust-to-gas ratio permits dynamical instability at longer radial wavelengths, which might be resolvable in global disk simulations.

The main application we envision for the dusty/adiabatic gas analogy is to develop physical interpretations of dust–gas drag instabilities, and to find dusty analogs of pure gas instabilities in protoplanetary disks, which is discussed in Section 7.3. One can also exploit the similarity to adapt existing hydrodynamic codes to simulate dusty protoplanetary disks (see Section 7.4).

It is important to keep in mind the assumptions used to develop our thermodynamic model of dusty gas. The terminal velocity approximation, ${{\boldsymbol{v}}}_{{\rm{d}}}-{{\boldsymbol{v}}}_{{\rm{g}}}={t}_{{\rm{s}}}{\rm{\nabla }}P/{\rho }_{{\rm{g}}}$, was employed from the outset. This is applicable to small particles with short stopping times that are strongly coupled to the gas. Generally we require ${t}_{{\rm{s}}}$ to be the shortest timescale in the physical problem. For example, to study dust settling, we require ${t}_{{\rm{s}}}\ll {t}_{\mathrm{settle}}\sim 1/{t}_{{\rm{s}}}{{\rm{\Omega }}}_{{\rm{K}}}^{2}$, or ${t}_{{\rm{s}}}{{\rm{\Omega }}}_{{\rm{K}}}\ll 1$.

However, the validity of the terminal velocity approximation may depend not only on ${t}_{{\rm{s}}}$. The value of the dust–gas ratio and the problem itself may also be important. For example, Table 1 shows that the dust-rich "linA" mode and dust-poor "linB" mode of the SI have the same ${t}_{{\rm{s}}}$, but the former is accurately captured by the simplified equations, while the latter is not. This suggests that for the SI the simplified equations are better suited for ${\rho }_{{\rm{d}}}/{\rho }_{{\rm{g}}}\gt 1$. The simplified equations also contain spurious modes in certain limits (see Appendix D.3).

Note that our thermodynamic interpretation of dust–gas drift does not actually depend on the terminal velocity approximation. Once the equation of state for the gas is fixed and the true energy equation deleted, the dust continuity equation can be converted to a new effective energy equation. As an example, for strictly isothermal dusty gas we obtain

$$\begin{eqnarray}&&\displaystyle \frac{{DP}}{{Dt}}=-P{\rm{\nabla }}\cdot {\boldsymbol{v}}+{c}_{{\rm{s}}}^{2}{\rm{\nabla }}\cdot [{f}_{{\rm{d}}}(1-{f}_{{\rm{d}}})\rho ({{\boldsymbol{v}}}_{{\rm{d}}}-{{\boldsymbol{v}}}_{{\rm{g}}})].\end{eqnarray} \tag{ 62 }$$

This equation does not use the terminal velocity approximation (which gives Equation (15)). Evaluating the right side generally requires solving an evolutionary equation for ${{\boldsymbol{v}}}_{{\rm{d}}}-{{\boldsymbol{v}}}_{{\rm{g}}}$ (Laibe & Price 2014). Nevertheless, Equation (62) shows that dust–gas relative drift can be interpreted as a heat flux within the dust–gas mixture. Thus, even without the terminal velocity approximation we can regard the dusty gas as a single ideal fluid subject to cooling.

We can extend the correspondence between dusty and adiabatic gas to other fixed equations of state. As an example, consider the locally polytropic disk

$$\begin{eqnarray}&&P=K(r,z){\rho }_{{\rm{g}}}^{{\rm{\Gamma }}},\end{eqnarray} \tag{ 63 }$$

where K is a prescribed function and Γ is the constant polytropic index. Then eliminating ${f}_{{\rm{d}}}$ from the dust Equation (6) gives

$$\begin{eqnarray}&&\displaystyle \frac{{DP}}{{Dt}}=-{\rm{\Gamma }}P{\rm{\nabla }}\cdot {\boldsymbol{v}}+P{\boldsymbol{v}}\cdot {\rm{\nabla }}\mathrm{ln}K+\displaystyle \frac{{\rm{\Gamma }}P}{{\rho }_{{\rm{g}}}}{\rm{\nabla }}\cdot ({f}_{{\rm{d}}}{t}_{{\rm{s}}}{\rm{\nabla }}P).\end{eqnarray} \tag{ 64 }$$

Thus the dusty gas behaves like a pure gas with adiabatic index Γ. The entropy is given by

$$\begin{eqnarray}&&{S}_{\mathrm{eff}}=\mathrm{ln}[{K}^{1/{\rm{\Gamma }}}(r,z)(1-{f}_{{\rm{d}}})].\end{eqnarray} \tag{ 65 }$$

The results of Sections 4.3 and 4.4 remain valid for the strictly polytropic disk with $K=$ constant, while one sets ${c}_{{\rm{s}}}^{2}\to K$ in Section 4.5.

7.3.1. Classic GIs

7.3.2. Rossby Wave Instability (RWI)

The RWI is a non-axisymmetric, 2D shear instability that operates in thin disks when it has radial structure (Lovelace et al. 1999; Li et al. 2000). These studies consider adiabatic pure gas and show that instability is possible if there is an extremum in the generalized potential vorticity ${{ \mathcal V }}_{{\rm{g}}}={\kappa }^{2}{{ \mathcal S }}_{{\rm{g}}}^{-2/\gamma }/2{\rm{\Omega }}{{\rm{\Sigma }}}_{{\rm{g}}}$, where ${{ \mathcal S }}_{{\rm{g}}}=P/{{\rm{\Sigma }}}_{{\rm{g}}}^{\gamma }$ is essentially the gas entropy. Here P should be interpreted as the vertically integrated pressure. The nonlinear result of the RWI is vortex formation (Li et al. 2001).

The dusty/adiabatic gas analogy implies that the condition for RWI in a polytropic dusty gas disk is an extremum,

$$\begin{eqnarray}&&{ \mathcal V }=\displaystyle \frac{{\kappa }^{2}}{2{\rm{\Omega }}{\rm{\Sigma }}}\displaystyle \frac{1}{{K}^{2/{\rm{\Gamma }}}}{\left(1+\displaystyle \frac{{{\rm{\Sigma }}}_{{\rm{d}}}}{{{\rm{\Sigma }}}_{{\rm{g}}}}\right)}^{2}.\end{eqnarray} \tag{ 66 }$$

Thus RWI may also be triggered by extrema in the dust-to-gas ratio, e.g., narrow dust rings/gaps. This may lead to direct formation of dusty vortices, as opposed to dust-trapping by a pre-existing gas vortex (Barge & Sommeria 1995; Lyra & Lin 2013).

7.3.3. Convective Overstability (ConO)

7.3.4. Zombie Vortices

The one-fluid model with the terminal velocity approximation, Equations (5)–(7), has been applied to simulate dusty protoplanetary disks (Dipierro et al. 2015; Ragusa et al. 2017). These studies explicitly evolve the dust density.

However, the equivalence of dusty and adiabatic gas identified in this paper means that one need not implement a separate module for simulating the dust component in locally isothermal/polytropic gas. The standard energy equation substitutes for the dust continuity equation. Thus pure gas dynamic codes can be used to simulate protoplanetary dusty disks. When the sound speed c_s (or K) is not constant and/or the dust–gas coupling is imperfect, ${t}_{{\rm{s}}}\ne 0$, one should add corresponding source terms in the energy equation (e.g., Equation (64)). The source term associated with dust–gas drag, ${ \mathcal C }$, is analogous to radiative diffusion (Price & Laibe 2015) or thermal conduction, which is also common in modern hydrodynamics codes.

We have taken advantage of the equivalence of dusty and adiabatic gas to convert the popular Pluto⁵ hydrodynamics code (Mignone et al. 2007) into a dusty gas dynamics code appropriate for simulating protoplanetary disks coupled to small dust. In fact, our code is unaware of the fact that it is modeling dusty gas. We will apply this modified code—dPluto—to study dusty disk–planet interaction (J.-W. Chen & M.-K. Lin, in preparation). Preliminary results show that our approach reproduces features such as dusty rings associated with planet-induced gaps similar to those obtained from explicit two-fluid simulations (e.g., Dong et al. 2017).

We can set $\zeta =0$ or consider ${c}_{{\rm{s}}}^{2}\to \infty$ to obtain the incompressible gas limit. Using ${\tau }_{{\rm{s}}}\equiv {t}_{{\rm{s}}}/(1-{f}_{{\rm{d}}})$, we obtain

$$\begin{eqnarray}&&{{if}}_{{\rm{d}}}{\tau }_{{\rm{s}}}{\sigma }^{4}-{\sigma }^{3}-{\tau }_{{\rm{s}}}[{{if}}_{{\rm{d}}}{\kappa }^{2}+{k}_{x}{F}_{r}(1-{f}_{{\rm{d}}})]{\sigma }^{2}\\ &&\quad +{\left(\displaystyle \frac{{k}_{z}\kappa }{| {\boldsymbol{k}}| }\right)}^{2}\sigma +{k}_{x}{\tau }_{{\rm{s}}}{F}_{r}(1-2{f}_{{\rm{d}}}){\left(\displaystyle \frac{{k}_{z}\kappa }{| {\boldsymbol{k}}| }\right)}^{2}=0,\end{eqnarray} \tag{ 97 }$$

as derived by Jacquet et al. (2011) and Laibe & Price (2014) for the SI. If $| \sigma | /{{\rm{\Omega }}}_{{\rm{K}}}=O({\tau }_{{\rm{s}}}{{\rm{\Omega }}}_{{\rm{K}}})$ and ${\tau }_{{\rm{s}}}{{\rm{\Omega }}}_{{\rm{K}}}\ll 1$, then the quartic term is small and may be neglected. In that case we obtain the cubic dispersion relation of Youdin & Goodman (2005).

Here we seek analytic solutions for the SI by examining limiting cases and with additional approximations. We will fix ${k}_{z},{f}_{{\rm{d}}}$ and maximize growth rates over k_x. It turns out that simple solutions exist in the dust-rich case with ${f}_{{\rm{d}}}$ near unity. We consider small stopping times, ${\tau }_{{\rm{s}}}\to 0$, and assume that the corresponding optimum ${k}_{x}\to \infty$.

We begin with the incompressible dispersion relation, Equation (97). We assume a Keplerian disk (${\rm{\Omega }}=\kappa ={{\rm{\Omega }}}_{{\rm{K}}}$), and low-frequency modes ($| \sigma | \ll {{\rm{\Omega }}}_{{\rm{K}}}$), which allows us to neglect the quartic term. (This is also necessary to avoid spurious modes, see Appendix D.3.) In dimensionless form, the dispersion relation is

$$\begin{eqnarray}&&{\nu }^{3}+\mathrm{St}[{{if}}_{{\rm{d}}}+2{K}_{x}{f}_{{\rm{g}}}(1-{f}_{{\rm{d}}})]{\nu }^{2}-{\left(\displaystyle \frac{{K}_{z}}{{K}_{x}}\right)}^{2}\nu \\ &&\quad -2\mathrm{St}{K}_{x}{f}_{{\rm{g}}}(1-2{f}_{{\rm{d}}}){\left(\displaystyle \frac{{K}_{z}}{{K}_{x}}\right)}^{2}=0,\end{eqnarray} \tag{ 98 }$$

where $\nu =\sigma /{{\rm{\Omega }}}_{{\rm{K}}}$, $\mathrm{St}={\tau }_{{\rm{s}}}{{\rm{\Omega }}}_{{\rm{K}}}$, ${f}_{{\rm{g}}}=1-{f}_{{\rm{d}}}$ (the gas fraction), and recall that ${K}_{x,z}=\eta {{rk}}_{x,z}$. We have used ${F}_{r}=2\eta {f}_{{\rm{g}}}r{{\rm{\Omega }}}^{2}$ and considered large K_x.

Analytic expressions for cubic roots are unwieldy. To keep the problem tractable, let us assume at this stage that the quadratic term can be neglected compared to the last term. That is,

$$\begin{eqnarray}&&\left|\displaystyle \frac{{{if}}_{{\rm{d}}}}{2{K}_{x}{f}_{{\rm{g}}}}+(1-{f}_{{\rm{d}}})\right|{| {K}_{x}\nu | }^{2}\ll {K}_{z}^{2}(2{f}_{{\rm{d}}}-1).\end{eqnarray} \tag{ 99 }$$

(The imaginary term can be neglected for fixed ${f}_{{\rm{d}}}$ but allowing ${K}_{x}\to \infty$.) We show in Appendix D.2.2 that this simplification is valid for dust-rich disks.

Equation (98) now becomes the depressed cubic

$$\begin{eqnarray}&&{\nu }^{3}-{ \mathcal P }\nu ={ \mathcal Q },\end{eqnarray} \tag{ 100 }$$

where ${ \mathcal P }$ and ${ \mathcal Q }$ can be read off Equation (98). This can be solved with Vieta's substitution

$$\begin{eqnarray}&&\nu =\mu +\displaystyle \frac{{ \mathcal P }}{3\mu },\end{eqnarray} \tag{ 101 }$$

then Equation (100) becomes a quadratic for ${\mu }^{3}$,

$$\begin{eqnarray}&&{\mu }^{6}-{ \mathcal Q }{\mu }^{3}+\displaystyle \frac{{{ \mathcal P }}^{3}}{27}=0.\end{eqnarray} \tag{ 102 }$$

The solution is

$$\begin{eqnarray}{\mu }^{3} & = & \displaystyle \frac{{ \mathcal Q }}{2}\left(1\pm \sqrt{1-\displaystyle \frac{4{{ \mathcal P }}^{3}}{27{{ \mathcal Q }}^{2}}}\right)\\ & = & \mathrm{St}\,{f}_{{\rm{g}}}(1-2{f}_{{\rm{d}}})\displaystyle \frac{{K}_{z}^{2}}{{K}_{x}}\\ & & \times \left\{1\pm \sqrt{1-\displaystyle \frac{{K}_{z}^{2}}{27{K}_{x}^{2}{[\mathrm{St}{K}_{x}{f}_{{\rm{g}}}(1-2{f}_{{\rm{d}}})]}^{2}}}\right\}.\end{eqnarray} \tag{ 103 }$$

At this point we assume that the term $\propto {K}_{z}^{2}$ inside the square root may be neglected. That is,

$$\begin{eqnarray}&&{K}_{z}^{2}\ll 27{({K}_{x}^{2}\mathrm{St})}^{2}{[{f}_{{\rm{g}}}(1-2{f}_{{\rm{d}}})]}^{2}.\end{eqnarray} \tag{ 104 }$$

We show in Appendix D.2.2 that this is readily satisfied. Then $\mu \simeq {{ \mathcal Q }}^{1/3}$.

Now, remembering that we are considering the dust-rich limit with $1-2{f}_{{\rm{d}}}\lt 0$, the explicit solution for μ is

$$\begin{eqnarray}&&\mu ={\left[2\mathrm{St}{f}_{{\rm{g}}}(2{f}_{{\rm{d}}}-1)\displaystyle \frac{{K}_{z}^{2}}{{K}_{x}}\right]}^{1/3}{e}^{i\pi /3}.\end{eqnarray} \tag{ 105 }$$

We have chosen the complex root for instability. Inserting this into Equation (101) gives the eigenfrequency ν:

$$\begin{eqnarray}&&\mathrm{Re}(\nu )=\displaystyle \frac{1}{2}({ \mathcal A }{\mathrm{St}}^{1/3}{K}_{x}^{-1/3}+{ \mathcal B }{\mathrm{St}}^{-1/3}{K}_{x}^{-5/3}),\end{eqnarray} \tag{ 106 }$$

$$\begin{eqnarray}&&\mathrm{Im}(\nu )=\displaystyle \frac{\sqrt{3}}{2}({ \mathcal A }{\mathrm{St}}^{1/3}{K}_{x}^{-1/3}-{ \mathcal B }{\mathrm{St}}^{-1/3}{K}_{x}^{-5/3}),\end{eqnarray} \tag{ 107 }$$

where

$$\begin{eqnarray}&&{ \mathcal A }={[2{f}_{{\rm{g}}}(2{f}_{{\rm{d}}}-1){K}_{z}^{2}]}^{1/3},\end{eqnarray} \tag{ 108 }$$

$$\begin{eqnarray}&&{ \mathcal B }=\displaystyle \frac{{K}_{z}^{2}}{3{ \mathcal A }}.\end{eqnarray} \tag{ 109 }$$

Maximizing growth rates over K_x by setting $\partial \mathrm{Im}(\nu )/\partial {K}_{x}=0$, we find the optimum wavenumber is given by

$$\begin{eqnarray}&&{K}_{x}={\left(\displaystyle \frac{5{ \mathcal B }}{{ \mathcal A }}\right)}^{3/4}{\mathrm{St}}^{-1/2}={\left(\displaystyle \frac{5}{3}\right)}^{3/4}\sqrt{\displaystyle \frac{{K}_{z}}{2{f}_{{\rm{g}}}(2{f}_{{\rm{d}}}-1)\mathrm{St}}}.\end{eqnarray} \tag{ 110 }$$

The maximum growth rate is then

$$\begin{eqnarray}\max [\mathrm{Im}(\nu )] & = & \displaystyle \frac{2\sqrt{3}{{ \mathcal A }}^{5/4}}{{5}^{5/4}{{ \mathcal B }}^{1/4}}\sqrt{\mathrm{St}}\\ & = & \displaystyle \frac{2\sqrt{2}\times {3}^{3/4}}{{5}^{5/4}}\sqrt{{f}_{{\rm{g}}}(2{f}_{{\rm{d}}}-1){K}_{z}\mathrm{St}}.\end{eqnarray} \tag{ 111 }$$

and the real frequency $\mathrm{Re}(\nu )=(\sqrt{3}/2)\mathrm{Im}(\nu )$. We see that as ${\tau }_{{\rm{s}}}\to 0$, the most unstable radial wavenumber diverges, ${K}_{x}\,\propto {\tau }_{{\rm{s}}}^{-1/2}\to \infty$ with a vanishing growth rate, $\max (s)\propto \sqrt{{\tau }_{{\rm{s}}}}\to 0$. This corresponds to instability at arbitrarily small radial length scales.

D.2.1. Finite Phase Lag as ${\tau }_{{\rm{s}}}\to 0$

An interesting property of the special solutions described above is that there is always a phase lag between the Lagrangian pressure and density perturbations. Since $\mathrm{Re}(\nu )\gt 0$, the phase lag may be defined as $\varphi =\arg ({\rm{\Delta }}P{\rm{\Delta }}{\rho }^{* })$. (Note that we may set ${\rm{\Delta }}\rho =1$ without loss of generality.) It may be shown that ${\rm{\Delta }}P\,\simeq i\delta {v}_{x}{\partial }_{r}P/\sigma$ in the limit of large K_x. To obtain $\delta {v}_{x}$ we use Equation (90) in the low-frequency limit. The expression for $\delta {v}_{x}$ involves the Eulerian pressure perturbation, $\delta P$, for which we use Equation (94) and assume incompressibility ($\zeta \to 0$).

With these additional simplifications the phase lag φ is given via

$$\begin{eqnarray}&&\tan \varphi =\displaystyle \frac{\mathrm{Im}({\rm{\Delta }}P{\rm{\Delta }}{\rho }^{* })}{\mathrm{Re}({\rm{\Delta }}P{\rm{\Delta }}{\rho }^{* })}\simeq \displaystyle \frac{{K}_{x}\mathrm{Im}(\nu )}{2{f}_{{\rm{g}}}^{2}\mathrm{St}{K}_{x}^{2}+{K}_{x}\mathrm{Re}(\nu )}.\end{eqnarray} \tag{ 112 }$$

From here it is clear for the above solutions, where ${K}_{x}\propto {\mathrm{St}}^{-1/2}$ and $\nu \propto {\mathrm{St}}^{1/2}$, that φ is a constant. Explicitly inserting the solutions gives

$$\begin{eqnarray}&&\tan \varphi =\displaystyle \frac{2\times {3}^{3/2}(2{f}_{{\rm{d}}}-1)}{16-7{f}_{{\rm{d}}}}.\end{eqnarray} \tag{ 113 }$$

For the case shown in Figure 4 with ${f}_{{\rm{d}}}=0.75$, we have $\varphi \simeq 26^\circ$, comparable to the $\sim 30^\circ$ obtained from numerical solutions.

D.2.2. Consistency Check

Here we check if the above explicit solutions are consistent with the assumptions used to obtain them. Inserting the solutions for the most unstable mode, we find Equation (99) becomes

$$\begin{eqnarray}&&{f}_{{\rm{d}}}\gg \displaystyle \frac{12}{17}\simeq 0.7.\end{eqnarray} \tag{ 114 }$$

Since ${f}_{{\rm{d}}}\leqslant 1$, this requirement can only be marginally satisfied. However, we find that this mostly gives errors in the real frequency (see Figure 4), while growth rates are still captured correctly. On the other hand, the assumption of Equation (104) becomes the trivial inequality

$$\begin{eqnarray}&&1\ll \displaystyle \frac{27}{4}{\left(\displaystyle \frac{5}{3}\right)}^{3}\simeq 30,\end{eqnarray} \tag{ 115 }$$

so Equation (104) is satisfied, which justifies the approximate solution for μ in Equation (105).

Ultimately, the validity of these assumptions is justified a posteriori by comparison with the solution to the full equations in Figure 4.

A caveat of the dispersion relations Equations (96) and (97) is that they admit spurious unstable modes with $| \sigma | \gtrsim {\rm{\Omega }}$. This violates the one-fluid approximation to model dusty gas. We demonstrate this below by considering the incompressible dispersion relation. (We checked numerically that compressibility has negligible effects on the modes examined.)

Consider the limit k_z = 0. Then Equation (97) becomes

$$\begin{eqnarray}&&i{\nu }^{2}-\displaystyle \frac{\nu }{{f}_{{\rm{d}}}\mathrm{St}}-\left[i+2{K}_{x}\displaystyle \frac{{f}_{{\rm{g}}}}{{f}_{{\rm{d}}}}(1-{f}_{{\rm{d}}})\right]=0\end{eqnarray} \tag{ 116 }$$

in dimensionless form. Neglecting the quadratic term leads to stability, $\mathrm{Im}(\nu )\lt 0$. This is consistent with a full two-fluid analysis (Youdin & Goodman 2005).

However, solving Equation (116) explicitly assuming $| {f}_{{\rm{d}}}{\tau }_{{\rm{s}}}\mathrm{Im}(\nu )| \ll 1$ would yield

$$\begin{eqnarray}&&\mathrm{Re}(\nu )\simeq 2{K}_{x}{(1-{f}_{{\rm{d}}})}^{2}\mathrm{St}\sim O(1)\end{eqnarray} \tag{ 117 }$$

$$\begin{eqnarray}&&\mathrm{Im}(\nu )\simeq {f}_{{\rm{d}}}\mathrm{St}[4{K}_{x}^{2}{(1-{f}_{{\rm{d}}})}^{4}{\mathrm{St}}^{2}-1].\end{eqnarray} \tag{ 118 }$$

Accordingly, growth is possible if

$$\begin{eqnarray}&&\mathrm{St}\gt \displaystyle \frac{1}{2{K}_{x}{(1-{f}_{{\rm{d}}})}^{2}}.\end{eqnarray} \tag{ 119 }$$

Alternatively, for fixed ${\tau }_{{\rm{s}}}$ growth is enabled by a sufficiently large K_x. Unlike the SI, which is strongly suppressed when ${f}_{{\rm{d}}}=1/2$, these modes can still grow at equal dust-to-gas ratio.

These growing epicycles with $| \omega | \gtrsim \kappa$ are absent from full two-fluid models (Youdin & Goodman 2005). The discrepancy lies in the fact that the one-fluid equations, to first order in ${\tau }_{{\rm{s}}}$, are only valid for low-frequency waves with $| \sigma | \lesssim {\rm{\Omega }}$. Thus only low-frequency modes should be retained from analyses based on Equations (5)–(8).

Figure 16 shows unstable modes found from Equation (96) as a function of K_z at fixed K_x = 1500, $\epsilon =2$, and ${\tau }_{{\rm{s}}}{{\rm{\Omega }}}_{{\rm{K}}}=0.01$. The growth rates of the (spurious) overstable dusty epicycles are weakly dependent on K_z and are well approximated by that in the limit K_z = 0, Equation (118). For the case considered in Figure 16 we find $s\simeq 0.067{{\rm{\Omega }}}_{{\rm{K}}}$, as observed. By contrast, the SI requires ${K}_{z}\gt 0$, and it dominates when ${K}_{z}\gtrsim 100$.

Zoom In Zoom Out Reset image size

Care must be taken if the first-order one-fluid equations are used to simulate dusty gas. For example, 2D, razor-thin disks would allow the spurious epicycles to dominate, since in that case the SI cannot operate. However, Figure 16 shows that the SI should dominate in realistic 3D disks where a range of K_z is allowed.

Thus simulations based on the first-order one-fluid equations should be set up to suppress these spurious epicycles. This might be achieved, for example, through physical or numerical viscosity to eliminate high-k_x modes, since for fixed disk/dust parameters these spurious epicycles operate only at sufficiently small radial wavelengths. Alternatively, one needs to ensure that the physical instabilities of interest have larger growth rates than the spurious epicycles.