ZOMBIE VORTEX INSTABILITY. I. A PURELY HYDRODYNAMIC INSTABILITY TO RESURRECT THE DEAD ZONES OF PROTOPLANETARY DISKS

After 40 years of intense theoretical and computational research, the answer has been a qualified—and unsatisfying—"maybe." It has long been known that molecular viscous diffusion is wholly inadequate to drive accretion because the transport timescale is several orders of magnitude longer than the age of the universe: ${\tau }_{\mathrm{visc}}\sim {R}^{2}/\nu \ \sim \;$ 100,000 Gyr for conditions typical at $R\sim 1$ AU in a protoplanetary disk (PPD) with kinematic viscosity $\nu \sim {v}_{\mathrm{th}}{\ell }$, thermal velocity ${v}_{\mathrm{th}}\sim 1$ km s⁻¹, molecular mean free path ${\ell }\;\sim$ 1 cm. Shakura & Sunyaev (1973) proposed that Keplerian differential rotation may be susceptible to hydrodynamic and/or MHD instabilities that would drive turbulent diffusion, for which they introduced an "eddy" viscosity: ${\nu }_{\mathrm{eddy}}=\alpha {c}_{{\rm{s}}}H$ (the idea being that turbulent "blobs" would move at subsonic speeds over distances of the order of the scale height, and the parameter α would subsume all our ignorance over the details of the transport efficiency). Because the specific angular momentum increases with radius in Keplerian differential rotation, it is often stated that Keplerian shear is stable to purely hydrodynamic instabilities. However, this centrifugal stability criterion (Rayleigh 1917; Synge 1933) is strictly true only for axisymmetric perturbations in an unstratified, inviscid flow, and therefore is of limited relevance for stratified astrophysical disks.

So, what is the nature of hydrodynamic instabilities that may lead to enhanced transport? Historically, by analogy with Poiseuille (pressure-driven) and Couette (viscous-driven) flows in both planar and cylindrical geometries, the go-to mechanism has been subcritical finite-amplitude instabilities that extract energy from the shear at very high Reynolds numbers (Richard & Zahn 1999). This has been, and remains, highly controversial. Balbus et al. (1996) and Hawley et al. (1999) have argued, with theory and computational analyses, that the Coriolis force stabilizes Keplerian shear and shuts off any way (in the absence of magnetic fields) for perturbations to extract energy from the differential rotation in order to power sustained turbulence. Others have argued that even the best to-date numerical simulations lack sufficient resolution to capture nonlinear triggering of instabilities at very high Reynolds numbers (Longaretti 2002; Richard 2003; Dubrulle et al. 2005; Lesur & Longaretti 2005). A related mechanism to generate shear turbulence starts with the transient amplification of a special class of initial conditions, which may trigger nonlinear feedback and continued growth if the amplitude of the transients peak above some critical threshold (Chagelishvili et al. 2003; Tevzadze et al. 2003, 2008; Afshordi et al. 2005; Mukhopadhyay et al. 2005; Umurhan & Shaviv 2005; Salhi & Cambon 2010; Salhi et al. 2013). However, none of these "by-pass to turbulence" approaches regenerate the initial perturbations required to sustain turbulence after the initial burst of turbulence has decayed.

With conflicting results in both theory and numerical work, one may hope that laboratory experiments may provide an independent test. Unfortunately, two recent Taylor–Couette flow experiments (water rotating between two concentric cylinders) have yielded opposite results: experiments by Paoletti & Lathrop (2011) claim a positive measurements of turbulence and outward transport of angular momentum; whereas, a similar set of experiments by Ji et al. (2006), Schartman et al. (2012) claim no such evidence of turbulence nor angular momentum transport. The measurements are difficult and sensitive to the ability of the experimentalists to control secondary flows from their apparatus (e.g., Ekman pumping). Direct numerical simulations of Taylor–Couette flows with laboratory boundary conditions seem to confirm that the axial boundaries are driving the instabilities and enhanced transport (Avila 2012).

Velikhov (1959) and Chandrasekhar (1960, 1961) showed that magnetic fields could destabilize otherwise stable rotational flows. Balbus and Hawley were the first to apply this magnetorotational instability (MRI) to Keplerian shear flows and demonstrate how the resulting turbulence efficiently transported angular momentum outward (Balbus & Hawley 1991; Hawley & Balbus 1991). However, there exist relatively dense, cool and nearly neutral regions in PPDs (approx. 1–10 AU) that likely lack sufficient coupling between matter and magnetic fields (Blaes & Balbus 1994), except perhaps in thin surface layers that have been ionized by cosmic rays or protostellar X-rays (Gammie 1996). The magnetically decoupled regions are known as "dead zones" and there is continued research into whether waves and turbulence from the magnetically coupled MRI-active regions can propagate into the dead zones of PPDs and drive motions that result in angular momentum transport or disrupt the settling of dust (Fleming & Stone 2003; Turner et al. 2007; Turner & Sano 2008; Oishi & Mac Low 2009).

In regions of a laminar near-Keplerian disk that are strongly coupled to magnetic fields, there is no question that of all the known relevant linear instabilities, MRI is the fastest to grow to large amplitude. However, the presence of dead zones in PPDs still motivates researchers to investigate other mechanisms for instability and transport, especially because the dead zone seems to be coincident with planet-forming regions within disks (Matsumura et al. 2007, 2009). A large class of potential instabilities involve global properties of the disk, specifically radial gradients of thermodynamic quantities. Extrema in the radial profile of potential vorticity (a.k.a. vortensity) is unstable to Rossby waves which eventually roll-up into vortices (Lovelace et al. 1999; Li et al. 2000, 2001; Varnière & Tagger 2006; Richard et al. 2013; Meheut et al. 2010, 2012a, 2012b). Radial entropy gradients that are stable according to the Solberg–Høiland criterion can still be unstable to a "subcritical baroclinic instability" (SBI) that is most effective when the thermal cooling time is of order an orbital period (Klahr & Bodenheimer 2003; Petersen et al. 2007a, 2007b; Lesur & Papaloizou 2010; Lyra & Klahr 2011). More recent work on the SBI by Klahr & Hubbard (2014) has shown that the instability is actually linear (i.e., does not require large amplitude perturbations) and more akin to a convective overstability rather than a traditional baroclinic instability. Still others have investigated vertical shear instabilities (VSI) related to the Goldreich–Shubert–Fricke instability in rotating stars (Goldreich & Schubert 1967; Fricke 1968; Nelson et al. 2013). While ZVI needs vertical stable stratification (e.g., entropy increasing in direction opposite of gravity and sufficiently long thermal relaxation times) the VSI needs vertically neutral stratification (e.g., no vertical entropy gradient, or sufficiently short thermal relaxation times). Thus depending on the thermodynamic structure of the disk and its optical properties, one or the other might locally dominate.

Of course, by direct observation, we know that some mechanisms exist within a PPD that leads to the transport of angular momentum in the way that is needed to complete star formation and that allows dust to accumulate and agglomerate into kilometer-sized planetesimals. However, the elucidation of that mechanism or mechanisms cannot be sloughed off as an unimportant "detail." A more realistic treatment of turbulence is essential to quantify the spatial and temporal distribution of dust in a PPD and therefore affects our modeling of the infrared emission from dust (relevant for ALMA and James Webb Space Telescope observations), as well as an initial-mass function (IMF) for planets.

In Barranco & Marcus (2005, hereafter BM05), we investigated the stability of three-dimensional vortices in the midplanes of PPDs. We created a model for a near-equilibrium vortex, which we initialized by hand in our spectral anelastic code that had specially tailored algorithms to handle shear, rotation and stratification (Barranco & Marcus 2006). The vortex suffered a linear instability which resulted in its destruction near the midplane, but we unexpectedly discovered that the stratified regions away from the midplane naturally filled with robust vortices. At the time, we speculated that the initial midplane vortex excited internal gravity waves, which propagated away from the midplane and deposited their energy in stratified regions where the shear rolled vorticity perturbations into new vortices.

