Molecular clouds as gravitational instabilities in rotating disks: a modified stability criterion

A major motivation for including the vertical dimension is to obtain a realistic description of molecular gas disks in which turbulent gas motions are three-dimensional and nearly isotropic and in which the embedded clouds are triaxial. As will be shown in $\S$ 3.4.1, the 3D dispersion relation derived in this work approaches the 2D Lin-Shu dispersion relation in the limit $\sigma_{z}\rightarrow 0$. For the more general scenario of interest here, both vertical and radial components of motion are allowed, each given the following 1D velocity dispersion (i.e. Chandrasekhar, 1951)

$$\sigma_{\rm i}^{2}=v_{s}^{2}+\sigma_{\rm turb,i}^{2}$$

(3)

which combines the sound speed in the gas $v_{s}$ with turbulent motions $\sigma_{\rm turb,i}$ in direction $i$. Although we allow that $\sigma_{z}$$\neq$$\sigma_{r}$, velocity dispersions in the gas are generally assumed to be isotropic.

The turbulent motions in the gas are envisioned as arising from two main sources that combine to yield an effective (non-thermal) pressure that places the disk in dynamical equilibrium. Star formation feedback (plus turbulent dissipation) is assumed to set a base pressure $p_{\rm FB}$. This combines with the effective pressure $p_{\rm eff}$ set up by the averaged kinematic response of many individual fluid elements to the force set up by the remainder between the gravitational force and the gradient in the baseline $p_{\rm FB}$. For a given $p_{\rm FB}$, $p_{\rm eff}$ is thus the pressure that maintains the gas disk in overall equilibrium. Thus, in what follows, dynamical equilibrium is applied even when feedback-driven turbulent pressure is either zero or very small due to the absence star formation.

For this equilibrium scenario, unless otherwise noted, the unperturbed vertical density distribution is envisioned as falling between the self-gravitating profile

$$\rho_{0}(z)=\rho_{c}sech^{2}(z/h),$$

(4)

where $h=\sigma_{z}/(2\pi G\rho_{c})^{1/2}$, and

$$\rho_{0}(z)=\rho_{c}e^{-z^{2}/h^{2}}$$

(5)

in the presence of a dominant external potential generated by the background distribution with density $\rho_{b}$, where $h=\sigma_{z}/(4\pi G\rho_{b})^{1/2}$. The equilibrium vertical velocity dispersion associated with these profiles is constant (independent of $z$).

Allowing the gas disk to have some thickness, we now wish to incorporate realistic perturbations that are compatible with the sorts of influences that a gas disk experiences. The goal with these perturbations is not to represent a specific process but rather to invoke a generic influence that is non-negligible away from the galactic mid-plane in a manner that satisfies the linearized equations of motion in 3D. Thus, we adopt wave perturbations of the form

$$\Phi_{1}(R,\phi,z,t)=Re[\mathcal{F}(R,z)e^{i(m\phi-\omega t)}e^{ikR}e^{ik_{z}z}]$$

(6)

or

$$\Phi_{1}(R,\phi,z,t)=iIm[\mathcal{F}(R,z)e^{i(m\phi-\omega t)}e^{ikR}e^{ik_{z}z}]$$

(7)

where the wavenumbers $k$=$2\pi/\lambda_{r}$ and $k_{z}$=$2\pi/\lambda_{z}$ describe the wavelengths of the perturbation in the radial and vertical directions, respectively.

As is typical, these perturbations are assumed to satisfy the WKB approximation in the radial direction, i.e. $\partial\mathcal{F}(R,z)/\partial R=\mathcal{F}(R,z)/R_{p}<<ik\mathcal{F}(R,z)$ where $R_{p}$ is the characteristic scale over which the amplitude of the perturbation varies in the radial direction. In practice this is adopted as the criterion $kR>>1$.

In the vertical direction, three different scenarios are considered. In the first two scenarios, the perturbations are assumed to be infinite, pervading every location in the disk, but either the wave nature is neglected in the vertical direction and so $k_{z}=0$ (the first scenario, as considered by GLB) or $k_{z}$ is assumed to be a constant and independent of height $z$ above the mid-plane (the second scenario). In both of these cases, the amplitude of the perturbation must vary faster in the vertical direction than the vertical density variation in the unperturbed disk, in order that the perturbation remains small with respect to $\rho_{0}$ everywhere ($\rho_{1}<<\rho_{0}$). Defining $Z_{d}=(d\textrm{ln}\Phi_{0}/dz)^{-1}$ and $Z_{p}$=$(d\textrm{ln}\mathcal{F}(R,z)/dz)^{-1}$, then for these infinite perturbations, $Z_{p}\lesssim Z_{d}$. As in $\S\S$ 3.2 and 3.3 this can be cast in terms of the vertical gradients in the unperturbed and perturbed densities, i.e. $z_{p}=z_{d}$ with $z_{d}=(d\textrm{ln}\rho_{0}/dz)^{-1}$ and $z_{p}$=$(d\textrm{ln}\rho_{1}/dz)^{-1}$.

The wave type of these infinite perturbations are not examined in the WKB approximation, given that they entail rapid variation in the perturbed amplitude. Thus, any infinitely extended 3D perturbation considered in this work satisfies the most general form of the perturbed Poisson’s equation

