The Dead Zones of Protoplanetary Disks are Not Dead

We show that the “dead” zone of a protoplanetary disk fills with robust 3D vortices from a purely hydrodynamic instability. This new instability is not linear and requires a weak finite-amplitude initial perturbation. The instability was not seen previously either due to a lack of numerical spatial resolution, or because many previous simulations either ignored vertical gravity or had initial flows with constant density. Our new finite-amplitude instability is due to a family of previously-unknown critical layers that form in rotating, shearing, vertically stratified flows like those in protoplanetary disks. Initial perturbations of white noise (with Mach numbers much less than unity), waves, or vortices can trigger the instability. A small-volume, small-amplitude initial vortex confined to one part of the disk can fill the disk with vortices by exciting a nearby critical layer. The critical layer produces an intense vortex layer that rolls-up to form vortices with large-amplitudes and volumes. This 1st generation of vortices then sheds waves that excite nearby critical layers, which in turn, create a 2nd generation of vortices with large amplitudes and volumes. The mechanism of exciting nearby critical layers and turning them into large vortices self-similarly, self-replicates until large vortices fill the disk at all radii.


Can Non-Magnetically-Coupled, Non-Self-Gravitating Keplerian Disks Have
Purely Hydrodynamic Instabilities? Can They Drive Angular Momentum Transport?
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 the task because the timescale is several orders of magnitude longer than the age of the universe: τ visc ∼ R 2 /ν ∼ 100, 000 Gyr for conditions typical at R ∼ 1 AU in a protoplanetary disk (PPD) with kinematic viscosity ν ∼ v th , thermal velocity v th ∼ 1 km/s, molecular mean free path ∼ 1 cm. Shakura & Sunyaev (1973) proposed that Keplerian differential rotation may be susceptible to shear instabilities that would drive turbulent diffusion, for which they introduced an "eddy" viscosity: ν eddy = αc s H (the idea being that turbulent "blobs" would move at subsonic speeds over distances of order the scale height, and the parameter α would subsume all our ignorance over the details of the transport efficiency). An oft-claimed counterpoint to the shear instability hypothesis is that Keplerian differential rotation is stable because its specific angular momentum increases with radius. 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 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 Tevzadze et al. 2003;Umurhan & Shaviv 2005;Afshordi et al. 2005;Mukhopadhyay et al. 2005;Tevzadze et al. 2008;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). Velikhov (1959) and Chandrasekhar (1960Chandrasekhar ( , 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 . However, relatively cool and nearly neutral PPDs 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. 2009(Matsumura et al. , 2007. 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. 2000Li et al. , 2001Varnière & Tagger 2006;Richard et al. 2013;Meheut et al. 2010Meheut et al. , 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,b;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 not necessarily baroclinic. Still others have investigated vertical shear instabilities related to the Goldreich-Shubert-Fricke instability in rotating stars (Nelson et al. 2013;Goldreich & Schubert 1967;Fricke 1968).
Of course, by direct observation, we know that some mechanism exists 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." The exact nature of the turbulent state may well be needed to understand future observations of PPDs, say with the James Webb Space Telescope. Moreover, if the main claim of this paper is correct, and laminar PPDs do not exist because they are always filled with large-amplitude turbulence, then the details of that turbulence may be needed to understand the late stages of star formation and the observed star formation rate, and the "universal" (if there is one -see Bastian et al. (2010)) stellar initial mass function (IMF). Moreover, turbulent, rather than laminar, PPDs could radically change the mechanisms that sweep and accumulate dust in a disk and therefore affect our modeling of the observation of the infrared emission from dust in a PPD, as well as an IMF for planets.

A New Instability That Leads to Space-Filling Robust "Zombie" Vortices & Turbulence
In Barranco & Marcus (2005, hereafter BM05), we investigated the stability of three-dimensional vortices in the midplanes of protoplanetary disks. 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 protoplanetary disk (e.g., spatially varying gravity and Brunt-Väisälä frequency) and numerically investigated simple Couette flow with constant gravity and constand 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 sheets 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 an "inverse 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 sheets (Marcus 1993). After a period of inverse 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 sheets (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.
The crucial new insight in MPJH13 was the recognition of critical layers as being sites that are receptive to perturbations. Formally, critical layers are special locations in the flow where the highest derivative in the linearized equations have coefficients that vanish, indicating that the eigenfunction solutions are singular at those locations (Maslowe 1986;Drazin & Reid 1981). We defer the detailed mathematics until Paper II, but note that the first such critical layers occur at distances (along the radial/cross-stream direction) ∆ ≡ N L y /(2πσ) on either side of the vortex, where N is the Brunt-Väisälä frequency, L y is the azimuthal size of the domain, and σ is the local shear rate. 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. The eigenmodes are neutrally stable (do not grow nor decay) to small perturbations, but the isolated anticyclone can act as a forcing and excite the critical layers. A dipolar vortex layer develops at the location of the 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 a new anticyclone will eventually be spawned at the location of the critical layer (panel e). This new anticyclone will then excite its neighboring critical layer (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 hyperviscosity and 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.

Paper Outline
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 §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 §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 sheets.
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 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, and would require at least 3 orders of magnitude more computational effort to get to the same effective resolution as the pseudo-spectral code.
In our discussion in §4, we explain why ZVI was not observed in previous studies by others, we review ZVI in the context of other hydrodynamic and magnetohydrodynamic instabilities in PPDs, and we consider the role of ZVI in angular momentum transport and planetesimal formation. In 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 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. (a) Domain is initially filled with Kolmogorov turbulence. In time, the background shear will stretch cyclonic regions into vortex sheets, while anticyclonic regions will roll-up into compact vortices. Smaller vortices will merge into larger vortices, a feature of "inverse cascades" in quasi-2D turbulence. (b) After a period of vortex mergers, a small region will have a few somewhat isolated anticyclones. Here, we show one dominant anticyclone. The dotted line indicates the location of the critical layer at a distance ∆ from the perturbing anticyclone. (c) Anticyclone acts as a source of forcing, exciting neighboring critical layers and creating a dipolar vortex sheet. 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. (d) Vortex sheets with same orientation as background shear are linearly unstable. (e) Anticyclonic vortex sheet rolls-up into a compact anticyclone. (f) New vortex excites neighboring critical layer and the process repeats. an appendix, we give a "recipe" for observing ZVI in Athena. In our companion Paper II, we shall confront ZVI on a more mathematical and theoretical level and address thresholds required to trigger ZVI. In particular, we shall show why extraordinarily weak Kolmogorov-like noise is so efficient in exciting ZVI.

HYDRODYNAMIC EQUATIONS IN A CARTESIAN DOMAIN
Consider a PPD in which the steady flow is only in the azimuthal direction with angular velocity where R is the cylindrical radius from the protostar. A value of q = 1.5 corresponds to Keplerian rotation. Now consider a three-dimensional box of gas centered at radius R 0 from the protostar that co-rotates with the gas with angular rate Ω 0 ≡ Ω(R 0 ). The box is sufficiently small that we can ignore curvature of the disk and map the cylindrical coordinates onto local Cartesian coordinates: (Hill 1878;Goldreich & Lynden-Bell 1965). In what follows, symbols with "hats" will denote unit vectors:x points in the local outward radial direction,ŷ points in the local azimuthal direction, andẑ points in the local vertical direction.

Compressible Euler Equations
Euler's equations (continuity, momentum and energy with no viscous, radiative or diffusive terms) for the gas in the rotating frame are: where 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Ω 0ẑ × v is the Coriolis force, and the term 2qΩ 2 0 xx is the tidal acceleration that arises from the difference between the inward central force that causes circular motion RΩ 2 (R) and the outward centrifugal force RΩ 2 0 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: where T is the gas temperature, R is the gas constant, and γ ≡ C P /C 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 γ = 5/3.
A time-independent equilibrium (denoted with overbars) to the above set of equations is: Because there is no radiative or diffusive effects in equations (2-2), there is a degeneracy of allowable equilibrium thermodynamic profiles, so the equilibrium temperatureT (z) corresponding to the equilibrium velocityv is not unique. In fact, one thermodynamic quantity,P (z),T (z),ρ(z), or¯ (z) can be arbitrarily specified, and the others will follow. In this work, we shall assumē where ρ 0 is the midplane density and z = 0 is the midplane of the disk.

Anelastic Approximation
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,b). We have previously used the anelastic approximation to study three-dimensional vortices in PPDs (Barranco et al. 2000;Barranco & Marcus 2005 and the Kelvin-Helmholtz instability of settled dust layers in PPDs (Barranco 2009;Lee et al. 2010a,b). 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: 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 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: where N (z) is the Brunt-Väisälä frequency, a measure of the convective stability of the atmosphere. In general, N 2 ∝ g · ∇s, wheres is the background entropy. Here, we focus on vertical entropy gradients, so: 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,P (z),T (z),ρ(z), or¯ (z) can be arbitrarily specified. As before, we arbitrarily set the background equilibrium temperature to be constant, 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.

Background Equilibrium States
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): 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 0 :

Boundary Conditions and Numerical Resolution
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 → ±∞ (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ρ(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. In the Godunov finite-volume calculations with Athena, the vertical velocity was set to zero at the top and bottom boundaries -see Appendix A for details.
The Godunov finite-volume computations with Athena use 256 grid points in each spatial dimension; whereas, the anelastic pseudo-spectral computations use 256 pseudo-spectral modes in each direction. Note that the Athena code has no explicit dissipation, but its Godunov finite-volume algorithms are inherently dissipative. The algorithms used in the pseudo-spectral calculations have no inherent dissipations, so an explicit, weak hyperviscosity is required.

Choices of Parameters and Dimensionless Numbers
There are three (reciprocal) timescales of interest: the rotation rate Ω 0 , the shear rate σ ≡ −qΩ 0 , and the Brunt-Väisälä frequency N 0 . We consider only Keplerian shear, q = 1.5 (with the exception of the calculations reported in Fig. 2). We define the dimensionless ratio β ≡ N 0 /Ω 0 . For a PPD with linear vertical gravity and no vertical temperature gradient, β(z) = 0.63|z|/H 0 . In BM05, we observed that vortices naturally formed in stratified regions above and below the midplane: 1 |z|/H 0 4, corresponding to 0.63 β 2.5. For constant gravity g 0 , one can show from Eq. (2-10) that: With constant gravity, the only independent length scale in either the fully compressible or anelastic equations is the gas pressure scale height H 0 . 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. Appendix B shows that for the anelastic and fully compressible flows presented here with constant vertical gravity: Readers who may be troubled by the fact that ∆ depends on L y , a seemingly arbitrary choice of our computational domain, rather than a physical length such as H 0 , should read Appendix B. The length scale ∆ is important because zombie vortices form at the excited critical layers, and they, in turn, perturb and excite new critical layers. In MPJH13 we showed that when an equilibrium flow was perturbed with a single vortex, the domain filled with zombie vortices that were evenly spaced in x with a separation of ∆. In §3 we show that when the equilibrium flow is perturbed with random noise, rather than a single vortex, the domain also fills with zombie vortices that are nearly evenly spaced in x and have a separation approximately equal to ∆. Furthermore, the volumes of the zombie vortices grow until they nearly run into each other. Thus, the characteristic x-diameter of mature zombie vortices is also ∼∆. In § 3.3, we explain why, for a fully compressible PPD with linear gravity, we expect the radial separation and radial size of zombie vortices will be limited by ∼H, rather than ∆.

EVIDENCE OF INSTABILITY IN PROTOPLANETARY DISKS WITH VERTICAL GRAVITY
Unless otherwise stated, in this section, all times are reported in units of orbital periods τ orb ≡ 2π/Ω 0 or local "years"; all lengths in units of ∆; all velocities in units of Ω 0 H 0 , and all energies per unit mass in terms of (Ω 0 H 0 ) 2 . Table 1

Temporal Growth and Decay of Initial Energy Fluctuations
One of the most cited pieces of evidence that PPDs are stable to purely-hydrodynamic instabilities is given by Fig. 1 in Balbus et al. (1996), which we copy below as Fig. 2a. 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 Eq. (2-1) and where the kinetic energy fluctuation is defined as |v(x, y, z, t) −v yŷ | 2 /2, wherev yŷ is the steady equilibrium flow in Eq. (2-4a). The initial-value calculations that produced Fig. 2a used the fully compressible equations (2-2)-(2-3), and were initialized with the steady equilibrium velocity in Eq. (2-4a). However, this flow was computed with no vertical gravity and was therefore initialized with the equilibrium that hadP = P 0 ,ρ = ρ 0 , andT = T 0 . The initial flow was perturbed with small-amplitude noise (with root-mean-square or 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 ≥ 2 and decrease in time for q < 2, Fig. 2a 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 Ω(R) ∝ 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 Fig. 2a corresponds to the case of pure Cartesian shear with the Coriolis and tidal acceleration terms dropped from Eq. (2-2b), which would be appropriate for fully compressible flow in a channel with cross-stream shear, but no rotation.
In Fig. 2b, we reproduce the results of Fig. 2a 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) ∝ 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 exact stability results for no vertical gravity with a pseudo-spectral code (vs. finite-volume code) which solve the anelastic equations (vs. the fully compressible equations).
Time Fig. 2.-(a) Left panel: Fig. 1 from Balbus et al. (1996), which shows the temporal evolution of the fluctuation kinetic energy per unit mass, for different values of q as defined by Eq. (2-1). These are fully compressible simulations with g = 0, N = 0, γ = 5/3. The calculation was done with the finite-difference Zeus code with spatial grid of 63 3 (127 3 for q = 1.5, 2) points. The equilibrium flow had uniform pressure, density, and temperature. The initial rms velocity amplitude of the noise was ∼ 0.1. (b) Right panel: Same as left panel, but computed with our pseudo-spectral code and the anelastic equations, rather than the fully compressible equations. The initial noise had an rms velocity of ∼ 0.02.
How does vertical stratification affect stability? The decaying dashed curve in Fig. 3a is the time evolution of the kinetic energy fluctuations for a unstratified flow similar to the one as given by the curve in Fig. 2b with q = 1.5 and g = N = β = 0, but plotted for a longer time and a smaller value of initial noise. The solid curve in Fig. 3a shows the evolution of kinetic energy fluctuations when vertical gravity is included; here, g = g 0 = 0 is constant, N = N 0 = 0, and β = 2. The growth of the kinetic energy fluctuations shows that the flow is unstable with q = 1.5 when vertical gravity is present. Fig. 3b shows the kinetic energy growth for the same flow as shown with the solid curve as in Fig. 3a, but computed with the fully compressible equations with the Godunov finite-volume code Athena. The rapid growth of the kinetic energy fluctuations shows that ZVI is not an artifact of the anelastic approximation or spectral codes.   The Keplerian background shear has stretched the relative vorticity and elongated it in the y direction; the flow's energy has cascaded to small scales where it is dissipated by hyperviscosity causing the amplitude of |Ro(x, y, z, t)| to decrease by an order of magnitude from its initial value. (c) t = 50.9 yrs: Mergers of small anticyclones in an inverse cascade of energy to larger scales has led to a few surviving, but spatially scattered, anticyclones. (d) t = 1370 yrs: Zombie turbulence and zombie vortices with Ro −0.3 fill the domain. The near spatial periodicity of the flow in x, with a wavenumber (in this case ∼7) slightly smaller than L x /∆ is one of the signatures that makes zombie turbulence easy to identify.   Fig. 5d. It is possible that at this time, the flow has not reached a statistically steady state, and the pixelated structures in panel (b) may indicate that the flow is slightly under-resolved spatially. However, the point we emphasize here is that zombie turbulence occurs in fully compressible flows computed with a finite-volume code, and the turbulence looks qualitatively similar to that in anelastic flows computed with a spectral code. See Appendix A for a precise recipe to produce ZVI with the Athena code.

The Vertical Vorticity of the Zombie Vortex Instability
The growth of kinetic energy fluctuations in Figs. 3 is evidence of instability, but spatial plots of the relative vorticity, defined as ω ≡ ∇×(v−v) = ∇×v+q Ω 0ẑ , 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 Ro(x, y, z, t) ≡ [ẑ · ω(x, y, z, t)]/(2Ω 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 §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. Fig. 4 shows 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.) Fig. 5 shows 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 Figs. 4 and 5 are similar to the flow shown by the solid curve in Fig. 3a with the difference that β = 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 (see Paper II). Thus, Figs. 4a and 5a 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 Figs. 4b and 5b. The stretching cascades the initial kinetic energy to spatial wavenumbers with higher Fourier modes where it is quickly dissipated by hyperviscosity (similar to Fig. 3a at t 50). In competition with the stretching, there is a weaker inverse cascade of energy from small to large length scales due to the fact the fluid is rapidly rotating. Here, the inverse 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 sheets (Marcus 1993). The shear in a PPD is anticyclonic, and therefore anticyclones merge, while cyclones are stretched. 2 The result of the competition between the forward energy cascade (resulting in energy dissipation) and the inverse energy cascade (resulting in the mergers of anticyclones) is visible in Figs. 4c and 5c at time t = 51 yrs; 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 Fig. 1.
By 51 yrs, ZVI is well underway in Figs. 4 and 5. The scattered anticyclones which emerged from the inverse cascade have now excited neighboring critical layers to form dipolar vortex sheets. 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 Ro(x, y, z, t). However, the vortex sheets that the critical layers create are easily identifiable in plots of Ro(x, y, z, t) because each dipolar sheet appears as a pair of oppositesigned "stripes" in the x-y plane (with the critical layer sandwiched between them - Fig. 1c). In an initial-value simulation with an initial perturbation consisting of noise, it is sometimes difficult to see the dipolar vortex sheets due to the fast instability of the anticyclonic sheet. However in carefully controlled simulations consisting of a single perturbing vortex, rather than noise, dipolar sheets can be found easily, as shown in Fig. 1a in MPJH13. As discussed in §1.2 (and elaborated on in Paper II), the sheets of cyclonic vorticity aligned in the stream-wise direction are linearly stable. In contrast, sheets of anticyclonic vorticity aligned in the stream-wise direction are linearly unstable and roll-up into stable anticyclones.
At late times (Figs. 4d and 5d), the flow has reached a statistically steady state of zombie turbulence. One can see by eye, especially in Fig. 5d, 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 2-dimensional sheets 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 /∆ = 3π 10, and at times between that of Figs. 5c and that of Fig. 5d, 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 /∆. 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 Ro(x, y, z, t) always look like the pattern in Figs. 4d and 5d. Our one fully compressible calculation (Fig. 6) also shows this same pattern at late times. In Fig. 6, β = 2 and L y = L x , so we would expect that L x /∆ 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 Fig. 6 has been reduced to ∼4 by vortex mergers and inverse energy cascades. Note that the two flows shown in Figs. 8 and 9, also have β = 2 and L y = L x , and although we have not included figures of their 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 sheets and elongated anticyclones aligned in the stream-wise direction of the background shear, and the pattern of sheets and anticyclones has a dominant x wavenumber slightly less than L x /∆.
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 (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 Eq. (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, Eq.
(3-1) shows that in a Keplerian disk, the aspect ratios χ of zombie vortices are between 4 and 5.

Space-Filling Properties of Zombie Turbulence
Figures 7 -10 show how the root-mean-square (rms) Rossby numbers Ro rms (t) and rms Mach number Ma rms (t) evolve in time for three anelastic flows and one fully compressible flow. The flows in Figs. 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 Ma rms (t), of the late-time flows in Fig. 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 Figs. 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.
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 rms (t)| ∼ 0.2 − 0.3 indicates large amplitude turbulence. We have no 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 rms | is consistent with the picture that we presented in MPJH13 and Fig. 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. In Paper II, we show that an initial perturbation must have |Ro rms | 0.2 in order to initiate ZVI in an unperturbed flow. Therefore, our numerical finding that zombie vortices all have strengths greater than 0.2 is consistent with the requirements of our picture of how ZVI spreads.
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 Ma 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 Ma rms (t) is only ∼ 10 −2 in our calculations. We now show that the dependence of Ma rms on the size of the computational box is easily explained, and that in a PPD we expect Ma rms to be equal to |Ro rms | at late times. Note that the rms Rossby number is approximately where L ω 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. 3 The rms Mach number is approximately 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 L ω because the velocity is dominated by the large length scales and the vorticity by the small length scales -see Paper II. 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 ω are both approximately the mean size of the vortices. Therefore, in a flow dominated by zombie vortices, 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 γ = 5/3), Eq. (2-11) determines c s , and Eq. (3-4) then shows that a flow with zombie vortices has Ma rms 2/5 Ro rms β −1 (L zv /H 0 ). (3)(4)(5) In all of the calculations presented in this paper, L zv ∆, so using Eq. where we used the fact that in our computations L y = H 0 . Equation (3-6) shows that Ma rms depends on the size of the computational box, is consistent with Figs. 7 -9 at late times, and explains why our calculations produce zombie turbulence with Ma rms 10 2 . 4 3 In the literature of turbulence, an "eddy" is not the same as a Fourier mode, rather it is a component of the velocity made from a band of Fourier modes. The kinetic energy per unit mass of an eddy is defined in terms of the differential kinetic energy spectrum E(k) of the flow, where k is the Fourier wavenumber. The kinetic energy per unit mass of the total flow is ∞ 0 E(k) dk, and the kinetic energy per unit mass of an eddy of diameter L is defined as The last equivalence is what is used to defined the eddy velocity V eddy (L). Note that an eddy of diameter L is made up of all velocity Fourier modes with wavenumbers 2π/(2L) < k < 2π/L. Kolmogorov turbulence has an energy spectrum E(k) ∝ k −5/3 , so in Kolmogorov turbulence, V eddy (L) is proportional to L 1/3see Paper II or Tennekes & Lumley (1972).
Equation (3-6) does 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 Eq. (2-9b) rather than Eq. (2-11). Therefore, for a zombie-vortex-dominated flow in a PPD, Eq. (3-4) becomes Ma rms (L zv /H) Ro rms . (3-7) The second reason that Eq. (3-6) does not apply to a PPD is we expect L zv H for zombie vortices in a PPD, rather than L zv ∆. 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 the vortices run into each other and stop growing when L zv ∆. (We explicitly chose the parameters in our studies to have L x > ∆ 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 rms Ω 0 L zv |. However, the characteristic change in the velocity of the total flow across the x-diameter of the vortex is ∼ |Ω 0 L zv | because the vortex is embedded in the Keplerian disk with a shear of |qΩ 0 |. For the velocity change across the diameter of the vortex to be less the the sound speed, L zv must be less than ∼ H. In a PPD, H ∆, 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 ∆. 5 Using L zv H in Eq. (3-7) makes the Mach numbers of the zombie vortices in a PPD of order |Ro rms |. Note that in some fully compressible flows that it is possible that H is greater than ∆. In fact, in §3.2 we showed that for the fully compressible flow examined in this paper that H 0 5∆. Therefore for the fully compressible flow examined in this paper, we expect, and have observed numerically, that L zv ∆, that vortex diameters are not restricted by the subsonic constraint, that the zombie vortices remain subsonic, and that Eq. (3-6) is valid.
Figs. 4d, 5d, 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 Ro (η, t) is the volume fraction of the computational domain that has |Ro(x, y, z, t)| ≥ η. Fig. 11 shows that for the anelastic flow in Fig. 8, approximately 10% of the flow is filled with vortices with Rossby numbers with magnitudes greater than 0.3; 30% with magnitudes greater 5 As discussed in Appendix B, for flows in which the Brunt-Väisälä frequency is not constant, but rather is a function of z, the distance ∆ between critical layers is still given by Eq. (2-12), but N0 must be replaced with N (z). The critical layer spacing is then ∆ = (1/3π) 2/5(Ly/H)|z| for a PPD. Setting Ly = 2πR, we obtain ∆ 0.42(R/H)(|z|/H). In general, R H, whereas for regions of interest, |z| ∼ H, so H ∆. 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 Fig. 11 are representative of all of our anelastic calculations. For the flow in Figs. 4 and 5, 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.  Fig. 7, but for the anelastic flow corresponding to the solid curve in Fig 3a. This flow has the same initial amplitude of the Kolmogorov noise as the flow in Fig. 7, but has β = 2 rather than 1.

M a rms (t)
Time t Time t Ro rms (t) Fig. 9.-Time evolution of Ma rms and Ro rms as in Fig. 8. The flow parameters are the same as in Fig. 8, but here the rms Mach number of the initial Kolmogorov noise is two-thirds the value in Fig. 8. After t 500 yrs, many of the statistical properties of the flows in Fig. 8 and Fig. 9 are nearly the same, which shows that the flows are being drawn to an attractor that is independent of the details of the initial conditions.

Summary
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 in the equations of motion. 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 sheets 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 sheets 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.

Why ZVI Was Not Seen in Previous Studies
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 vertical density stratification. These ingredients are all present in a PPD, but many previous studies neglected vertical density 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.
• Sufficient numerical resolution and/or control over numerical dissipation is required to resolve critical layers. Figs. 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.
• Sufficient integration time is needed for vorticity at small length scales to inverse 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 Fig 1b) to excite critical layers. It then takes additional time for dipolar vortex sheets 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., Fig. 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 Fig. 2), then ZVI would not be observed.

ZVI in Protoplanetary Disks
ZVI was not seen in previous disk studies that were initiated with flows that were far from equilibrium. These types of initial-value experiments cannot be used to determine whether an equilibrium is stable but are useful to investigate whether there are multiple statistically-steady states of the flow. For example, an initial disk flow and magnetic field that are not in equilibrium might evolve to a steady disk with Keplerian velocity, indicating that laminar Keplerian flow is a strong, possibly unique, attractor. On the other hand, the flow could evolve to a statisticallysteady, non-Keplerian flow with sustained turbulence, indicating that there are multiple solutions to the equations of motion and that a state of sustained turbulence is a robust attractor. Constructing numerical experiments to test whether the magnetorotational instability (MRI) overtakes a competing ZVI or vice versa, or whether the ZVI overtakes a competing streaming instability is more complicated than comparing the linear growth rates of eigenmodes and requires a careful consideration of the initial perturbations, not only because the ZVI requires a finite (albeit small -see Paper II) amplitude initial perturbation, but also because the winner of such a competition also depends upon the full nonlinear dynamics. There are situations where it is important to know which instability dominates. For example, if ZVI rapidly creates a turbulent PPD whose turbulence prevents the streaming instability, then it would be impossible to create planetesimals via clumping from the streaming instability. However, the concept that there is a competition among the various instabilities of an equilibrium PPD may possibly be misleading. Numerical simulations by Krumholz et al. (2007) show that when a pre-stellar core collapses and a PPD begins to form, the disk is highly turbulent. Currently, there is no evidence that the late-time disk formed from collapse does not remain in a state of sustained turbulence or that the flow ever evolves to a steady, completely laminar Keplerian disk. It is also possible that while a disk is still forming from a collapsing core that the large spatial fluctuations in its velocity, temperature, and density evolve into features that resemble an MRI instability, a zombie vortex, or the signature of some other instability of a steady Keplerian disk. There is no evidence that the instability that overtakes its competition in a forming turbulent disk is the same instability that overtakes its competition in a steady disk.

α and Angular Momentum Transport
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 α ≡ ρv x v y /(c 2 s ρ ), was a factor of 10-100 times smaller than the values needed for the observed rate of star formation from a protostar, where Q means the volume-averaged value of Q over the entire domain. The value of α is of order: where c is the correlation between v x and v y , f Ma is the time-averaged spatial fill factor of the flow with Ma rms , and f ρ 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 ρv x v y 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 produces 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 also 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). In a turbulent fully compressible gas, there is evidence that there is an equi-partition of energy among the three orthogonal velocity components. In particular, Kritsuk et al. (2007Kritsuk et al. ( , 2009 and Lemaster & Stone (2009) showed numerically that when a fully compressible fluid in a disk is forced by sustained rotational turbulence, the potential flow becomes turbulent, and the energy of the potential flow comes into equi-partition with the two rotational components of the turbulence. Therefore, sustained zombie turbulence in a fully compressible flow could, by energy equi-partition, create irrotational turbulence with acoustic waves with Mach numbers of ∼ 0.2 -0.3. (N.B. We could not test these ideas in the simulation of the zombie turbulence in the fully compressible flow in Fig. 6 Raettig et al. (2013) have shown that acoustic waves in a PPD transport angular momentum radially outward with values of α of order Ma 2 rms . Acoustic waves have correlations c of order unity. For acoustic waves launched by zombie vortices, we estimate c 1, Ma rms 0.2, f Ma 0.2, and f ρ 0.3, which yields a value of α ∼ 2 × 10 −3 , consistent with values inferred from observed rates of star formation. Our estimate that the fill factor f ρ is small is based on the numerical simulations by BM05 and Lesur & Papaloizou (2009) that found that vortices can exist above and below the the disk midplane, but the midplane itself at |z| < 0.1H is devoid of laminar vortices with |Ro| > 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. Laminar PPDs do not exist. Whatever mechanisms that radially transport angular momentum and that accumulate dust will be different in a turbulent disk than they would be in a steady laminar disk. This suggests that numerical and theoretical investigations of star and planet formation in laminar disks might be unrelated to the real dynamics of astrophysical disks.
PSM is supported in part by NSF grants AST-0905801 and AST-1009907, and by NASA PATM grants NNX10AB93G and NNX13AG56G. Part of the computational work used an allocation of computer resources from the Extreme Science and Engineering Discovery Environment (XSEDE), which was supported by National Science

A. APPENDIX A -How to Create ZVI with the Athena Code
Athena is a Godunov finite-volume code that can solve the fully compressible hydrodynamic equations in a shearing box Stone et al. 2008;Stone & Gardiner 2010). We have confirmed that ZVI can be simulated with Athena, so ZVI is not an artifact of the pseudospectral method, the anelastic approximation, nor hyperviscosity (which is necessary to dissipate energy at high wavenumbers). The Athena simulation presented in this paper (cf., Fig. 6) has a computational domain of size H 3 0 resolved with 256 3 grid points. The background state is initially isothermal T 0 with uniform vertical gravity g 0 , which yields exponential hydrostatic density and pressure profiles. Units are chosen so that the gas pressure scale height H 0 = 1, gas density ρ 0 = 1 at the center of the box, and Keplerian angular speed Ω 0 = 1. In these units, we choose the gravitational acceleration g 0 = 10 and buoyancy frequency N 0 = 2. The isothermal atmosphere is perturbed with random velocity fluctuations, which are chosen to be incompressible and follow the Kolmogorov power spectrum E(k) ∝ k −5/3 for the whole range of k in the simulation. The overall normalization of the perturbation is chosen such that the initial kinetic energy is 1.9 × 10 −4 , with an rms velocity of 0.019 (i.e., 0.6% of the sound speed).
We use a third-order spatial reconstruction (although Athena is second-order both spatially and temporally), with the HLLC flux solver and CTU unsplit integrator. The timestep is calculated with a CFL number of 0.4. The boundary conditions are periodic in the y direction, and shearingperiodic in the x direction. We employ reflecting boundary condition in the z direction for all variables except the energy. We set the energy by extrapolating the pressure field within ghost zones in the z direction assuming hydrostatic balance. The implementation of the z boundary conditions using ghost zones is as follows: (i) the field data is copied from the ith grid cell to the ith ghost zone (where i varies from one to three); (ii) the sign of the vertical velocity is flipped; (iii) (2i − 1)∆zρg/(γ − 1) is added (subtracted) from the ghost zones below (above) the domain, where ∆z is the grid spacing in the z direction.
The simulation evolved for ∼64 orbits before showing signs of zombie vortices. We did several runs in the same atmosphere as above, but with a vortex as an initial condition as in MPJH13. At a resolution of 128 3 , the vortex formed critical layers, but no new vortices formed. When we increased the resolution to 256 in the radial direction (keeping a resolution of 128 in the vertical and azimuthal directions), the critical layers formed new vortices, which formed their own critical layers, etc. Thus, it appears that a radial resolution of ∆x = H 0 /128 is insufficient to form zombie vortices in Athena.

B. APPENDIX B -Determination of ∆
In MPJH13, we analytically showed that for a Boussinesq linear shear flow with shear rate σ, constant vertical gravity g 0 , and constant Brunt-Väisälä frequency N 0 , that the distance (in the cross-stream direction x) between a perturbation and the critical layer that it excites is δ ≡ N 0 /(k y |σ|), where k y is the wavenumber in the stream-wise direction y of the marginally-stable, linear, singular eigenmode (in which the singularity forms the critical layer). We assume spatial periodicity in the y direction, so k y ≡ 2π|m|/L y , where m is a non-zero integer. Thus, the distance between a perturbation and its resulting critical layer is: where ∆ is defined in Eq. (2-12). Our analytic calculations included both unbounded domains in x and bounded domains with rigid boundaries such that v x = 0 at the boundaries. Although we did not perform analytic calculations for critical layers with a shearing box boundary condition in x because the equations of motion are not autonomous in time, numerical calculations in MPJH13 show that Eq. (B1) appears to also be valid in a shearing box. Furthermore, using WKB theory, one can show that Eq. (B1) is also a good approximation in fully compressible and anelastic flows, and when vertical gravity and/or the Brunt-Väisälä frequency are functions of z. (See Paper II for details.) What physics selects the value of m in Eq. (B1), and, in particular, why do our simulations with shearing box boundary conditions always have |m| = 1 at late times? To answer these questions, we performed a series of simulations with different domain sizes, with both rigid boundary conditions in x and with shearing box boundary conditions. We observed that the critical layers that develop first are those with the smallest values of δ (or largest values of |m|) that are not so close to the perturbing source as to interfere with it -see MPJH13. However as time progressed, the critical layers (and their accompanying dipolar vortex sheets) with smaller values of |m| began to dominate the flow. Although calculations with shearing box boundaries were always eventually dominated by |m| = 1, we found that calculations with rigid boundary conditions in x sometimes had dominant critical layers with |m| > 1.
In our shearing box calculations, we had expected that δ would be nearly independent of the size of the computational domain and that the value of |m| would scale with domain size so as to keep δ nearly fixed. In the calculations with rigid boundary conditions, the critical layer that grew to dominate the early evolution of the flow was the one with the largest value of δ (or, smallest non-zero value of |m|) that "fit" into the computational domain, i.e., with δ ≤ L x , or equivalently with the smallest value of |m| such that |m| > ∆/L x . Why did the calculations with the shearing box give a different result? The answer is that the calculations did not give a different result! In the shearing box calculation, if the source that excites a critical layer is located at x = x 0 , then the symmetry associated with the shearing box creates an infinite number of "phony" image sources located at x = x 0 + nL x , for all integers n = 0. There is always some image source with an integer n such that the distance between that image source and some location within the computational domain is equal to ∆. That location is where the dominant critical layer appears in our calculations, and it always has |m| = 1. Therefore in the shearing box calculations, the critical layer that grew to dominate the early evolution of the flow was the one with the largest value of δ (which is δ = ∆) that "fit" into the computational domain, where now the definition of "fit" must be modified to include the possibility of an infinite number of image sources. This result concerning the "fit" is easily tested because it predicts that the location of the dominant critical layer should be at x = x 0 ± ∆ + nL x , with |m| = 1, and with a value of n such that 0 ≤ x < L x . The distance between the real source and critical layer should be equal to |nL x ± ∆|, rather than δ. This prediction has been numerically verified. It must be noted that unless n = 0, calculations with the shearing box boundary conditions are not physically relevant and are artifacts of the image sources created by the boundary conditions. To force n to be zero, we require ∆ < L x /2, or equivalently that β(L y /L x ) < 3π. Our choices of β and (L y /L x ) in §2.5 meet this requirement. Therefore, all of our initial-value calculations have |m| = 1 and δ = ∆.