However, our original explanation was not entirely complete nor satisfying. Marcus et al. (2013, hereafter MPJH13) correctly diagnosed the true mechanism for the formation of these vortices and in the process identified a new instability, which we now call the "Zombie Vortex Instability" or ZVI. In order to get at the essential nature of the phenomenon, MPJH13 stripped out complicating features of a PPD (e.g., spatially varying gravity and Brunt–Väisälä frequency) and numerically investigated simple Couette flow with constant gravity and constant Brunt–Väisälä frequency in the limit of the Boussinesq approximation. In this far simpler system, MPJH13 found that a small perturbing vortex could trigger an instability in a rapidly rotating, strongly stratified flow, yielding vigorous, space-filling vortices, vortex layers and turbulence.

Figure 1 is a cartoon picture summary of the mechanism outlined in MPJH13. Each panel shows a small patch of a PPD in frame rotating with the gas; upward is the forward azimuthal direction, rightward is increasing radius from the protostar. In the rotating frame, the background flow appears as an anticyclonic shear (as indicated with the arrows). Color illustrates vertical vorticity (aligned with the rotation axis of the disk), with red indicating cyclonic (out of the page) vorticity and blue indicating anticyclonic (into the page) vorticity. The domain is initially filled with Kolmogorov turbulence on small scales (panel (a)). At low Rossby number (when the Coriolis effect is dominant over advection), the turbulence can undergo a "reverse cascade" in which energy transfers to larger scales as small regions of vorticity merge with other regions of vorticity with the same sense of rotation. Regions of vorticity that have the same sense of rotation as the background shear tend to roll-up into compact vortices, whereas vorticity patches with the opposite sense of rotation as the background shear tend to get stretched out into thin layers (Marcus 1993). After a period of reverse cascades and mergers, a region of the disk may only be populated with a few isolated anticyclones that move with the shear and interspersed with cyclonic vortex layers (panel (b)). In our simplified cartoon, we consider a small region with only one remaining anticyclone. This vortex acts as a source of perturbations and waves.

Zoom In Zoom Out Reset image size

The crucial new insight in MPJH13 was the recognition of "baroclinic critical layers" as being sites that are receptive to perturbations. Critical layers are special locations in a shear flow where the coefficients of the highest derivatives of the linearized equations vanish, indicating that the neutrally stable eigenmodes are singular there (Drazin & Reid 1981; Maslowe 1986). Because our focus in this paper is to demonstrate the existence of ZVI, we will defer the mathematics of critical layers to a future paper. For our purposes here, we note that there exist multiple critical layers on either side of the source vortex. Let $(x,y,z)$ be local Cartesianized radial, azimuthal, and vertical directions within a small box rotating with the gas. The cross-stream distance $\delta (m)$ between the source and the mth critical layer is given by (MPJH13):

$$\begin{eqnarray}&&\delta (m)\equiv \displaystyle \frac{N}{| \sigma | }\displaystyle \frac{1}{{k}_{y}}=\displaystyle \frac{N}{| \sigma | }\displaystyle \frac{{L}_{y}/2\pi }{m}\equiv \displaystyle \frac{{\rm{\Delta }}}{m},\end{eqnarray} \tag{ 1.1 }$$

where N is the Brunt–Väisälä frequency, σ is the local shear rate, k_y is a spatial wavenumber associated with a Fourier component of the source vortex, L_y is the domain size in the azimuthal direction, and m is some non-zero integer. The critical layer closest to the source, without intersecting the source itself, has the largest integer m_max. At first glance, it may appear that the closest critical layer depends on the arbitrary choice of the domain size in the azimuthal direction, but in numerical experiments, the physical location of the first excited critical layer is always just outside the perturber, and the value of the integer m_max will scale with L_y. Critical layers for smaller integers $m\leqslant {m}_{\mathrm{max}}$ are also excited. The critical layer farthest from the source generally has m = 1 and is at a distance $\delta (1)\equiv {\rm{\Delta }}$. However, if a boundary of the domain in the cross-stream direction is too close to the source, or if there are physics that prohibit vortex coherence over distances as large Δ, then the excited critical layer farthest from the perturbing source will have $m\gt 1$. In the calculations presented here the boundaries are chosen so that the farthest excited critical layer is 1, and the equations of motion contain no physics that inhibit large vortices from forming, so the farthest distance is Δ. However, as we discuss Section 3.3, the physics of the compressible gas in PPDs will generally limit the size of vortices such that the farthest critical layers do not have m = 1.

While the singular eigenmodes are neutrally stable (do not grow nor decay) to small perturbations, they are susceptible to being forced and excited by the nearby anticyclone. It is important to note that although the cartoon makes it appear that the excited critical layers are at the same height as the perturber, they are in fact located somewhat above and below—the instability is inherently 3D. A dipolar vortex layer (two juxtaposed oppositely signed layers of vorticity) develops at the location of each critical layers (panel (c)). Vortex layers that have the same sense of rotation as the background shear are linearly unstable (panel (d)) and roll-up into compact vortices (Marcus 1993), and so new anticyclones will eventually be spawned at the locations of the critical layers with m vortices at the mth critical layer (panel (e)). The newly spawned vortices from the different critical layers tend to merge into one larger vortex with that vortex located at the critical layer location farthest from the source (generally with m = 1 at a distance Δ). The newly created large vortex then excites its neighboring critical layers (panel (f)) and the process continues. In numerical simulations, one can observe perturbations spawning new vortices farther and farther from the original perturbation vortex.

We named the instability the Zombie Vortex Instability not only because it occurs in the dead zones of PPD, but also because of the way one zombie vortex "infects" its neighboring critical layer, spawning new zombie vortices, which "infect" their neighboring critical layers, etc. In numerical simulations, one can see the entire domain rapidly fill with zombie vortices from one initial vortex. The instability is not an artifact of the numerical method as we have observed it with spectral codes and finite-volume codes (e.g., Athena), with fully compressible, anelastic and Boussinesq treatments of the continuity condition, with and without the shearing box, and with either hyperviscosity or real molecular viscosity. One may ask, if it is so robust, how was it missed in previous numerical calculations? As will be elucidated in this series of papers, ZVI requires vertical stratification, high resolution to resolve the narrow critical layers, a broad spectrum of perturbations (i.e., Kolmogorov, but not Gaussian-peaked), and enough simulation time to allow the critical layers to amplify perturbations.

In this paper, we consider ZVI in the astrophysical context, rather than the Boussinesq approximation appropriate to laboratory Couette flow as we did in MPJH13. In Section 2, we review the hydrodynamic equations of motion and outline the parameter space of interest for ZVI. Because the primary tuning parameter for ZVI is the Brunt–Väisälä frequency N(z), we shall replace linear gravity with a spatially constant value so that N itself is also constant (for vertically isothermal backgrounds). Our previous numerical experiments (including those in BM05) show that ZVI is present even when N(z) is spatially varying. However, those experiments also convinced us that dealing with a constant N makes it much easier to interpret the results and provide a pedagogical explanation for them. In Section 3, we shall study ZVI and its development into turbulence. We shall also investigate triggers of ZVI other than an initial vortex (which was used as the initial condition in BM05 and MPJH13), namely initial isotropic, homogeneous Kolmogorov-like noise. We shall characterize properties of zombie turbulence: its magnitude, its space-filling nature, and its inhomogeneity and anisotropy due to persistent, but turbulent, large anticyclonic vortices and undulating cyclonic vortex layers.

Most of our numerical experiments presented here use the anelastic approximation (which filters acoustic waves) with our pseudo-spectral code that has specially tailored algorithms to handle shear, rotation and stratification (Barranco & Marcus 2006). We have also done simulations with the fully compressible, Godunov finite-volume code Athena (Gardiner & Stone 2008; Stone et al. 2008; Stone & Gardiner 2010), and show that qualitatively and quantitatively we observe similar results. However, our timing tests show that Athena is approximately 100 times more computationally expensive for these numerical experiments when using the same number of grid points as spectral modes in the pseudo-spectral code.