$$\Phi_{1}=\frac{4\pi G\rho_{1}}{-k_{\textrm{plane}}^{2}-k_{z}^{2}+T^{2}}$$

(8)

where $k_{\textrm{plane}}^{2}\equiv k^{2}+m^{2}/R^{2}$ and

	$\displaystyle T^{2}$	$\displaystyle\equiv$	$\displaystyle(\nabla_{z}^{2}\Phi_{1})/\Phi_{1}+k_{z}^{2}$		(9)
		$\displaystyle=$	$\displaystyle\nabla_{z}\left(\frac{1}{Z_{p}}\right)+\frac{1}{Z_{p}^{2}}+\frac{2ik_{z}}{Z_{p}}.$

In the third scenario, wave perturbations are restricted to a finite extent above and below the mid-plane. Such perturbations can be studied without the requirement that $Z_{p}=Z_{d}$ (or $z_{p}=z_{d}$), since they are not at risk of becoming non-negligible as the unperturbed disk density drops towards $z\rightarrow\pm\infty$. These perturbations could vary arbitrarily in amplitude in the $z$ direction (as long as this variation is negligible with respect to the perturbation’s phase variation, in the WKB approximation) and would thus not need to satisfy $k_{z}z_{d}>>1$, or necessarily share the overall variation of the unperturbed disk. The amplitude of the perturbation might instead vary much more slowly than $\rho_{0}$, perhaps tied to the density distribution of an embedding disk or a process active therein, for instance. In this scenario, the WKB approximation is invoked as $k_{z}z_{p}>>1$ (or $T\approx 0$).¹¹1These finite perturbations could be selected to also satisfy $k_{z}z_{d}>>1$ (although it would not be necessary be design). This might be equivalent to $k_{z}z>>1$, in the case of a flattened logarithmic potential $\Phi_{0}\propto\ln(R^{2}+z^{2}/q^{2})$, for example, or $k_{z}h>>1$ in the case of an exponential vertical distribution with scale length $h$.

In practice, finite (WKB) wave perturbations are described in what follows by introducing a truncation at some height $h_{1}$, above which the density becomes zero. (Note that, beyond $h$, $\rho_{1}$ is still assumed to be considerably less than $\rho_{0}$.) Such finite WKB perturbations are maximally flexible as they require only $k_{z}z_{p}>>1$ and not $k_{z}z_{d}>>1$ or $k_{z}h>>1$.

Introducing a truncation in the perturbation at some height $h_{1}$ comes with one important additional requirement. To make the perturbation physical, it must satisfy both Poisson’s equation and Laplace’s equation beyond $h_{1}$. This introduces a strict boundary condition at the interface $|z|=h_{1}$, which places restrictions on the relationship between the vertical and radial perturbations. This condition is determined below by matching the solution to Poisson’s equation with the solution to Laplace’s equation at $z=h_{1}$ and requiring a similar matching of the gravitational force $\partial\Phi_{1}/\partial z$ at the interface (so that the gravitational force remains smooth). Here $h_{1}$ is taken to be the disk scale height or greater.

To satisfy Laplace’s equation

$$-k^{2}\Phi_{1}-\frac{m^{2}}{R^{2}}\Phi_{1}+\frac{\partial^{2}\Phi_{1}}{\partial z^{2}}=0$$

(10)

above and below the perturbation (at and beyond the vertical extent), the solution must have $k_{z}$=$ik_{\textrm{plane}}$. Outside the perturbed part of the disk the potential thus becomes a decaying function $\Phi_{1}$$\propto e^{-|k_{\textrm{plane}}z|}$.

Meanwhile, over the extent of the perturbation, Poisson’s equation

$$-k_{\textrm{plane}}^{2}\Phi_{1}+\frac{\partial^{2}\Phi_{1}}{\partial z^{2}}=4\pi G\rho_{1}$$

(11)

implies that

$$\Phi_{1}=\frac{-4\pi G\rho_{1}}{k_{\textrm{plane}}^{2}+k_{z}^{2}}$$

(12)

in the WKB approximation, where again $k_{\textrm{plane}}^{2}=k^{2}+m^{2}/r^{2}$ and now $\partial^{2}\Phi_{1}/\partial z^{2}\approx-k_{z}^{2}\Phi_{1}$. For Lin-Shu perturbations that are confined to an infinitssimally thin sheet $\Phi_{1}=-2\pi G\Sigma_{1}/|k_{\textrm{plane}}|$, in contrast (Toomre, 1964).

Below both even and odd density and potential wave perturbations are considered, although even perturbations are the focus of the remainder of the paper. (Odd wave perturbations can be shown to yield a consistent view of the main features of 3D disk instability.) As discussed later in $\S$ 3.4.1, the 2D Lin-Shu dispersion relation and Toomre criterion are retrieved adopting an even wave perturbation, which can be envisioned as an over-density in the galaxy mid-plane. This is the nominal perturbation for describing, i.e., the propagation of density waves in the disk. The amplitudes of all perturbations considered in this work (whether even or odd) are assumed to be even functions of distance from the mid-plane. Perturbations with amplitudes that are odd functions of $z$ have been shown to be stable by GLB.