In our discussion in Section 4, we explain why ZVI was not observed in previous studies by others, we review ZVI in the context of other hydrodynamic and MHD instabilities in PPDs, and we consider the role of ZVI in angular momentum transport and planetesimal formation. In a future paper, we will confront ZVI on a more mathematical and theoretical level and address thresholds required to trigger ZVI. In particular, we will show why extraordinarily weak Kolmogorov-like noise is so efficient in exciting ZVI.

Euler's equations (continuity, momentum and energy with no viscous, radiative or diffusive terms) for the gas in the rotating frame are:

$$\begin{eqnarray}&&\displaystyle \frac{\partial \rho }{\partial t}+{\boldsymbol{\nabla }}\cdot (\rho {\boldsymbol{v}})=0,\end{eqnarray} \tag{ 2.2a }$$

$$\begin{eqnarray}\displaystyle \frac{\partial {\boldsymbol{v}}}{\partial t}+({\boldsymbol{v}}\cdot {\boldsymbol{\nabla }}){\boldsymbol{v}} & = & -\displaystyle \frac{1}{\rho }{\boldsymbol{\nabla }}P-g(z)\hat{{\boldsymbol{z}}}-2{{\rm{\Omega }}}_{0}\;\hat{{\boldsymbol{z}}}\times {\boldsymbol{v}}\\ & & +2q{{\rm{\Omega }}}_{0}^{2}\;x\;\hat{{\boldsymbol{x}}},\end{eqnarray} \tag{ 2.2b }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial \epsilon }{\partial t}+({\boldsymbol{v}}\cdot {\boldsymbol{\nabla }})\epsilon =-\displaystyle \frac{P}{\rho }({\boldsymbol{\nabla }}\cdot {\boldsymbol{v}}),\end{eqnarray} \tag{ 2.2c }$$

where ${\boldsymbol{v}}(x,y,z,t)$ is the gas velocity with respect to the center of the box that co-rotates with the PPD. P, ρ, and are, respectively, the gas pressure, density and specific internal (thermal) energy. The vertical acceleration of gravity $-g(z)$ does not include any contribution from the self-gravity of the gas itself. The term $-2{{\rm{\Omega }}}_{0}\;\hat{{\boldsymbol{z}}}\times {\boldsymbol{v}}$ is the Coriolis force, and the term $2q{{\rm{\Omega }}}_{0}^{2}\;x\hat{{\boldsymbol{x}}}$ is the tidal acceleration that arises from the difference between the inward central force that causes circular motion $R{{\rm{\Omega }}}^{2}(R)$ and the outward centrifugal force $R{{\rm{\Omega }}}_{0}^{2}$ due to being in a rotating frame. The term on the right-hand side of the thermal energy equation is the rate at which pressure forces do work due to volume changes of fluid parcels. To close the system of equations, we assume an ideal gas for the equation of state:

$$\begin{eqnarray}&&P=(\gamma -1)\rho \epsilon =\rho {\mathcal{R}}T,\end{eqnarray} \tag{ 2.3 }$$

where T is the gas temperature, ${\mathcal{R}}$ is the gas constant, and $\gamma \equiv {C}_{{\rm{P}}}/{C}_{{\rm{V}}}$ is the adiabatic index, or ratio of specific heats at constant pressure and constant volume. Throughout this paper, we set ratio of specific heats $\gamma =5/3$.

A time-independent equilibrium (denoted with overbars) to the above set of equations is:

$$\begin{eqnarray}&&\bar{{\boldsymbol{v}}}=\left({\bar{v}}_{x},{\bar{v}}_{y},{\bar{v}}_{z}\right)=\left(0,-q{{\rm{\Omega }}}_{0}\;x,0\right),\end{eqnarray} \tag{ 2.4a }$$

$$\begin{eqnarray}&&\displaystyle \frac{d\bar{P}(z)}{{dz}}=-\bar{\rho }(z)g(z).\end{eqnarray} \tag{ 2.4b }$$

Because there is no radiative or diffusive effects in Equations (2.2), there is a degeneracy of allowable equilibrium thermodynamic profiles, so the equilibrium temperature $\bar{T}(z)$ corresponding to the equilibrium velocity $\bar{{\boldsymbol{v}}}$ is not unique. In fact, one thermodynamic quantity, $\bar{P}(z)$, $\bar{T}(z)$, $\bar{\rho }(z)$, or $\bar{\epsilon }(z)$ can be arbitrarily specified, and the others will follow. In this work, we shall assume $\bar{T}(z)={T}_{0}=\mathrm{constant}$, so:

$$\begin{eqnarray}&&\displaystyle \frac{\bar{P}(z)}{{\mathcal{R}}{T}_{0}}=\bar{\rho }(z)={\rho }_{0}\mathrm{exp}\left[-{\displaystyle \int }_{0}^{z}\displaystyle \frac{g(z\prime )}{{\mathcal{R}}{T}_{0}}\;{dz}\prime \right],\end{eqnarray} \tag{ 2.5 }$$

where ${\rho }_{0}$ is the midplane density and z = 0 is the midplane of the disk.

The anelastic approximation has been extensively used in the study of deep, subsonic convection in planetary atmospheres (Ogura & Phillips 1962; Gough 1969; Bannon 1996) and stars (Gilman & Glatzmaier 1981; Glatzmaier & Gilman 1981a, 1981b). We have previously used the anelastic approximation to study three-dimensional vortices in PPDs (Barranco et al. 2000; Barranco & Marcus 2005, 2006) and the Kelvin–Helmholtz instability of settled dust layers in PPDs (Barranco 2009; Lee et al. 2010a, 2010b). The basic idea is that there may be large variations in the background pressure and density in hydrostatic equilibrium, but that at any height in the atmosphere, the fluctuations of the pressure and density are small compared to the background values at that height. Mathematically, all flow variables are expanded in powers of the Mach number (presumed to be small), and only lowest order terms are retained in the equations of motion. In general, the anelastic approximation results in a valid description of the flow when:

$$\begin{eqnarray}&&\displaystyle \frac{| {\boldsymbol{v}}(x,y,z,t)| }{{c}_{{\rm{s}}}}\equiv {Ma}\ll 1,\end{eqnarray} \tag{ 2.6a }$$

$$\begin{eqnarray}\displaystyle \frac{| \rho (x,y,z,t)-\bar{\rho }(z)| }{\bar{\rho }(z)} & \sim & \displaystyle \frac{| P(x,y,z,t)-\bar{P}(z)| }{\bar{P}(z)}\\ & \sim & {{Ma}}^{2}\ll 1,\end{eqnarray} \tag{ 2.6b }$$

where c_s is the local sound speed. Diagnostically, it is trivial to check every step during a numerical calculation whether the conditions (2.6) are satisfied.

There are a number of variations of the anelastic approximation, all equivalent to the same order of the Mach number, but varying in whether they retain or drop certain smaller terms or in exactly how they linearize the equation of state, often with an eye toward conserving some quantity such as energy; see references above as well as Brown et al. (2012), Vasil et al. (2013). For the special case of uniform temperature equilibrium background, the most straight-forward application of the anelastic approximation yields energy-conserving equations, so we shall defer going into more detail until future papers in which we consider different background temperature profiles. The Euler equations with the anelastic approximation for a vertically uniform background temperature are:

$$\begin{eqnarray}&&{\boldsymbol{\nabla }}\cdot [\bar{\rho }(z){\boldsymbol{v}}]=0,\end{eqnarray} \tag{ 2.7a }$$

$$\begin{eqnarray}\displaystyle \frac{\partial {\boldsymbol{v}}}{\partial t} & = & -({\boldsymbol{v}}\cdot {\boldsymbol{\nabla }}){\boldsymbol{v}}-2{{\rm{\Omega }}}_{0}\;\hat{{\boldsymbol{z}}}\times {\boldsymbol{v}}+2q{{\rm{\Omega }}}_{0}^{2}\;x\;\hat{{\boldsymbol{x}}}\\ & & -{\boldsymbol{\nabla }}\left[\displaystyle \frac{P-\bar{P}(z)}{\bar{\rho }(z)}\right]+\left[\displaystyle \frac{T-{T}_{0}}{{T}_{0}}\right]g(z)\;\hat{{\boldsymbol{z}}},\end{eqnarray} \tag{ 2.7b }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial T}{\partial t}=-({\boldsymbol{v}}\cdot {\boldsymbol{\nabla }})T-\left[\displaystyle \frac{{N}^{2}(z)}{g(z)}\right]{{Tv}}_{z},\end{eqnarray} \tag{ 2.7c }$$

$$\begin{eqnarray}&&\left[\displaystyle \frac{P-\bar{P}(z)}{\bar{P}(z)}\right]=\left[\displaystyle \frac{\rho -\bar{\rho }(z)}{\bar{\rho }(z)}\right]+\left[\displaystyle \frac{T-{T}_{0}}{{T}_{0}}\right],\end{eqnarray} \tag{ 2.7d }$$

where N(z) is the Brunt–Väisälä frequency, a measure of the convective stability of the atmosphere. In general, ${N}^{2}\propto {\boldsymbol{g}}\cdot {\boldsymbol{\nabla }}\bar{s}$, where $\bar{s}$ is the background entropy. Here, we focus on vertical entropy gradients, so:

$$\begin{eqnarray}&&{N}^{2}(z)\equiv (g/{C}_{{\rm{P}}})(d\bar{s}/{dz}),\end{eqnarray} \tag{ 2.8a }$$

$$\begin{eqnarray}&&={g}^{2}/({C}_{{\rm{P}}}{T}_{0})(\mathrm{for}\;\mathrm{an}\;\mathrm{isothermal}\;\mathrm{hydrostatic}\;\mathrm{background}).\end{eqnarray} \tag{ 2.8b }$$

One of the main dynamical consequences of this approximation is that the total density is replaced by the time-independent mean density in the mass continuity equation, which has the effect of filtering acoustic waves, but allowing slower wave phenomena such as internal gravity waves, Rossby waves, etc.

Similarly to the fully compressible equations, the dissipationless anelastic equations contain a degeneracy, so any one thermodynamic quantity, $\bar{P}(z)$, $\bar{T}(z)$, $\bar{\rho }(z)$, or $\bar{\epsilon }(z)$ can be arbitrarily specified. As before, we arbitrarily set the background equilibrium temperature to be constant, $\bar{T}(z)={T}_{0}$. We note that the isothermal equilibrium solutions to the fully compressible equations are the same isothermal equilibrium solutions to the anelastic equations.

In a PPD, the vertical component of the protostellar gravity varies linearly with distance from the midplane, which yields Gaussian profiles of pressure and density for isothermal vertical profiles, and a linear profile for the Brunt–Väisälä frequency (see BM05):

$$\begin{eqnarray}&&g(z)={{\rm{\Omega }}}_{0}^{2}\;z\quad \mathrm{for}\;| z| \ll {R}_{0},\end{eqnarray} \tag{ 2.9a }$$

$$\begin{eqnarray}\bar{P}(z)/({\mathcal{R}}{T}_{0}) & = & \bar{\rho }(z)={\rho }_{0}\mathrm{exp}(-{z}^{2}/2\;{H}^{2}),\\ H\equiv \sqrt{{\mathcal{R}}{T}_{0}}/{{\rm{\Omega }}}_{0} & = & {c}_{{\rm{s}}}/{{\rm{\Omega }}}_{0},\end{eqnarray} \tag{ 2.9b }$$

$$\begin{eqnarray}&&N(z)=\sqrt{{\mathcal{R}}/{C}_{{\rm{P}}}}\;{{\rm{\Omega }}}_{0}^{2}| z| /{c}_{{\rm{s}}}=\sqrt{(\gamma -1)/\gamma }\;{{\rm{\Omega }}}_{0}| z| /H,\end{eqnarray} \tag{ 2.9c }$$

where ${c}_{{\rm{s}}}\equiv \sqrt{{\mathcal{R}}{T}_{0}}$ is the isothermal sound speed.⁵

As we shall demonstrate, the essential tuning parameter for the ZVI is the Brunt–Väisälä frequency, so for this paper, rather than a vertical gravity that is linear in $| z|$, we use a constant vertical gravity, which yields exponential profiles of pressure and density for isothermal vertical profiles, and a constant Brunt–Väisälä frequency N₀:

$$\begin{eqnarray}&&g={g}_{0},\end{eqnarray} \tag{ 2.10a }$$

$$\begin{eqnarray}\bar{P}(z)/({\mathcal{R}}{T}_{0}) & = & \bar{\rho }(z)={\rho }_{0}\mathrm{exp}(-z/{H}_{0}),\\ {H}_{0}\equiv {\mathcal{R}}{T}_{0}/{g}_{0} & = & {c}_{{\rm{s}}}^{2}/{g}_{0},\end{eqnarray} \tag{ 2.10b }$$

$$\begin{eqnarray}&&{N}_{0}=\sqrt{{\mathcal{R}}/{C}_{{\rm{P}}}}\;{g}_{0}/{c}_{{\rm{s}}}=\sqrt{[(\gamma -1)/\gamma ]({g}_{0}/{H}_{0})}.\end{eqnarray} \tag{ 2.10c }$$

In BM05, we previously showed that ZVI is present even when N(z) is spatially varying. In this paper, by focusing on constant N₀, we will be able to probe what thresholds of stable stratification are required to trigger ZVI. However, this restriction prevents us from exploring how gravity waves and turbulence generated by ZVI will propagate into the nearly unstratified midplane; we defer that important issue to future work.

All of the simulations in this paper are periodic in the azimuthal or stream-wise direction y, and use shearing box boundary conditions in the radial or cross-stream direction x (Goldreich & Lynden-Bell 1965; Marcus & Press 1977; Rogallo 1981). In the pseudo-spectral calculations, we have performed computations in which the top and bottom boundary conditions are treated by (a) mapping the boundaries to $z\to \pm \infty$ (as in BM05), (b) setting the vertical velocity to zero (i.e., rigid boundary conditions), or (c) forcing all variables to be spatially periodic in z. When choice (c) is used, it is necessary to include a small artificial damping layer adjacent to the top and bottom boundaries to smooth out the inherent non-periodicity of $\bar{\rho }(z)$ and to prevent spurious reflections of inertio-gravity waves and other unphysical effects. All three of these boundary conditions produce quantitatively similar solutions. However, option (c) leads to a triply periodic and computationally fast code, and thus was used for exploring parameter space for almost all the pseudo-spectral simulations reported in this paper. The pseudo-spectral domain size was ${H}_{0}^{3}$ and was resolved with 256³ spectral modes. The algorithms used in the pseudo-spectral calculations have no inherent dissipations, so an explicit, weak hyperviscosity is required.

Athena is a Godunov finite-volume code that can solve the fully compressible hydrodynamic equations in a shearing box (Gardiner & Stone 2008; Stone et al. 2008; Stone & Gardiner 2010). The Athena simulations used a computational domain of size ${H}_{0}^{3}$ resolved with 256³ grid points. The vertical velocity was set to zero at the top and bottom boundaries. Note that the Athena code has no explicit dissipation, but its finite-volume algorithms are inherently dissipative. We found that lower resolution simulations with Athena (e.g., 128³ grid points) did not develop ZVI.