In the case that the potential perturbation in the disk is an even function $\Phi_{1}\propto\cos{k_{z}z}$ and symmetric about the mid-plane with extent $h_{1}$, we obtain the following matching conditions

	$\displaystyle Ae^{-k_{\textrm{plane}}h_{1}}$	$\displaystyle=$	$\displaystyle Bcos(k_{z}h_{1})$
	$\displaystyle-k_{\textrm{plane}}Ae^{-k_{\textrm{plane}}h_{1}}$	$\displaystyle=$	$\displaystyle-k_{z}Bsin(k_{z}h_{1})$

at the interface $h_{1}$, where $A$ is the amplitude of the potential perturbation beyond $|z|=h_{1}$ and $B$ is the amplitude within $h_{1}$. These conditions yields the relation (see also Griv & Gedalin, 2012):

$$\arctan{\frac{k_{\textrm{plane}}}{k_{z}}}=k_{z}h_{1}$$

(13)

In the standard long wavelength scenario with $k_{\textrm{plane}}<<k_{z}$, eq.(13) reduces to $k_{\textrm{plane}}/k_{z}\approx k_{z}h_{1}$. Taking $h_{1}$=$h$, this yields the approximation

$$\frac{\rho_{0}k_{\textrm{plane}}^{2}}{k_{\textrm{plane}}^{2}+k_{z}^{2}}\approx\frac{\Sigma_{0}k_{\textrm{plane}}}{2(1+k_{\textrm{plane}}h)}$$

(14)

in terms of the unperturbed disk volume density $\rho_{0}$ and surface density $\Sigma_{0}\approx 2h\rho_{0}$. This approximation is analogous to corrections incorporated into 2D dispersion relations (of single or multi-component disks) to account for weakened self-gravity when finite thickness is assumed (Toomre, 1964; Vandervoort, 1970; Jog & Solomon, 1984; Romeo, 1992; ylin05; Elmegreen, 2011).

This paper is also interested in the short wavelength regime where the above approximation is not valid. Here ‘short’ is used in relation the vertical wavelength, not the disk scale height. (Instability is indeed still restricted below the Jeans length.) Most relevant in this scenario is the extent of the perturbation $h_{1}$. In the limit $k_{\textrm{plane}}>>k_{z}$, the boundary condition in eq. (13) requires $k_{z}\approx(\pi/2)h_{1}^{-1}$ to lowest order, or that $\lambda_{R}$$<<$$\lambda_{z}\approx h_{1}$. This scenario is consistent with $\lambda_{z}>h$ when we envision the perturbation’s vertical edges at $h_{1}>h$, and it can thus be used to probe the instability regime in which $\lambda_{R}$ is brought down near the size of disk scale height.

In the case that the potential perturbation in the disk is an odd function, $\Phi_{1}\propto\sin{k_{z}z}$, the matching conditions at $h_{1}$ above the plane become

	$\displaystyle Ae^{-k_{\textrm{plane}}h_{1}}$	$\displaystyle=$	$\displaystyle Bsin(k_{z}h_{1})$
	$\displaystyle-k_{\textrm{plane}}Ae^{-k_{\textrm{plane}}h_{1}}$	$\displaystyle=$	$\displaystyle k_{z}Bcos(k_{z}h_{1})$

requiring that

$$\arctan{\frac{-k_{z}}{k_{\textrm{plane}}}}=k_{z}h_{1}.$$

(15)

Similarly, below the plane at -$h_{1}$, the boundary condition requires

$$\arctan{\frac{k_{z}}{k_{\textrm{plane}}}}=k_{z}h_{1}.$$

(16)

Thus in the long wavelength scenario $|k_{\textrm{plane}}|<<|k_{z}|$, allowed perturbations have a vertical wavenumber $k_{z}\sim(\pi/2)h_{1}^{-1}$ and radial wavenumbers are restricted to $|k_{\textrm{plane}}|$$<<$$1/h_{1}$ or $|k_{\textrm{plane}}|$$<$$1/h$ when $h_{1}>h$. Perturbations in the short-wavelength limit $|k_{\textrm{plane}}|>>|k_{z}|$ have $|k_{\textrm{plane}}|\approx 1/h_{1}$ and $k_{z}<<1/h_{1}$, which again can correspond to the case $|kh|\lesssim 1$ and $|k_{z}h|<<1$ where $h_{1}>h$.

This section introduces the 3D perturbations from the previous section (corresponding to infinite non-periodic perturbations, infinite waves and finite WKB waves) into the linearized 3D equations of motion to solve for the perturbed motions in the radial, azimuthal and vertical directions. The 3D dispersion relation is then obtained using the continuity equation, which couples these motions to the time evolution of the perturbed density.

Using either the 3D version of the dispersion relation derived below ($\S$ 2.3) or the 2D version ($\S$ 3), the conditions for stability can be easily determined according to the expectation that a stable, non-growing mode (with Real $\omega$) must have $\omega^{2}$$>$0. (Thus the line of stability is usually taken to be $\omega^{2}$=0; e.g. Binney & Tremaine 1987) In the interest of diagnosing the basic stability of disks to 3D perturbations, sections 2.3 and on examine in detail the axisymmetric scenario with $m$=0, in which case $k_{\textrm{plane}}=k$. The calculations presented in what immediately follows, though, adopt an arbitrary $m$.