There are three (reciprocal) timescales of interest: the rotation rate ${{\rm{\Omega }}}_{0}$, the shear rate $\sigma \equiv -q{{\rm{\Omega }}}_{0}$, and the Brunt–Väisälä frequency N₀. We consider only Keplerian shear, q = 1.5 (with the exception of the calculations reported in Figure 2). We define the dimensionless ratio $\beta \equiv {N}_{0}/{{\rm{\Omega }}}_{0}$. For a PPD with linear vertical gravity and no vertical temperature gradient, $\beta (z)=0.63| z| /{H}_{0}$. In BM05, we observed that vortices naturally formed in stratified regions above and below the midplane: $1\lesssim | z| /{H}_{0}\lesssim 4$, corresponding to $0.63\lesssim \beta \lesssim 2.5$. For constant gravity g₀, one can show from Equation (2.10) that:

$$\begin{eqnarray}&&{c}_{{\rm{s}}}=\sqrt{\gamma /(\gamma -1)}\;\beta \;{H}_{0}\;{{\rm{\Omega }}}_{0}=\sqrt{5/2}\;\beta \;{H}_{0}\;{{\rm{\Omega }}}_{0}.\end{eqnarray} \tag{ 2.11 }$$

With constant gravity, the only independent length scale in either the fully compressible or anelastic equations is the gas pressure scale height H₀. Therefore, the size of the computational domain $({L}_{x},{L}_{y},{L}_{z})$ introduces three more dimensionless parameters. In this study, we set ${L}_{x}/{H}_{0}={L}_{y}/{H}_{0}={L}_{z}/{H}_{0}=1$.

Zoom In Zoom Out Reset image size

One of the most cited pieces of evidence that PPDs are stable to purely hydrodynamic instabilities is given by Figure 1 in Balbus et al. (1996), which we copy below as Figure 2(a). The figure shows the growth or decay of kinetic energy fluctuations as a function of time for various values of q, where q is defined in Equation (2.1) and where the kinetic energy fluctuation is defined as $| {\boldsymbol{v}}(x,y,z,t)-{\bar{v}}_{y}\;\hat{{\boldsymbol{y}}}{| }^{2}/2$, where ${\bar{v}}_{y}\;\hat{{\boldsymbol{y}}}$ is the steady equilibrium flow in Equations (2.4). The initial-value calculations that produced Figure 2(a) used the fully compressible Equations (2.1)–(2.3), and were initialized with the steady equilibrium velocity in Equation (2.4a). However, this flow was computed with no vertical gravity and was therefore initialized with the equilibrium that had $\bar{P}={P}_{0}$, $\bar{\rho }={\rho }_{0}$, and $\bar{T}={T}_{0}$. The initial flow was perturbed with small-amplitude noise (with rms velocity of 0.1) that had a Gaussian-like energy spectrum as a function of wavenumber. The equations were solved with periodic boundary conditions in y and z, and with shearing box boundaries in x. Because the energy fluctuations all increase in time for $q\geqslant 2$ and decrease in time for $q\lt 2$, Figure 2(a) is used to support the hypothesis that these flows of ideal gases are stable (unstable) to all purely hydrodynamic perturbations when the angular momentum per unit mass of the flow, ${R}^{2}\;{\rm{\Omega }}(R)\propto {R}^{(2-q)}$, increases (decreases) in the radially outward direction. This hypothesis is consistent with Rayleigh's centrifugal stability theorem; however, it must be noted that Rayleigh's theorem was proved only for the case of constant density fluids and not for stratified ideal gases (Rayleigh 1917; Synge 1933), and therefore, it may not be applicable to astrophysical flows in disks. In particular, initial conditions in which the density is constant prohibit baroclinic instabilities such as ZVI (MPJH13). The curve labeled with "shr" in Figure 2(a) corresponds to the case of pure Cartesian shear with the Coriolis and tidal acceleration terms dropped from Equation (2.2b), which would be appropriate for fully compressible flow in a channel with cross-stream shear, but no rotation.

In Figure 2(b), we reproduce the results of Figure 2(a) with our pseudo-spectral numerical code (Barranco & Marcus 2006) using the anelastic Equations (2.7) with no vertical gravity. The initial noise that we added to the equilibrium velocity in all of our numerical experiments was random, homogeneous, and isotropic with a Kolmogorov spectrum so that the kinetic energy per unit mass was $E(k)\propto {k}^{-5/3}$, where k is the wavenumber. Our primary reason for reproducing this figure from Balbus et al. (1996) is to show that we get the same stability results for no vertical gravity with a pseudo-spectral code (versus finite-volume code) which solve the anelastic equations (versus the fully compressible equations).

How does vertical stratification affect stability? The decaying dashed curve in Figure 3(a) is the time evolution of the kinetic energy fluctuations for a unstratified flow similar to the one as given by the curve in Figure 2(b) with q = 1.5 and $g=N=\beta =0$, but plotted for a longer time and a smaller value of initial noise. The solid curve in Figure 3(a) shows the evolution of kinetic energy fluctuations when vertical gravity is included; here, $g={g}_{0}\ne 0$ is constant, $N={N}_{0}\ne 0$, and $\beta =2$. The growth of the kinetic energy fluctuations shows that the flow is unstable with q = 1.5 when vertical gravity is present. Figure 3(b) shows the kinetic energy growth for the same flow as shown with the solid curve as in Figure 3(a), but computed with the fully compressible equations with Athena. The rapid growth of the kinetic energy fluctuations shows that ZVI is not an artifact of the anelastic approximation or spectral codes.

The growth of kinetic energy fluctuations in Figure 3 is evidence of instability, but spatial plots of the relative vorticity, defined as $\omega \equiv {\boldsymbol{\nabla }}\times ({\boldsymbol{v}}-\bar{{\boldsymbol{v}}})={\boldsymbol{\nabla }}\times {\boldsymbol{v}}+q\;{{\rm{\Omega }}}_{0}\;\hat{{\boldsymbol{z}}}$, are more useful in illustrating ZVI and its nonlinear evolution. The point-wise dimensionless Rossby number, defined as the ratio of relative vorticity to the global Keplerian vorticity $\mathrm{Ro}(x,y,z,t)\equiv [\hat{{\boldsymbol{z}}}\cdot \omega (x,y,z,t)]/(2{{\rm{\Omega }}}_{0})$, is a better diagnostic of late-time zombie turbulence than, say, the Mach number, because the amplitude of the Mach number at late-times is sensitive to the values of ${L}_{y}/{L}_{x}$ and ${H}_{0}/{L}_{x}$ (which are chosen arbitrarily in our simulations), whereas the Rossby number is not—see Section 3.3. The simulations in this section are initialized by random perturbations which are closer to the expected conditions in a PPD rather than by an ad hoc isolated initial vortex as was done in the simulations in BM05 and MPJH13. Figure 4 shows $\mathrm{Ro}(x,y,z,t)$ for an anelastic flow in the $x-z$ plane at four different times and at an arbitrary stream-wise location in y. (Because the equations, boundary conditions and initial conditions are invariant under translation in y, the flow in all $x-z$ planes is statistically the same for all time.) Figure 5 shows $\mathrm{Ro}(x,y,z,t)$ for the same flow in the $x-y$ plane at z = 0, which is midway between the upper and lower boundaries. The parameter values, initial conditions and boundary conditions of the flow in Figures 4 and 5 are similar to the flow shown by the solid curve in Figure 3(a) with the difference that $\beta =1$, rather than 2. The initial perturbing noise has a power-law spectrum with random velocity phases, so there are no inherent length scales or coherent features of any kind in the initial perturbation. For Kolmogorov noise, the energy is dominated by Fourier modes with the smallest spatial wave numbers, while the relative vorticity field is dominated by Fourier modes with the largest spatial wavenumbers. Thus, Figures 4(a) and 5(a) are dominated by the smallest length scales (which are equal to the spatial resolution of the calculation, ${L}_{x}/256$ in each direction). Much of the initial vorticity is quickly stretched and elongated in the stream-wise direction by the Keplerian shear, as shown in Figures 4(b) and 5(b). The stretching cascades the initial kinetic energy to spatial wavenumbers with higher Fourier modes where it is quickly dissipated by hyperviscosity (similar to Figure 3(a) at $t\simeq 50$). In competition with the stretching, there is a weaker reverse cascade of energy from small to large length scales due to the fact the fluid is rapidly rotating. Here, the reverse cascade manifests itself via vortex mergers. In a rapidly rotating flow without a background shear, like-signed vortices of both signs merge, but when a strong background shear is also present, vortices with the same sign as the background shear merge into larger vortices with shapes that are elongated in the stream-wise y direction, while vortices with a sign opposite to the background shear flow are stretched into thin stream-wise-oriented filaments and layers (Marcus 1993). The shear in a PPD is anticyclonic, and therefore anticyclones merge, while cyclones are stretched.⁶ The result of the competition between the forward energy cascade (resulting in energy dissipation) and the reverse energy cascade (resulting in the mergers of anticyclones) is visible in Figures 4(c) and 5(c) at time t = 51 years; the flow is dominated by a few scattered anticyclones. Incidentally, this flow looks very much like the initial condition of a single isolated anticyclone that we studied in MPJH13 and illustrated in Figure 1.