To describe motions in our rotating gas disk, we adopt the Euler equations of motion in cylindrical coordinates, with $z$ oriented parallel to the axis of rotation. We then introduce a small perturbation. Writing all quantities as the sum of perturbed and small unperturbed components (i.e. $\rho=\rho_{0}+\epsilon\rho_{1}$ and $v_{R}=v_{R,0}+\epsilon v_{R,1}$, etc., where $\epsilon$ is small) and keeping only terms to first order in small quantities, the linearized versions of the equations of motion are obtained (see Binney & Tremaine, 1987). These are satisfied in this work by the perturbations introduced in $\S$ 2.1.2 with the form $\Phi_{1}(R,\phi,z,t)=\Phi_{a}(R,z)e^{i(m\phi-\omega t)}$ where $\Phi_{a}(R,z)=\mathcal{F}(R,z)e^{ikR+ik_{z}z}$ and the radial gradient of $\mathcal{F}(R,z)$ is neglected in the WKB approximation. Through Poisson’s equation the density perturbation has a similar dependence, i.e. $\rho_{1}=\rho_{a}(R,z)e^{i(m\phi-\omega t)}$ where $\rho_{a}(R,z)=\mathcal{R}(R,z)e^{ikR+ik_{z}z}$. Solutions to the linearized equations of motion (eq. (2)) thus also have the form $v_{R,1}=v_{R,a}(R,z)e^{i(m\phi-\omega t)}$, $v_{\phi,1}=v_{\phi,a}(R,z)e^{i(m\phi-\omega t)}$ and $v_{z,1}=v_{z,a}(R,z)e^{i(m\phi-\omega t)}$ where $v_{R,a}(R,z)$, $v_{R,a}(R,z)$ and $v_{r,a}(R,z)$ are all $\propto e^{ikR+ik_{z}z}$.

Substituting the density and potential perturbations into the perturbed radial and azimuthal equations of motion, it can be shown (adopting the convention in Binney & Tremaine, 1987) that

	$\displaystyle v_{R,a}$	$\displaystyle=$	$\displaystyle-\frac{(\Phi_{a}+\sigma^{2}\frac{\rho_{a}}{\rho_{0}})}{\Delta}\left(k(\omega-m\Omega)+i\frac{2m\Omega}{R}\right)$		(17)
		$\displaystyle-$	$\displaystyle v_{z,a}\frac{2\Omega}{\Delta}\frac{dV_{c}}{dz}$

	$\displaystyle v_{\phi,a}$	$\displaystyle=$	$\displaystyle-\frac{(\Phi_{a}+\sigma^{2}\frac{\rho_{a}}{\rho_{0}})}{\Delta}\left(2Bik+\frac{m(\omega-m\Omega)}{R}\right)$		(18)
		$\displaystyle+$	$\displaystyle iv_{z,a}\frac{(\omega-m\Omega)}{\Delta}\frac{dV_{c}}{dz}$

where $V_{c}=\Omega R$ and

	$\displaystyle B$	$\displaystyle=$	$\displaystyle-\Omega-\frac{1}{2}R\frac{d\Omega}{dR}$		(19)
	$\displaystyle\Delta$	$\displaystyle=$	$\displaystyle\kappa^{2}-(m\Omega-\omega)^{2}$
	$\displaystyle\kappa$	$\displaystyle=$	$\displaystyle-4B\Omega.$

In the vertical direction, the linearized equation of motion

$$i(-\omega+m\Omega)v_{z,1}=-\nabla\Phi_{1}-\sigma^{2}\left(\frac{\nabla\rho_{1}}{\rho_{0}}\right)$$

(20)

implies

$$v_{z,a}=-\frac{k_{z}(\Phi_{a}+\sigma^{2}\frac{\rho_{a}}{\rho_{0}})}{(m\Omega-\omega)}+i\frac{(\nabla\mathcal{F}+\sigma^{2}\frac{\nabla\mathcal{R}}{\rho_{0}})}{(m\Omega-\omega)}e^{ikR+ik_{z}z}.$$

(21)

where, at this stage, no assumption has been made about the relative sizes of the vertical perturbation’s phase and amplitude variations. Different choices for $k_{z}$ in relation to the perturbation amplitude and the unperturbed disk will be examined later in this work.

These expressions for $v_{r,a}$, $v_{\phi,a}$ and $v_{z,a}$ are based on the assumption of equilibrium in the unperturbed disk, such that $v_{r,0}$=0, $v_{z,0}$=0 and $v_{\phi,0}$$\approx(Rd\Phi_{0}/dR)^{1/2}=\Omega R$, neglecting the pressure term $(\partial p_{0,\phi}/\partial R)/\rho_{0}$ since $\sigma_{\phi,0}$$<<$$\Omega R$.

Adopting the WKB approximation in the radial direction leads to further simplification. This work focuses on scenarios in which the radial variation in the amplitude of the potential perturbation is comparable to (and no less than) the radial gradient in the unperturbed disk (as discussed in 2.1.2). As is typical, then, the factors proportional to $1/R$ in eqs. (17) and (18) are neglected relative to those that are proportional to $k$ (e.g. Binney & Tremaine, 1987). The WKB condition $kR_{p}>>1$ ($\S$ 2.1.2) is satisfied by $kR>>$1 assuming $k$ increases towards small $R$.

The terms proportional to $v_{z,1}$ in the expressions for $v_{r,1}$ and $v_{\phi,1}$ are similarly neglected in the set-up of interest here, since the rotational lag

$$\frac{dV_{c}}{dz}\approx\frac{1}{2\Omega}\frac{d}{dz}\frac{d}{dR}\Phi_{0}$$

(22)

(again assuming that the radial pressure gradient is negligible) contains a factor considerably smaller than $k\Phi_{a}$.²²2Appendix A identifies the precise set of perturbations for which the lag term is negligible. .

With an identical vertical perturbation specifically in the WKB approximation, Griv & Gedalin (2012) arrive at a different expression for the perturbed vertical velocity $v_{z,1}$, as the disk in their scenario of interest is out of hydrostatic equilibrium. This introduces factors proportional to $v_{z,0}$, such that the numerator in eq. (21) includes include a term proportional to $\nu$, the vertical epicyclic frequency.

Now taking $z_{d}=(\partial\ln\rho_{0}/\partial z)^{-1}$, substituting in the expression for $v_{z,1}$ and setting $k$=0, eq. (24) becomes

	$\displaystyle 0$	$\displaystyle=$	$\displaystyle(m\Omega-\omega)^{2}\rho_{1}$		(25)
		$\displaystyle-$	$\displaystyle k_{z}^{2}\left(\Phi_{1}\rho_{0}+\sigma^{2}\rho_{1}\right)+ik_{z}\Phi_{1}\rho_{0}\left(\frac{1}{z_{d}}\right)$
		$\displaystyle+$	$\displaystyle\mathcal{A}e^{i(m\Omega-\omega)t}$

where

	$\displaystyle\mathcal{A}$	$\displaystyle=$	$\displaystyle e^{ikr+ik_{z}z}\Big{[}\rho_{0}\nabla^{2}\mathcal{F}+\sigma^{2}\nabla^{2}\mathcal{R}$		(26)
		$\displaystyle+$	$\displaystyle\left(\frac{1}{z_{d}}\right)\rho_{0}\nabla\mathcal{F}$
		$\displaystyle+$	$\displaystyle 2ik_{z}\left(\rho_{0}\nabla\mathcal{F}+\sigma^{2}\nabla\mathcal{R}\right)\Big{]}$

In the non-wave scenario ($k_{z}=0$) the dispersion relation reads

	$\displaystyle 0$	$\displaystyle=$	$\displaystyle(m\Omega-\omega)^{2}\rho_{1}+\Big{[}\rho_{0}\nabla^{2}\mathcal{F}+\sigma^{2}\nabla^{2}\mathcal{R}$		(27)
		$\displaystyle+$	$\displaystyle\left(\frac{1}{z_{d}}\right)\rho_{0}\nabla\mathcal{F}\Big{]}e^{ikr+i(m\Omega-\omega)t}$

whereas in a WKB scenario,

	$\displaystyle 0$	$\displaystyle=$	$\displaystyle(m\Omega-\omega)^{2}\rho_{1}-k_{z}^{2}\left(\Phi_{1}\rho_{0}+\sigma^{2}\rho_{1}\right)$		(28)
		$\displaystyle+$	$\displaystyle ik_{z}\Phi_{1}\rho_{0}\left(\frac{1}{z_{d}}\right)$

The imaginary, out-of-phase term in eq. (28) is notably negligible when $k_{z}z_{d}>>1$. This would be equivalent to the condition required to keep an infinite perturbation consistent with WKB approximation (using the requirement $z_{p}$=$z_{d}$ to keep the perturbation small with respect to the unperturbed disk as $|z|~{}\rightarrow~{}\infty$). Likewise, the second factor in the term in square brackets in eq. (27) drops when $k_{z}z_{d}>>1$, considerably simplifying the expression. However, in the case of the unperturbed Gaussian vertical profile (for which $1/z_{d}=z/h^{2}$), $k_{z}z_{d}>>1$ applies only very near the galactic mid-plane, i.e. where $z<<(k_{z}h)h$. Thus, both the in-phase and out-of-phase terms are relevant for the overall evolution of extended perturbations.

Indeed, in the special case highlighted in the next section in which the perturbation extends to $\pm$ infinity and tracks the decrease in $\rho_{0}$ with increasing $z$ (as in our equilibrium disks) then integration over the vertical direction from $-\infty$ to $\infty$ yields zero when all terms in eqs. (25), (27) and (28) are included. In other words, $\rho_{0}v_{z,1}|_{-\infty}^{\infty}=0$. This is the ‘no mass flux at infinity’ requirement invoked by GLB. As a result, the vertical-only 2D dispersion relation in this ’no mass flux’ scenario reads $(m\Omega-\omega)^{2}=0$ signifying that the vertical direction is neutrally stable to infinite axisymmetric perturbations and stable to all non-axisymmetric perturbations (since $\omega^{2}=m^{2}\Omega_{p}^{2}$ when $k$=0).

This neutral stability characteristic of the vertical direction is leveraged when calculating the 2D dispersion relation later in $\S$ 3, following GLB. Below, it will first be useful to examine how the terms in eqs. (27) and (28) proportional to $1/z_{d}$ contribute to this neutral vertical stability.