By 51 years, ZVI is well underway in Figures 4 and 5. The scattered anticyclones which emerged from the reverse cascade have now excited neighboring critical layers to form dipolar vortex layers. The critical layers themselves are easy to identify in plots of ${v}_{z}(x,y,z,t)$ (not shown), but difficult to observe in plots of $\mathrm{Ro}(x,y,z,t)$. However, the vortex layers that the critical layers create are easily identifiable in plots of $\mathrm{Ro}(x,y,z,t)$ because each dipolar sheet appears as a pair of opposite-signed "stripes" in the x–y plane (with the critical layer sandwiched between them—Figure 1(c)). In an initial-value simulation with an initial perturbation consisting of noise, it is sometimes difficult to see the dipolar vortex layers due to the fast instability of the anticyclonic sheet. However in carefully controlled simulations consisting of a single perturbing vortex, rather than noise, dipolar layers can be found easily, as shown in Figure 1(a) in MPJH13. As discussed in Section 1.2, the layers of cyclonic vorticity aligned in the stream-wise direction are linearly stable. In contrast, layers of anticyclonic vorticity aligned in the stream-wise direction are linearly unstable and roll-up into stable anticyclones.

At late times (Figures 4(d) and 5(d)), the flow has reached a statistically steady state of zombie turbulence. One can see by eye, especially in Figure 5(d), that the flow has formed a pattern with cross-stream or x wavenumber of ∼7. The relative vorticity, although very turbulent, has developed spatial coherence and is neither homogeneous nor isotropic. The cyclonic vorticity has formed approximately two-dimensional layers that are approximately aligned in the y–z planes. Between these planes are approximately ellipsoidal turbulent anticyclones. At intermediate times in our calculation, we would expect a dominant x-wavenumber of ${L}_{x}/{\rm{\Delta }}=3\pi \simeq 10$, and at times between that of Figure 5(c) and that of Figure 5(d), the dominant x-wavenumber is indeed ∼10. However, at later times, the vortices grow in size, run into each other in the x direction, and become larger. Due to mergers, the late-time x-diameters of the vortices and the average spacing in x between the cyclonic sheers both become slightly larger than Δ. Therefore, the dominant wavenumber in x becomes slightly less than ${L}_{x}/{\rm{\Delta }}$. We have carried out initial-value calculations of zombie turbulence in anelastic flows and Boussinesq flows (MPJH13) for a wide variety of parameters, and the $\mathrm{Ro}(x,y,z,t)$ always look like the pattern in Figures 4(d) and 5(d). Our one fully compressible calculation (Figure 6) also shows this same pattern at late times. In Figure 6, $\beta =2$ and L_y = L_x, so we would expect that ${L}_{x}/{\rm{\Delta }}\simeq 5$, which argues that the x wavenumber that initially dominates the flow is 5. This wavenumber agrees with this simulation at intermediate times. At late times, the wavenumber in Figure 6 has been reduced to ∼4 by vortex mergers and reverse energy cascades. Note that the two flows shown in Figures 8 and 9, also have $\beta =2$ and L_y = L_x, and although we have not included figures of their $\mathrm{Ro}(x,y,z,t)$, these flows also have dominant x-wavenumbers of ∼5 at intermediate times and ∼4 at late times. Thus, the main signature of zombie turbulence that differentiates it from other types of turbulence in shear flows that are triggered by finite-amplitude instabilities is that the turbulence in ZVI is never homogeneous or isotropic; it is dominated by cyclonic layers and elongated anticyclones aligned in the stream-wise direction of the background shear, and the pattern of layers and anticyclones has a dominant x wavenumber slightly less than ${L}_{x}/{\rm{\Delta }}$.

The aspect ratio χ (defined as the stream-wise diameter of a vortex in the y direction divided by its cross-stream diameter in x) of the zombie vortices is approximately the same as the laminar vortices studied by Moore & Saffman (1971) and investigated in BM05 in the PPD context, where χ is given implicitly by

$$\begin{eqnarray}&&-\displaystyle \frac{{\omega }_{z}}{q\;{{\rm{\Omega }}}_{0}}\equiv \displaystyle \frac{-2{Ro}}{q}=\displaystyle \frac{(\chi +1)}{(\chi -1)}\displaystyle \frac{1}{\chi }.\end{eqnarray} \tag{ 3.1 }$$

The Moore–Saffman relation was derived for a steady two-dimensional vortex with uniform relative vorticity embedded in flow with uniform shear. However, the turbulent zombie vortices have aspect ratios similar to that in Equation (3.1) because the relation is the quantification of the fact that a large relative vorticity tends to make a vortex "round" and a large background shear tends to elongate a vortex is its stream-wise direction. At late times the characteristic magnitude of Ro of the anticyclones in zombie turbulence is always ∼−0.25, so regardless of the parameters of the flow, Equation (3.1) shows that in a Keplerian disk, the aspect ratios χ of zombie vortices are between 4 and 5.

Figures 7–10 show how the rms Rossby numbers Ro_rms(t) and rms Mach number ${\mathrm{Ma}}_{\mathrm{rms}}(t)$ evolve in time for three anelastic flows and one fully compressible flow. The flows in Figures 8 and 9 have the same flow parameters, with the only difference being the rms value of their initial noise. All of the statistical properties, including Ro_rms(t) and ${\mathrm{Ma}}_{\mathrm{rms}}(t)$, of the late-time flows in Figures 8 and 9 are nearly the same, which suggests that they evolve to a common attracting solution that is independent of the initial conditions. For all four flows in Figures 7–10, the values of Ro_rms(t) initially plummet due to small-scale hyperviscosity (or numerical diffusion in Athena), but then grow after ZVI develops, and eventually reach a plateau with Ro_rms(t) between 0.2 and 0.3.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The Ro_rms(t) is the ratio of the strength of the vorticity in the turbulence to that of the unperturbed Keplerian flow, so $| {{Ro}}_{\mathrm{rms}}(t)| \sim 0.2-0.3$ indicates large amplitude turbulence. We have no rigorous explanation for why the Rossby numbers evolve to these values in all of our simulations that produce zombie turbulence, regardless of parameter values. However, this range of values of $| {{Ro}}_{\mathrm{rms}}|$ is consistent with the picture that we presented in MPJH13 and Figure 1 of how ZVI spreads throughout the computational domain after it starts at one or more discrete locations. We showed that after the first critical layer is excited, it creates a vortex dipolar layer that rolls up into the first zombie vortex. Then, that first zombie vortex triggers an instability in an adjacent critical layer in an unperturbed region of the flow, which then launches a new vortex layer and zombie vortex that triggers the instability in its adjacent critical layer, and so on ad infinitum. However, this scenario requires that the vortices that make up the zombie turbulence have sufficiently large amplitudes to trigger the finite-amplitude ZVI in unperturbed regions of the flow.