In the case of the non-wave perturbation, the third term in eq. (27) dominates away from the mid-plane and

$$(-\omega+m\Omega)^{2}\rho_{1}\approx-\frac{\rho_{0}}{z_{d}}\nabla\Phi_{1}$$

(29)

where $\Phi_{1}=\mathcal{F}$. This corresponds to stability ($\omega^{2}>0$) everywhere since the right-hand side is always positive when $\nabla\Phi_{1}>0$ and $z_{d}<0$, as it is in the equilibrium disks under consideration.

Stablility far above the mid-plane is also a feature of periodic wave perturbations. Consider eq. (28) in a scenario in which the perturbation is extended but finite, for example.⁴⁴4The perturbation must be finite or it will not satisfy the WKB approximation assumed for the present exercise. This requirement is not invoked in other sections unless noted. In the regime $k_{z}z_{d}<<1$, or at heights much larger than $h$ (far away from the mid-plane), the vertical-only continuity equation reads

$$-\frac{\omega^{2}(k^{2}+k_{z}^{2}-T^{2})h^{2}}{f_{g}\nu^{2}k_{z}z}=\textrm{cot}(k_{z}z+kr+(m\Omega-\omega)t)$$

(30)

adopting our Gaussian vertical density profile and letting $\nu^{2}=4\pi G\rho_{0}f_{g}$. Since the arccotangent of the left hand side is $\pm\pi/2$ for all $z>>h$, $\omega$ is always real. It remains real as $z$ approaches nearer to $h$, where the arccotangent of the left hand side is a small positive or negative quantity.

The stability away from the mid-plane is in contrast to the situation very near the mid-plane. In the case of wave perturbations, the dispersion relation where $k_{z}z_{d}>>1$ (and adopting $m$=0 for simplicity) becomes:

$$\omega^{2}=-4\pi G\rho_{0}\frac{k_{z}^{2}}{k_{z}^{2}-T^{2}}+k_{z}^{2}\sigma^{2}$$

(31)

or

$$\omega^{2}=-4\pi G\rho_{0}+k_{z}^{2}\sigma^{2}$$

(32)

in the limit $k_{z}>>T$. (The stability of non-wave perturbations near the mid-plane is examined in $\S$ 2.3.3.) In this situation, perturbations have the opportunity for growth as long as $k_{z}<(4\pi G\rho_{0})/\sigma_{z}^{2}=k_{J}$. In other words, very near to the roughly constant-density mid-plane, instability in the vertical direction proceeds in a Jeans-like manner, unaffected by rotation (Chandrasekhar, 1954) and restricted to similar scales.

As exemplified by the ‘no-mass-flux at infinity’ case described above and considered in detail in $\S\S$ 3.2 and 3.3, the combination of instability near the plane with stability away from the plane results in a neutrally-stable disk.⁵⁵5It is worth noting that the vertical variation in $\omega$ discussed here has not been made explicit in writing eqs. (25), (27) and (28). This corresponds to the assumption that the perturbation’s amplitude and/or phase variations are faster, i.e. $|T|>>(d\omega/dz)t$ or $|k_{z}|>>(d\omega/dz)t$ everwhere. The disk is thus neither as unstable as predicted where $k_{z}z_{d}>>1$ or as stable as predicted where $k_{z}z_{d}<<1$.

Still, eq. (32) does suggest that an avenue to avoid stability would be to perturb the disk over a limited extent, very near the mid-plane. As demonstrated later in $\S$ 3.4, perturbations extending a finite distance above and below $z=0$ are defined by the non-zero mass flux they entail over their extent, with the consequence that the 2D dispersion relation retains terms associated with the vertical direction. As the perturbation’s height decreases, stability approaches the behavior predicted in the limit $k_{z}z_{d}>>1$ near the mid-plane.

Setting $m$=0 and $k_{z}$=0 in eq. (24) and substituting in the expression for $v_{r,1}$, the axisymmetric radial-only dispersion relation reads

$$\omega^{2}=\kappa^{2}-4\pi G\rho_{0}\frac{k^{2}}{k^{2}-T^{2}}+\sigma_{r}^{2}k^{2}\\$$

(33)

or

$$\omega^{2}=\kappa^{2}-4\pi G\rho_{0}+\sigma_{r}^{2}k^{2}\\$$

(34)

in the limit that $k>>T$. This is a restatement of the finding that wave perturbations perpendicular to the axis of rotation (in this case, the radial direction) can be stabilized by rotation, since the scales of instability are pushed over the Jeans length. This was first found by Chandrasekhar (1954) in the case of uniform rotation, then generalized by Bel & Schatzman (1958) for non-uniform rotation (as considered here) and then confirmed to apply in the presence of vertical flattening (Safronov, 1960).

This scenario resembles the case of ‘no vertical mass flux at infinity’ perturbations and 3D perturbations near the mid-plane and so a discussion of the instability scale is postponed until $\S\S$ 2.3.3 and 3.2 . For now it should be noted that simply omitting a vertical perturbation very clearly does not retrieve the 2D Lin-Shu dispersion relation.⁶⁶6As discussed later in $\S$ 3.4.1, to retrieve the Lin-Shu dispersion relation in the long-wavelength limit starting from the 3D dispersion relation requires taking the limit in which the disk is an infinitesimally thin sheet with $\sigma_{z}\rightarrow 0$). Instabilities instead have a Jeans-like quality even in the presence of rotation. (Eq. [34] indeed approaches the condition for Jeans instability in the limit $\kappa\rightarrow 0$.)