It is more common among astrophysicists to measure the amplitude of turbulence with the Mach number, rather than the Rossby number, so it is somewhat disconcerting that the late-time values of ${\mathrm{Ma}}_{\mathrm{rms}}(t)$ in our calculations depend on the values of ${H}_{0}/{L}_{y}$, and that the rms Mach numbers can be made big or small by adjusting the size of the computational box. It is also discouraging (to those who might believe that zombie turbulence can transport angular momentum in PPDs) that ${\mathrm{Ma}}_{\mathrm{rms}}(t)$ is only ∼10⁻² in our calculations. We now show that the dependence of ${\mathrm{Ma}}_{\mathrm{rms}}$ on the size of the computational box is easily explained, and that in a PPD we expect ${\mathrm{Ma}}_{\mathrm{rms}}$ to be equal to $| {{Ro}}_{\mathrm{rms}}|$ at late times. Note that the rms Rossby number is approximately

$$\begin{eqnarray}&&{{Ro}}_{\mathrm{rms}}\simeq {V}_{\mathrm{eddy}}({L}_{\omega })/({L}_{\omega }\;{{\rm{\Omega }}}_{0}),\end{eqnarray} \tag{ 3.2 }$$

where ${L}_{\omega }$ is the characteristic length scale of the flow where the vorticity has its maximum value and V_eddy(L) is the characteristic velocity of a turbulent eddy of diameter L.⁷ The rms Mach number is approximately

$$\begin{eqnarray}&&{\mathrm{Ma}}_{\mathrm{rms}}\simeq {V}_{\mathrm{eddy}}({L}_{v})/{c}_{{\rm{s}}},\end{eqnarray} \tag{ 3.3 }$$

where L_v is the characteristic length scale of the flow where the velocity has its maximum value. In a Kolmogorov energy spectrum, ${L}_{v}\gg {L}_{\omega }$ because the velocity is dominated by the large length scales and the vorticity by the small length scales. However, zombie turbulence does not have a Kolmogorov spectrum and is dominated by large vortices. In flows dominated by large vortices, L_v and ${L}_{\omega }$ are both approximately the mean size of the vortices. Therefore, in a flow dominated by zombie vortices,

$$\begin{eqnarray}&&{\mathrm{Ma}}_{\mathrm{rms}}\simeq ({{\rm{\Omega }}}_{0}{L}_{{zv}}/{c}_{{\rm{s}}})\;{{Ro}}_{\mathrm{rms}},\end{eqnarray} \tag{ 3.4 }$$

where L_zv is the characteristic diameter of a zombie vortex. In a flow with constant vertical gravity (either an anelastic or a fully compressible flow with $\gamma =5/3$), Equation (2.11) determines c_s, and Equation (3.4) then shows that a flow with zombie vortices has

$$\begin{eqnarray}&&{\mathrm{Ma}}_{\mathrm{rms}}\simeq \sqrt{2/5}\;{{Ro}}_{\mathrm{rms}}\;{\beta }^{-1}\;({L}_{{zv}}/{H}_{0}).\end{eqnarray} \tag{ 3.5 }$$

In all of the calculations presented in this paper, ${L}_{{zv}}\simeq {\rm{\Delta }}$, so using Equation (1.1) for Δ, Equation (3.5) becomes

$$\begin{eqnarray}&&{\mathrm{Ma}}_{\mathrm{rms}}\simeq \displaystyle \frac{\sqrt{2/5}}{3\pi }\;\displaystyle \frac{{L}_{y}}{{H}_{0}}\;{{Ro}}_{\mathrm{rms}}\;\simeq 0.067\;{{Ro}}_{\mathrm{rms}},\end{eqnarray} \tag{ 3.6 }$$

where we used the fact that in our computations ${L}_{y}={H}_{0}$. Equation (3.6) shows that ${\mathrm{Ma}}_{\mathrm{rms}}$ depends on the size of the computational box, is consistent with Figures 7–9 at late times, and explains why our calculations produce zombie turbulence with ${\mathrm{Ma}}_{\mathrm{rms}}\simeq {10}^{-2}$.⁸

Equations (3.5) and (3.6) do not apply to a PPD for two reasons. First, due to the fact that the vertical component of the protostellar gravity is linear in z, rather than constant, the sound speed is given by Equation (2.9b) rather than Equation (2.11). Therefore, Equation (3.5) is replaced by

$$\begin{eqnarray}&&{\mathrm{Ma}}_{\mathrm{rms}}\simeq ({L}_{{zv}}/H)\;{{Ro}}_{\mathrm{rms}}.\end{eqnarray} \tag{ 3.7 }$$

The second reason that Equation (3.6) does not apply to a PPD is we expect ${L}_{{zv}}\simeq H$ for zombie vortices in a PPD, rather than ${L}_{{zv}}\simeq {\rm{\Delta }}$. This expectation follows from the observation that in our numerical simulations, zombie vortices grow until they run into each other in the x direction (or they run into the x-boundary of the computational domain). Because the spacing of the dominant critical layers in x is Δ, the spacing of the initial vortices is also Δ and the vortices run into each other and stop growing when ${L}_{{zv}}\simeq {\rm{\Delta }}$. (We explicitly chose the parameters in our studies to have ${L}_{x}\gt {\rm{\Delta }}$ so that L_zv would be determined by Δ and not by L_x.) However in a fully compressible flow, it is well known that vortex diameters are often constrained by the fact that the vortices must be subsonic, otherwise the vortices rapidly dissipate their energy via acoustic radiation. In particular, the change in the velocity from one side of the vortex to the opposite side must be less than the sound speed. Therefore, in fully compressible flows, vortices stop growing when they become large enough to be supersonic, which can occur while the vortices are still too small to run into each other. By definition of the Rossby number, the characteristic relative velocity of a zombie vortex is $| {{Ro}}_{\mathrm{rms}}\;{{\rm{\Omega }}}_{0}\;{L}_{{zv}}|$. However, the characteristic change in the velocity of the total flow across the x-diameter of the vortex is ∼ $| {{\rm{\Omega }}}_{0}{L}_{{zv}}|$ because the vortex is embedded in the Keplerian disk with a shear of $| q{{\rm{\Omega }}}_{0}|$. For the velocity change across the diameter of the vortex to be less the sound speed, L_zv must be less than ∼H. In a PPD, $H\ll {\rm{\Delta }}$, so we would expect that the diameters of zombie vortices have their growth limited by the constraint that the vortices remain subsonic, and not by the size Δ. Using ${L}_{{zv}}\simeq H$ in Equation (3.7) makes the Mach numbers of the zombie vortices in a PPD of order $| {{Ro}}_{\mathrm{rms}}|$. Note that in some fully compressible flows that it is possible that H is greater than Δ. In fact, in Section 3.2 we showed that for the fully compressible flow examined in this paper that ${H}_{0}\simeq 5{\rm{\Delta }}$. Therefore for the fully compressible flow examined in this paper, we expect, and have observed numerically, that ${L}_{{zv}}\simeq {\rm{\Delta }}$, that vortex diameters are not restricted by the subsonic constraint, that the zombie vortices remain subsonic, and that Equation (3.6) is valid.