Jog & Solomon (1984) pointed out this resemblance to Jeans instability by taking the 2D dispersion relation (see eq. (80), $\S$ 3.4.1) in the small wavelength (large $k$) limit, opposite to the standard long wavelength regime. As examined in $\S$ 2.3.3 (and later in $\S$ 3.4.2), this 3D quality signifies a change in the disk stability threshold compared to the value required for stability to Lin-Shu density-wave perturbations.

This section describes stability and fragmentation from a fully 3D perspective embedded within the gas disk, near to the galactic mid-plane. Later this view is traded for a 2D perspective that can be used to assess the overall stability of the disk (including all material out to $z=\pm\infty$).

Including both radial and vertical perturbations (with $k_{z}\neq 0$ and $k\neq 0$) and substituting in the expression for $v_{r,1}$, equation (24) can be rewritten as

		$\displaystyle\rho_{1}$	$\displaystyle\bigg{[}(m\Omega-\omega)+\left(-\frac{4\pi G\rho_{0}}{k_{\textrm{plane}}^{2}+k_{z}^{2}-T^{2}}+\sigma_{r}^{2}\right)\frac{k^{2}(m\Omega-\omega)}{\Delta}$		(35)
		$\displaystyle+$	$\displaystyle\frac{C_{z}}{(m\Omega-\omega)}\bigg{]}=0,$

where $C_{z}$ represents vertical stability and is equated with either the second term on the right hand side of eq. (27) or the sum of last two terms on the right side of eq. (28), specifically including both in-phase and out-of-phase terms. Notice that when $C_{z}$ is positive (negative) in eqs. (27) or (28) the vertical direction is unstable (stable).

The 3D dispersion relation in eq.(35) is quadratic in $\omega^{2}$, with solutions

$$\omega^{2}=\frac{\omega_{min}^{2}}{2}\left(1\pm\sqrt{1+\frac{4C_{z}\kappa^{2}}{\omega_{min}^{4}}}\right)$$

(36)

where

$$\omega_{min}^{2}=\kappa^{2}+\left(\frac{-4\pi G\rho}{k^{2}+k_{z}^{2}-T^{2}}+\sigma_{r}^{2}\right)k^{2}-C_{z}$$

(37)

in the case that $m$=0. It is straightforward to show that the condition $\omega^{2}$$<$0 can be met when both

$$\omega_{min}^{2}<0$$

(38)

and $C_{z}>0$, corresponding to vertical instability. (When $\omega_{min}^{2}$ is positive, there is a limited range of conditions under which one of two branches of $\omega^{2}$ can still become negative. But we neglect such a scenario here, considering that the criterion in eq. (38) is readily met.)

Notice that when $\omega_{min}^{2}<0$ and $C_{z}>0$, then $\omega_{min}^{2}$ is the minimum that $\omega^{2}$ can reach. In what follows, eq. (38) is used as the condition for instability, with the understanding that growth may happen faster than indicated by $\omega_{min}$. Below, conditions on $k$ (and/or $k_{z}$) for instability specifically near the mid-plane are obtained from eq. (38) in the case of both wave and non-wave 3D perturbations.

For wave perturbations near the mid-plane (i.e. in the limit $k_{z}z_{d}>>1$) that are also assumed to locally satisfy the WKB approximation ($k_{z}>>T$) for illustration purposes, the 3D dispersion relation is written

$$\kappa^{2}+\left(\frac{-4\pi G\rho_{c}}{k^{2}+k_{z}^{2}}+\sigma_{r}^{2}\right)k^{2}+\left(\frac{-4\pi G\rho_{c}}{k^{2}+k_{z}^{2}}+\sigma_{z}^{2}\right)k_{z}^{2}<0.$$

(39)

using that

$$C_{z}=-k_{z}^{2}\left(-\frac{4\pi G\rho_{c}}{k_{\textrm{plane}}^{2}+k_{z}^{2}}+\sigma_{z}^{2}\right)$$

(40)

in this limit (see previous section) with $\rho_{c}=\rho(z\rightarrow 0)$.

Now, setting $k_{S}^{2}=k^{2}+k_{z}^{2}$ (with $S$ denoting ’shell’), instability is found to require

$$\kappa^{2}-4\pi G\rho+\sigma^{2}k_{S}^{2}+k_{z}^{2}(\sigma_{z}^{2}-\sigma_{r}^{2})<0$$

(41)

or

$$k_{S}^{2}<k_{J}^{2}\left(1-\frac{\kappa^{2}}{4\pi G\rho_{c}}\right)$$

(42)

assuming that the velocity dispersion is isotropic ($\sigma_{z}$=$\sigma_{r}$=$\sigma$). Stability within the roughly constant density region staddling the mid-plane thus takes on a Jeans-like quality, though rotation succeeds in increasing the size of stable fragments.

Rotation can also eventually suppress instabilities above a threshold

$$Q_{M}\equiv\kappa^{2}/(4\pi G\rho_{c})>1.$$

(43)

It is notable that the form of this threshold resembles the 3D threshold $\kappa^{2}/(\pi G\rho_{c})\approx 0.3$ determined for the overall disk by GLB better than it matches $Q_{T}$. As discussed in detail later in $\S$ 3.2, the difference in the numerical value of the threshold is a consequence of the vertical extent of the perturbed region.

The threshold $Q_{M}=1$ also corresponds to higher stability threshold than $Q_{T}=1$. In the case of weakly self-gravitating disks (with $\Sigma=\rho_{c}\sqrt{2\pi}h$),

$$Q_{M}=\frac{\pi\alpha^{2}f_{g}Q_{T}^{2}}{8},$$

(44)

while in the case of fully self-gravitating disks (with $\Sigma=\rho_{c}2h$)

$$Q_{M}=\frac{\alpha^{2}f_{g}Q_{T}^{2}}{4}.$$

(45)

Thus, $Q_{M}=1$ is equivalent to $Q_{T}\approx 2$, signifying that disks are more susceptible to partial 3D instability (endemic to the mid-plane) than to total destablization described by the 2D Toomre criterion, as discussed more in $\S$ 4.3.

It is also noteworthy that, as a criterion specifically on the radial $k$ wavenumber, eq. (39) implies

$$k^{2}<k_{J}^{2}\left(1-Q_{M}-k_{z}^{2}h^{2}\right)$$

(46)

with a stability threshold $Q_{M}=1-k_{z}^{2}h^{2}$. Here the radial Jeans length $k_{J}=4\pi G\rho/\sigma_{r}^{2}$. Since $h$ is roughly equivalent to the effective Jeans length (applicable in the presence of thermal and non-thermal motion), it is only when the disk is perturbed on scales larger than the Jeans length that radial fragmentation is seeded. That is, the largest perturbations, with $k_{z}h<<1$, correspond to the highest threshold and thus most easily seed radial fragmentation.

From a more qualitative perspective, the onset of this ‘mid-plane’ Jeans-like instability can be described as follows: At the mid-plane where the density is approximately constant, gas pressure applies a negligible force and only the perturbed pressure force is left to compete with self-gravity. The vertical component of this force is negligible when the wavelength of the perturbation is large, i.e. $k_{z}h<1$ or when the disk is perturbed above the vertical (effective) Jeans length. As a result, the primary competition against self-gravity comes from the pressure force in the plane. Since the pressure force is scale-dependent while self-gravity at the mid-plane is not, the result is that the disk is able to destabilize, but only on scales larger than the radial Jeans length (lengthened by rotation).

It is worth noting that, although this instability is described as occurring ‘at the mid-plane’, it is limited by pressure to scales larger than the vertical Jeans length. Thus the scale height sets the minimum vertical extent of the region that becomes unstable. Indeed, in either the radial or vertical directions, the disk is stabilized by pressure below the Jeans length. According to eq. (46), rotation also contributes to stability on the largest scales.

Non-wave ($k_{z}<<T$) vertical perturbations that are infinite (and satisfy $1/z_{p}=1/z_{d}$) exhibit almost identical behavior near the mid-plane. For these,

$$C_{z}=\rho_{0}\nabla_{z}^{2}\Phi_{1}+\sigma^{2}\nabla_{z}^{2}\rho_{1}+\left(\frac{1}{z_{d}}\right)\left(\nabla_{z}\mathcal{F}+\sigma^{2}\frac{\nabla_{z}\mathcal{R}}{\rho_{0}}\right)$$

(47)

(see previous section). In the limit, $z<<h$ where the unperturbed and perturbed densities are roughly constant and $1/z_{d}=z/h^{2}<<1/h$ (for our adopted Gaussian equilibrium vertical profile), the second and third (pressure) terms drop. Thus, substituting eq. (47) into eq. (35), the condition for instability becomes

$$\kappa^{2}+\frac{\rho_{0}}{\rho_{1}}k^{2}\Phi_{1}+\sigma_{r}^{2}k^{2}-\frac{\rho_{0}}{\rho_{1}}\nabla_{z}^{2}\Phi_{1}<0$$

(48)

or

$$\kappa^{2}-4\pi G\rho_{0}+\sigma_{r}^{2}k^{2}<0$$

(49)

using Poisson’s equation for the substitution $\nabla_{z}^{2}\Phi_{1}=4\pi G\rho_{1}+k^{2}\Phi_{1}$. Specifically near the mid-plane, instability in the presence of an arbitrary infinite perturbation is possible as long as

$$k^{2}<k_{J}^{2}(1-Q_{M})$$

(50)

with stability once again setting in above the threshold $Q_{M}=1$. The minimum instability scale is thus the radial Jeans length in the radial direction and effectively the scale height in the vertical direction (or, more precisely, the extent of the region where the disk density is approximately constant).

As discussed in $\S$ 2.3.1, the factors proportional to $1/z_{d}$ that were neglected here become important away from the mid-plane and serve to lower the threshold for the overall stability of the disk to infinitely extended perturbations. This was previously determined by GLB, who calculated a total (non-wave) perturbation-weighted threshold $\bar{Q}_{M}<Q_{M}=1$ from the 2D dispersion relation derived by integrating over the vertical direction from $-\infty$ to $+\infty$. The next section examines this further, expanding the calculation to include the infinite and finite wave perturbations considered in this work.