Figures 4(d), 5(d), and 6 show that at late times zombie turbulence fills the computational domain. We can quantify this space-filling property by defining a spatial filling factor for the turbulent vorticity: ${f}_{\mathrm{Ro}}(\eta ,t)$ is the volume fraction of the computational domain that has $| {Ro}(x,y,z,t)| \geqslant \eta$. Figure 11 shows that for the anelastic flow in Figure 8, approximately 10% of the flow is filled with vortices with Rossby numbers with magnitudes greater than 0.3; 30% with magnitudes greater than 0.2; and almost 60% with magnitudes greater than 0.1. These fill factors may still be growing in time. The fill factors in Figure 11 are representative of all of our anelastic calculations. For the flow in our fully compressible Athena simulations, the fill factors are slightly larger with approximately 20% of the flow is filled with vortices with Rossby numbers with magnitudes greater than 0.3; 45% with magnitudes greater than 0.2; and almost 75% with magnitudes greater than 0.1.

Zoom In Zoom Out Reset image size

We have shown that Keplerian disks are unstable to the purely hydrodynamic ZVI instability when a vertical component of gravity from the central star or protostar is included and the vertical thermodynamic structure is stably stratified. It may seem counterintuitive that introducing stable stratification makes a flow more unstable. Yet we have demonstrated that vertical stratification introduces baroclinic critical layers into the flow. In this work, we have explored ZVI in the limit of infinite cooling time. ZVI has a fast growth rate, requiring only a few hundred local "years" to produce large-magnitude turbulence. Even if the instability is initially spatially confined to a small region, the entire domain fills rapidly with turbulence with rms Rossby numbers between 0.2 and 0.3. We have argued that in a fully compressible PPD, we expect that the Mach number of the space-filling zombie turbulence will be approximately equal to its Rossby number. When turbulence created by ZVI reaches a statistically steady state, it is neither homogeneous nor isotropic, and the pattern of the turbulence is highly asymmetric with respect to cyclonic and anticyclonic vorticity. The flow is dominated by large, elongated, anticyclonic vortices and narrow cyclonic layers of vorticity that are aligned in the stream-wise direction (i.e., the azimuthal direction in a PPD). This pattern makes zombie turbulence easily identifiable and distinct from other forms of turbulence. In a PPD, the vortex size and the spacing between the layers would be expected to be of order the vertical pressure scale height. Our calculations show that zombie turbulence is robust and does not decay. In MPJH13, the permanence of the turbulence was shown to be due to the fact that it draws its energy from the energy stored in the Keplerian shear. Disks filled with statistically steady zombie turbulence appear to be attracting solutions that are independent of the details of the initial conditions that trigger the turbulence; properties, such as the energy of the non-Keplerian velocity, the energy spectrum as a function of wavenumber, and the probability distributions of the Mach and Rossby numbers are independent of the triggering perturbations.

The robustness of ZVI prompts the question of how it was missed in previous studies. Part of the answer is that ZVI has four necessary ingredients: rotation, a component of gravity along the (vertical) rotation axis, horizontal shear, and stable vertical stratification. These ingredients are all present in a PPD, but many previous studies neglected stable vertical stratification. Any stability study that does not include all four ingredients—whether a numerical calculation or a laboratory experiment that uses a constant-density fluid—cannot produce ZVI.

However, there are several others reasons that previous studies did not find ZVI.

• 1.  
  Sufficient numerical resolution and/or control over numerical dissipation is required to resolve critical layers. Figures 4–6 show that these layers are very thin. Simulations using the Godunov finite-volume code Athena show that ZVI is present only when more than 128 grid points per vertical pressure scale height are used in the radial direction. Thus the instability could not be observed in the simulations of vertically stratified disks by Fromang & Papaloizou (2006) who used 30 radial points per pressure scale height or by Fleming & Stone (2003), Oishi & Mac Low (2009) who used 64 or fewer.
• 2.  
  Sufficient integration time is needed for vorticity at small length scales to reverse cascade to larger scales, for the cyclonic vorticity to homogenize, and for the anticyclones to merge into vortices that are sufficiently large and spatially scattered (as in Figure 1(b)) to excite critical layers. It then takes additional time for dipolar vortex layers to form and several more vortex turn-around times for them to roll up and create the first generation of zombie vortices. In our studies (cf., Figure 3), it takes at least 50 years for the zombie vortices to fill the domain, for energy to be drawn from the Keplerian shear by the vortices, and for the non-Keplerian energy component to increase substantially. If a calculation were terminated early (e.g., after 8 years, as in Figure 2), then ZVI would not be observed.

In BM05, we hypothesized that vortices themselves could directly transport angular momentum radially outward in a PPD. However, we found that the rate of transport, as parameterized by $\alpha \equiv \langle \rho {v}_{x}\;{v}_{y}\rangle /({c}_{{\rm{s}}}^{2}\langle \rho \rangle )$, was a factor of ∼100 times smaller than the values needed for the observed rate of star formation from a protostar, where $\langle Q\rangle$ means the volume-averaged value of Q over the entire domain. The value of α is of order:

$$\begin{eqnarray}&&\alpha \sim {\mathrm{Ma}}_{\mathrm{rms}}^{2}\;{f}_{\mathrm{Ma}}\;{f}_{\rho }\;c,\end{eqnarray} \tag{ 4.1 }$$

where c is the correlation between v_x and v_y, ${f}_{\mathrm{Ma}}$ is the time-averaged spatial fill factor of the flow with ${\mathrm{Ma}}_{\mathrm{rms}}$, and ${f}_{\rho }$ is the ratio of the Mach-number-weighted average density of the gas to the density of the gas in the mid-plane of the PPD. In BM05, the small values of α were due to the small values of the correlation c because of the inherent symmetry of a vortex. Vortices in the disk are approximately mirror symmetric with respect to the x–z plane that passes through the vortex center; thus, the values of $\langle \rho {v}_{x}\;{v}_{y}\rangle$ on either side of the symmetry plane nearly cancel. We therefore believe that zombie vortices cannot directly produce a sufficiently large α for star formation.

We now hypothesize that zombie turbulence in a fully compressible fluid indirectly may produce much larger values of α. A fully compressible fluid has three orthogonal velocity components (Chandrasekhar 1961): two are rotational and divergence-free (i.e., the poloidal and toroidal components, which anelastic flows have), and one is an irrotational component with a non-zero divergence (i.e., the curl-free potential flow, which anelastic flows do not have). We suspect that the rotational divergence-free zombie turbulence will excite acoustic modes (i.e., equipartition of energy among the poloidal, toroidal and potential components). Johnson & Gammie (2005), Shen et al. (2006), Lesur & Papaloizou (2010), Lyra & Klahr (2011), Raettig et al. (2013) have shown that acoustic waves in a PPD transport angular momentum radially outward with values of α of order ${\mathrm{Ma}}_{\mathrm{rms}}^{2}$. Acoustic waves have correlations c of order unity. For acoustic waves launched by zombie vortices, we estimate $c\simeq 1$, ${\mathrm{Ma}}_{\mathrm{rms}}\simeq 0.2$, ${f}_{\mathrm{Ma}}\simeq 0.2$, and ${f}_{\rho }\simeq 0.3$, which yields a value of $\alpha \sim 2\times {10}^{-3}$, consistent with values inferred from observed rates of star formation. Our estimate that the fill factor ${f}_{\rho }$ is small is based on the numerical simulations by BM05 and Lesur & Papaloizou (2009) that found that vortices can exist above and below the disk midplane, but the midplane itself at $| z| \lt 0.1\;H$ is devoid of laminar vortices with $| {Ro}| \gt 0.1$. The calculation of α in a fully compressible calculation of a PPD will be the subject of future papers. However, in this paper we have shown, regardless of the value of α, that ZVI destabilizes a PPD and creates space-filling turbulence with rms Rossby numbers between 0.2 and 0.3.

Other outstanding issues include: (i) varying different profiles of stable vertical stratification, (ii) determining how close to the nearly unstratified midplane can turbulence from ZVI can penetrate, (iii) including explicit cooling terms, (iv) studying ZVI in global simulations, (v) exploring how ZVI interacts and competes with MRI, (vi) investigating how ZVI interacts and competes with streaming instabilities, and (vii) computing dust capture inside zombie vortices.