Cooling Requirements for the Vertical Shear Instability in Protoplanetary Disks
aa r X i v : . [ a s t r o - ph . E P ] J u l Draft version October 14, 2018
Preprint typeset using L A TEX style emulateapj v. 5/2/11
COOLING REQUIREMENTS FOR THE VERTICAL SHEAR INSTABILITY IN PROTOPLANETARY DISKS
Min-Kai Lin & Andrew N Youdin Department of Astronomy and Steward Observatory,University of Arizona, 933 North Cherry Avenue, Tucson, AZ 85721, USA
Draft version October 14, 2018
ABSTRACTThe vertical shear instability (VSI) offers a potential hydrodynamic mechanism for angular momen-tum transport in protoplanetary disks (PPDs). The VSI is driven by a weak vertical gradient in thedisk’s orbital motion, but must overcome vertical buoyancy, a strongly stabilizing influence in colddisks, where heating is dominated by external irradiation. Rapid radiative cooling reduces the effec-tive buoyancy and allows the VSI to operate. We quantify the cooling timescale t c needed for efficientVSI growth, through a linear analysis of the VSI with cooling in vertically global, radially local diskmodels. We find the VSI is most vigorous for rapid cooling with t c < Ω − h | q | / ( γ −
1) in terms ofthe Keplerian orbital frequency, Ω K ; the disk’s aspect-ratio, h ≪
1; the radial power-law temperaturegradient, q ; and the adiabatic index, γ . For longer t c , the VSI is much less effective because growthslows and shifts to smaller length scales, which are more prone to viscous or turbulent decay. Weapply our results to PPD models where t c is determined by the opacity of dust grains. We find thatthe VSI is most effective at intermediate radii, from ∼ ∼ ∼
30 local orbital periods. Growth is suppressed by long cooling times both in the opaqueinner disk and the optically thin outer disk. Reducing the dust opacity by a factor of 10 increasescooling times enough to quench the VSI at all disk radii. Thus the formation of solid protoplanets, asink for dust grains, can impede the VSI. INTRODUCTION
Understanding how disks transport mass and an-gular momentum underlies many problems in astro-physics, including star and planet formation (Armitage2011). The turbulence associated with many transportmechanisms is particularly important for dust evolutionand planetesimal formation (Youdin & Lithwick 2007;Chiang & Youdin 2010).Magneto-hydrodynamic (MHD) turbulencedriven by the magneto-rotational instability (MRI,Balbus & Hawley 1991) has long been the most promis-ing transport mechanism in low mass disks with weakself-gravity. However, many parts of protoplanetarydisks (PPDs) are cold, have low levels of ionization,and do not support the MRI (Blaes & Balbus 1994;Salmeron & Wardle 2003). Recent simulations suggestthat significant portions of PPDs fail to develop MHDturbulence (e.g. Simon et al. 2013; Lesur et al. 2014; Bai2015; Gressel et al. 2015).A purely hydrodynamic mechanism could circumventdifficulties with non-ideal MHD, but must overcome thestrong centrifugal stability imposed by the positive ra-dial specific angular momentum gradient in nearly Ke-plerian disks (Balbus et al. 1996). One possible routeto hydrodynamic turbulence is the vertical shear in-stability (VSI, Urpin & Brandenburg 1998; Urpin 2003;Nelson et al. 2013, hereafter N13). The basic mechanismof the VSI in disks was first identified in the context ofdifferentially rotating stars (Goldreich & Schubert 1967,hereafter GS67, Fricke 1968). The VSI arises when ver-tical shear, i.e. a variation in orbital motion along therotation axis, destabilizes inertial-gravity waves, which [email protected], [email protected] Steward Theory Fellow are oscillations with rotation and buoyancy as restoringforces. Vertical shear occurs wherever the disk is baro-clinic, i.e. when constant density and constant pressuresurfaces are misaligned. Baroclinicity, and thus verticalshear, is practically unavoidable in astrophysical disks,except at special locations like the midplane.To overcome centrifugal stabilization, the VSI triggersmotions which are vertically elongated and radially nar-row. Vertical elongation taps the free energy of the ver-tical shear (Umurhan et al. 2013), but is also subject tothe stabilizing effects of vertical buoyancy if the disk isstably stratified. To overcome vertical buoyancy, the VSIrequires a short cooling time (GS67; N13). Rapid radia-tive cooling, i.e. a short thermal relaxation timescale,brings a moving fluid element into thermal equilibriumwith its surroundings, thereby diminishing buoyancy. Anisothermal equation of state implies instantaneous ther-mal relaxation, and is the ideal context for studying theVSI (Urpin 2003; N13; McNally & Pessah 2014, here-after MP14; Barker & Latter 2015, hereafter BL15).Alternatively, vertical buoyancy can be eliminated bystrong internal heating, i.e. by the onset of convection,so that the disk becomes neutrally stratified in the ver-tical direction (N13; BL15). However, realistic PPDsshould be vertically stably stratified in the outer re-gions, beyond ∼ t c . Non-linearhydrodynamical simulations with a prescribed t c findthat VSI turbulence in stably stratified disks requiresrapid cooling with t c shorter than orbital timescales M-K. Lin, A. N. Youdin(N13). When the cooling time is short enough, the VSIcan drive moderately strong transport, with Reynoldsstresses of α ∼ − times the mean pressure (N13).Simulations with realistic radiative transfer, in lieu ofa fixed t c , find that the VSI in irradiated disks drivetransport with α ∼ − in a ∼ § § §
4. Section 5 contains ourmain analytic results, leading to the derivation of thecritical cooling time in § § §
7, with coolingtimes derived from dust opacities. We discuss caveatsand extensions in §
8. We conclude in § WHY MUST COOLING BE SO FAST?
In a stably stratified thin disk, the VSI requires a ther-mal timescale t c significantly shorter than the disk dy-namical timescale, t c ≪ Ω − , (1)where Ω K is the Keplerian frequency (N13). This require-ment — that the cooling time be much shorter than anorbital period, which in turn is much smaller than therelevant oscillation period of vertically elongated grav-ity waves — is quite stringent. It highlights the fact thatvertical buoyancy is strongly stabilizing. Rapid radiativedamping is thus required to weaken buoyancy and allowthe weak vertical shear to drive instability.To roughly quantify the required smallness of t c , weconsider a vertically isothermal disk with aspect-ratio h ,radial temperature profile T ∝ r q and adiabatic index γ > h ∼ . q ∼ − . γ ≃ . h ≪ z , from the midplane as r ∂ Ω ∂z ≃ hq Ω K (cid:16) z H (cid:17) , (2)where H = hr is the characteristic disk scale height.This destabilizing shear competes with the stabilizingvertical buoyancy frequency, N z . In a thin disk, N z ≃ (cid:18) γ − γ (cid:19) (cid:16) zH (cid:17) Ω . (3) Fig. 1 maps the vertical shear rate and buoyancy fre-quency (as well as the gas density) in a fiducial PPDmodel, see § § t c be for vertical shear to prevail? Westart with the Richardson number Ri = N z / ( r∂ z Ω) ,a ratio of buoyant to shear energies. Though not pre-cisely applicable to our problem, non-rotating shear flowsare stable if Ri > / t c | r∂ z Ω | ≤ t c . | r∂ Ω /∂z | N z . (4)If we interpret Eq. 4 as a local criterion, then we seethat the maximum cooling time which permits instabilitydecreases with height as 1 / | z | ; the stabilizing effect ofvertical buoyancy increases more rapidly away from themidplane than the destabilizing effect of vertical shear.This finding is relevant to the radiative damping of socalled ‘surface modes’, which we consider in § t c < h | q | γ − − (5)for VSI growth that we derive in § | z | = γH/ § ∼
10 – 150AU. The innerdisk is complicated by the fact that the optically thickmidplane has longer cooling times. Thus the cooling cri-terion fails in the midplane but is satisfied above. Thiscase requires the more detailed analysis of § k x = 10 /H and 50 /H . In the inner disk, the higherwavenumber mode cools faster, due to radiative diffusion.Thus in the inner disk, shorter wavelength VSI modes aremore likely to grow. In the outer disk, the cooling timeis the same for different k x and z , because optically thincooling is independent of both lengthscale and density.Thus the outer barrier to growth, beyond ∼
150 AU inthis model, is independent of wavelength. This exam-ple highlights the fact that optically thin cooling sets alower limit to the cooling time; a fact that is sometimesertical shear instability 3
Figure 1.
The radial and vertical structure of a ‘minimum mass’ PPD model showing gas density (left, relative to the midplane densityat 1 AU), vertical shear rate (middle) and buoyancy frequency (right). z / H -3 -2 -1 Cooling times in Ω K-1 for k x H=10 t c Ω K < β crit z / H -3 -2 -1 Cooling times in Ω K-1 for k x H=50 t c Ω K < β crit Figure 2.
Estimate of the thermal timescale t c in a fiducial PPDmodel, the Minimum Mass Solar Nebula. Here, β crit = h | q | / ( γ − obscured when the radiative diffusion approximation ismade at the outset. GOVERNING EQUATIONS AND DISK MODELS
An inviscid, non-self-gravitating disk orbiting a cen-tral star of mass M ∗ obeys the three-dimensional fluidequations: ∂ρ∂t + ∇ · ( ρ v ) = 0 , (6) ∂ v ∂t + v · ∇ v = − ρ ∇ P − ∇ Φ ∗ , (7) ∂P∂t + v · ∇ P = − γP ∇ · v − Λ , (8)where ρ is the mass density, v is the velocity field (therotation frequency being Ω = v φ /r ), P is the pressure, γ is the constant adiabatic index, and Φ ∗ = − GM ∗ ( r + z ) − / is the gravitational potential of the central star with G as the gravitational constant. The cylindrical co-ordinates ( r, φ, z ) are centered on the star. In the energyequation (Eq. 8) the sink term Λ includes non-adiabaticcooling and (negative) heating. The gas temperature T obeys the ideal gas equation of state, P = R µ ρT, (9)where R is the gas constant and µ is the mean molecularweight. Thermal relaxation
We model radiative effects as thermal relaxation witha timescale t c , i.e., Λ ≡ ρ R µ T − T eq t c , (10)where T eq = T eq ( r, z ) is the equilibrium temperature. Wedefine the dimensionless cooling time β ≡ Ω K t c (11)relative to the Keplerian frequency, Ω K = p GM ∗ /r .We use the terms ‘thermal relaxation’ and ‘cooling’ in-terchangeably throughout this paper.For most of this paper we will take β as a constantinput parameter, so that we can vary the thermodynamicresponse of the disk from isothermal ( β ≪
1) to adiabatic( β ≫
1) in a controlled manner. However, in § β in PPDs with a dustopacity (as in Fig. 2). Baroclinic disk equilibria
The equilibrium disk is steady and axisymmetric withdensity, pressure and rotation profiles ρ ( r, z ), P ( r, z ) andΩ( r, z ), respectively. These functions satisfy0 = 1 ρ ∂P∂z + ∂ Φ ∗ ∂z , (12a) r Ω = 1 ρ ∂P∂r + ∂ Φ ∗ ∂r . (12b) M-K. Lin, A. N. YoudinA unique solution requires additional assumptions aboutdisk structure. We choose ρ ( r, ≡ ρ ( r ) = ρ (cid:18) rr (cid:19) p , (13a) T ( r, ≡ T ( r ) = T (cid:18) rr (cid:19) q , (13b) P ( r, z ) = K ( r ) ρ Γ ( r, z ) . (13c)where r is a fiducial radius and p and q are the standardpower-law indices for midplane density and temperature,respectively. Without loss of generality we take q <
0, asis typical in PPDs. The polytropic index Γ parametrizesthe vertical stratification with Γ = 1 (Γ = γ ) describingvertically isothermal (adiabatic) disks, respectively. Theideal gas law requires K = ρ − Γ0 R T /µ .We further define a modified sound speed c s ≡ p Γ P/ρ, (14)which in general differs from the isothermal, c iso = c s / √ Γ, and adiabatic c ad = c s p γ/ Γ, sound speeds.We also introduce the characteristic scale-height, H = c s ( r, / Ω K and disk aspect-ratio h ( r ) ≡ Hr ∝ r (1+ q ) / (15)We are interested in thin disks with h ≪ Density structure
The equilibrium density field follows from vertical hy-drostatic equilibrium (Eq. 12a). The solution for Γ = 1is (cid:18) ρρ (cid:19) Γ − = 1 + (Γ − h p z /r − ! . (16)Note that for Γ > h there is a disk surface H s where ρ ( r, H s ) = 0; and H s ≃ p / (Γ − H for h ≪ → ρ ( r, z ) ρ ( r ) = exp " h p z /r − ! , (17a) ≃ exp (cid:18) − z H (cid:19) for | z | ≪ r. (17b)Eq. 17b is the approximate form of the density field inthe thin-disk limit. We will primarily focus on verticallyisothermal disks, the relevant case for passively irradi-ated PPDs (Chiang & Goldreich 1997). Rotation and vertical shear profiles
The equilibrium rotation field, Ω( r, z ) follows from thedensity field and centrifugal balance (Eq. 12b), givingΩ ( r, z )Ω ( r ) = 1 + p + q Γ h ( r ) + s Γ − p z /r ! , (18) for all Γ, where s ≡ q + p (1 − Γ). We also refer to theepicyclic frequency κ defined via κ ≡ r ∂∂r (cid:0) r Ω (cid:1) = 1 r ∂j ∂r , (19)where j is the specific angular momentum.The vertical shear rate follows from Eq. 18, r ∂ Ω ∂z = sz Γ r Ω (1 + z /r ) / . (20)For Γ = 1, notice that s/ Γ = q is the only disk parameterthat affects vertical shear. Vertical shear increases lin-early with height near the midplane, and | ∂ z Ω | is max-imized at min( r/ √ , H s ), which is expected to limit themaximum VSI growth rate.A more general expression for vertical shear holds forvertical polytropes (which satisfy Eq. 13c), with no as-sumed radial structure: r ∂ Ω ∂z = − ρ Γ − ∂ ln ρ∂z dKdr = − g z Γ d ln Kdr , (21)where g z = ∂ z P/ρ . Entropy gradients and vertical buoyancy
The gradients of specific entropy, S ≡ C P ln P /γ /ρ , inour disk models are ∂S∂r = C P (cid:20) sγr + (cid:18) Γ γ − (cid:19) ∂ ln ρ∂r (cid:21) , (22a) ∂S∂z = C P (cid:18) Γ γ − (cid:19) ∂ ln ρ∂z = C P g z c s Γ − γγ , (22b)where C P is the heat capacity at constant pressure.The vertical buoyancy frequency N z is N z ≡ − C − P g z ∂S∂z = g z c s γ − Γ γ . (23)We only consider convectively stable disks with with Γ <γ and N z ≥ Stability without cooling
In the absence of cooling, axisymmetric stability is en-sured if both Solberg-Hoiland criteria are satisfied: κ − ρC P ∇ P · ∇ S > , (24a) − ∂P∂z (cid:18) ∂j ∂r ∂S∂z − ∂j ∂z ∂S∂r (cid:19) > z ), we consider z > ∂ z P < ∂ r j ≃ r Ω and using Eq. 21 for ∂ z j we find that γ − Γ > h Γ (cid:18) ρρ (cid:19) Γ − (cid:20) s − s ( γ − Γ) ∂ ln ρ∂ ln r (cid:21) (25)implies stability. For typical model parameters, theright hand side is O ( h ) . − . The left hand sideis positive and order unity in our disk models (withertical shear instability 5 γ − Γ ≥ .
1) implying strong stability to convection, seeEq. 23. This stable stratification is expected for irradi-ated disks (Chiang & Goldreich 1997). Thus, adiabaticdisturbances in a standard, irradiated disk are stable toa disk’s vertical shear, explaining why the VSI requirescooling. LINEAR PROBLEM
This section presents the two sets of equations we useto study the linear development of the VSI. Both sets areradially local and vertically global with a finite coolingtime. The first set, presented in § § Radially local approximation
We consider axisymmetric perturbations to the abovedisk equilibria. The growth of the VSI is strongest forshort radial wavelengths, as compared to the disk ra-dius (N13; BL15). We thus perform a standard two stepprocess to obtain linearized equations in the radially lo-cal approximation (e.g., Goldreich & Schubert 1967, whoalso consider vertically localized perturbations).We first expand all fluid variables into the equilibriumvalue plus an Eulerian perturbation denoted by δ , e.g. δρ for the perturbed density, and drop all non-linearperturbations. Second, we perform a Taylor expansionabout a fiducial radius r with the local radial coordi-nate x = r − r , keeping only the leading order terms in x/r . We also relabel (trivially for axisymmetric pertur-bations) the azimuthal direction as y .Perturbations take the form of a radial plane wave witharbitrary vertical dependence, e.g. δρ → δρ ( z ) exp (i k x x − i υt ) , (26)where υ = ω + i σ is a (generally) complex frequency,with ω and σ being the real frequency and growth rate,respectively, and k x is a real radial wavenumber. Wetake k x > k x r ≫ ∂ z Ω, refer toequilibrium values at ( r , z ).We further define the perturbation variables W ≡ δP/ρ and Q ≡ c s δρ/ρ . With this procedure the lin-earized system of equations arei υ Qc s = i k x δv x + dδv z dz + ∂ ln ρ∂z δv z + ζ ∂ ln ρ∂r δv x , (27a)i υδv x = i k x W − δv y − ζ ρ ∂P∂r Qc s (27b)i υδv y = r ∂ Ω ∂z δv z + κ δv x , (27c)i υδv z = dWdz + ∂ ln ρ∂z ( W − Q ) , (27d)i υW = c s ∂ ln ρ∂z δv z + c s γ Γ (cid:18) i k x δv x + dδv z dz (cid:19) + 1 t c (cid:18) W − Q Γ (cid:19) + ζ ρ ∂P∂r δv x . (27e) This is a set of ordinary differential equations (ODEs) in z . Solutions to these are presented in § §
7. Thecoefficient ζ = 1 is introduced simply to label terms withexplicit radial gradients of the equilibrium state, whichagain are evaluated at ( r , z ). For clarity, hereafter wedrop the subscript 0. We will consider the effects of ig-noring these radial gradients below. Reduced equation for vertically isothermal disks
Our simplified model starts with Eqs. 27 and makesthe following additional simplifications:1. We set Γ = 1, focusing on vertically isothermaldisks.2. We set ζ = 0, neglecting terms with an explicit de-pendence on the radial structure of the equilibriumdisk. Vertical shear, which implicitly depends onthe radial temperature gradient, is retained. Thisfully-radially-local approximation is also made inthe ‘vertically global shearing box’ of MP14.3. We make the low frequency approximation, as-suming | υ | ≪ κ , Ω . Similar to the incom-pressible (GS67) or anelastic (N13; BL15) approx-imations, the low frequency approximation filtersacoustic waves in favor of inertial-gravity waves(Lubow & Pringle 1993).4. We make the Keplerian approximation, settingΩ → Ω K and κ → Ω K , but retaining the verti-cal dependence in the crucial vertical shear term, ∂ z Ω.5. We consider thin disks with h ≪ z = z/H, ˆ k = k x H, ˆ υ = ˆ ω + iˆ σ = υ/ Ω K , (28)where ˆ ω and ˆ σ are real, the above approximations leadto a single second order ODE, δv ′′ z − zAδv ′ z + ( B − C ˆ z ) δv z = 0 , (29)where ′ denotes d/d ˆ z and A ≡ hq ˆ k, (30a) B ≡ ˆ υ (cid:16) χ + ˆ k (cid:17) − (cid:16) χ + i hq ˆ k (cid:17) , (30b) C ≡ (1 − χ ) (cid:16) ˆ k − i hq ˆ k (cid:17) , (30c)where χ = 1 − iˆ υβ − iˆ υβγ . (31)The derivation of Eq. 29 is detailed in Appendix A.In Appendix B we discuss the limits of these approx-imations. In particular we point out that the fully-radially-local approximation is only valid for short cool-ing times, β . O (1), which is the regime in which the VSIoperates, as demonstrated below. On the other hand, for β & O (1), this approximation can introduce artificial in-stability due to the non-self-consistent neglect of globalradial gradients (see § B.1). M-K. Lin, A. N. Youdin ANALYTIC RESULTS WITH FINITE COOLING TIMES
Previous analytic studies of the VSI have largely fo-cused on isothermal perturbations, with infinitely rapidcooling, as discussed in the introduction. We further ex-plore this β = 0 limit in Appendix C both to furtherdevelop intuition for this idealized case and to establisha connection with previous works. However our maininterest is the effect of finite cooling times, which we ex-plore analytically in this section.In § β = 0, slows the growthrate of the VSI. In § § Introducing a small but finite cooling time
We are interested in finite thermal relaxationtimescales β >
0, but it is instructive to first ask themore analytically tractable question: how do the eigen-frequencies and eigenfunctions change when we change β from zero to a small but finite value? For sufficientlysmall β we expect the solution to only differ slightly froma case with β ≡
0. We thus perturb a solution for β = 0to see the effect of finite cooling.For definiteness, let us consider the simple solution β ≡ , δv z = 1 , ˆ υ = 1 + i hq ˆ k k , (32)which solves Eq. 29 since δv z =constant and χ = 1 sothat B = C = 0. This is the lowest order VSI mode, the‘fundamental corrugation mode’, where the vertical ve-locity is constant throughout the disk column. Hereafterwe shall simply call it the fundamental mode.The fundamental mode has been observed to dominatenumerical simulations (Stoll & Kley 2014, N13), and arethe ones we find to typically dominate in numerical cal-culations with increasing β ( § k , as well as in PPDs with a realistic estimate for ther-mal timescales ( § β ≡ β → δβ, δv z → δv z + δv z , ˆ υ → ˆ υ + δ ˆ υ, (33)which implies χ → δχ = 1 + iˆ υ ( γ − δβ, (34)with δv z and ˆ υ given by Eq. 32. We may then seek δv z = a ˆ z + b, (35)where a , b are constants.We insert Eq. 33 into Eq. 29, keeping only first orderterms, and solve for δ ˆ υ using the above expressions for δχ and ˆ υ . We find imaginary part of δ ˆ υ is δ ˆ σ = − ( γ −
1) ˆ k (cid:16) ˆ k − h q ˆ k − h q (cid:17) (cid:16) h q ˆ k (cid:17) (cid:16) k (cid:17) δβ. (36) Since h ≪
1, introducing finite cooling δβ > δ ˆ σ <
0, i.e. stabilization, since γ > z dependence in δv z makes sense because it becomessignificant at large | ˆ z | , i.e. away from the midplane wherethe effect of buoyancy first appears as β is increased fromzero (since N z ∝ z for a thin disk, see Eq. 23). Explicit solutions and dispersion relation
We now solve Eq. 29 explicitly. We first write δv z (ˆ z ) = g (ˆ z ) exp (cid:18) λ ˆ z (cid:19) , (37)where λ is a constant to be chosen for convenience. In-serting Eq. 37 into Eq. 29 gives0 = g ′′ − ˆ z ( A − λ ) g ′ + ( B + λ ) g + (cid:0) λ − λA − C (cid:1) ˆ z g. (38)We choose λ to make the coefficient of ˆ z g vanish, andimpose the vertical kinetic energy density ρ | δv z | ∝| g | exp (Re λ − /
2) ˆ z to remain finite as | ˆ z | → ∞ .Then assuming g (ˆ z ) is a polynomial, we requireRe λ < , (39)which amounts to choosing λ = 12 (cid:16) A − p A + 4 C (cid:17) . (40)Eq. 29 becomes0 = g ′′ − ˆ z ( A − λ ) g ′ + ( B + λ ) g. (41)We seek polynomial solutions g (ˆ z ) = M X m =0 b m ˆ z m , (42)which requires B (ˆ υ ) + λ (ˆ υ ) = M [ A − λ (ˆ υ )] . (43)For ease of analysis, we re-arrange Eq. 43 and square itto obtain a polynomial in ˆ υ ,0 = X l =0 c l ˆ υ l , (44)where the coefficients c l are given in Appendix D.1. Thedispersion relation ˆ υ = ˆ υ (ˆ k ) is complicated and gener-ally requires a numerical solution. However, simple re-sults may be obtained in limiting cases which we considerbelow. The critical cooling time derived
Here we estimate the maximum thermal relaxationtimescale β c that allows growth of the VSI. In AppendixD.2 we derive the following relation between β c and thewave frequency ˆ ω c , assuming marginal stability to theertical shear instability 7VSI: (cid:16) ˆ ω c ˆ k (cid:17) = (cid:16) ˆ ω c ˆ k (cid:17) + M ( M + 1) (cid:16) − h q ˆ k (cid:17) , (45a) (cid:16) ˆ ω c ˆ k (cid:17) = β c hq ( γ − M ) (cid:16) ˆ ω c ˆ k (cid:17) (45b) − M ( M + 1) . These equations consider ˆ k ≫ | ˆ ω ˆ k | is O (1) and β c ≪
1. Note that Eq. 45a reduces to thedispersion relation for low-frequency inertial waves in theabsence of vertical shear (see Appendix C.2.1 for details).We are most interested in the longest cooling timewhich allows growth for any M . In Appendix D.3, wefind that if | hq ˆ k | ≤
1, so that marginal stability ex-ists, then β c decreases with increasing M . The VSIcriterion is then set by β c for the fundamental mode( M = 0), for which the exact solution to Eq. 45a —45b is ˆ ω c ˆ k = β c ( γ − /hq = −
1. We thus find the cool-ing criterion for VSI growth is β < β crit ≡ h | q | γ − . (46)The thin disk approximation, h ≪
1, indicates that β crit ≪
1, as assumed to obtain Eq. 45a—45b. Weheuristically explain β crit in §
2, and numerically test itsvalidity in § NUMERICAL RESULTS FOR FIXED COOLING TIMES
This section investigates the linear growth of the VSIwith fixed cooling times, β Ω − , that are independent ofheight. We solve the radially local model of Eqs. 27 with ζ = 1, i.e. retaining radial gradients of the backgrounddisk structure. The equilibrium disk structure given is in § § T l up to l = 512and discretizing the equations on a grid with z ∈ [ − z max , z max ]. Our standard boundary conditions at thevertical boundaries is a free surface,∆ P ≡ δP + ξ · ∇ P = 0 at z = ± z max , (47)where ξ is the Lagrangian displacement with meridionalcomponents ξ x,z = i δv x,z /υ . In some cases we impose arigid boundary with δv z ( ± z max ) = 0.The above discretization procedure converts the linearsystem of differential equations to a set of algebraic equa-tions, which we solve using LAPACK matrix routines .Unless otherwise stated, we adopt fiducial disk parame-ters ( γ, Γ) = (1 . , . p, q, h ) = ( − . , − , . z max = 5 H . This setup is effectivelyvertically isothermal since c s ( z max ) ≃ . c s (0). In Fig.3 we plot the corresponding equilibrium rotation profile,which shows that Ω decreases slightly away from the mid-plane. From Eq. 46, we expect that β < β crit = 0 .
125 isneeded for VSI growth in our fiducial model. Available at . While this model has a zero density surface at z = ± H s ≃± . H (see Eq. 18), the surface is outside the numerical domain. -4 -2 0 2 4z/H0.9850.9900.9951.000 Ω / Ω ( ) Figure 3.
Equilibrium rotation profile Ω( z ), normalized by itsmid-plane value, for the fiducial disk model with Γ = 1 .
011 and( p, q, h ) = ( − . , − , . ω /(h Ω K )0.00.20.40.60.81.01.21.4 σ / ( h Ω K ) f z max =5H, ∆ P=0z max =5H, δ v z =0z max =3H, ∆ P=0z max =7H, ∆ P=0k x H=10, β =10 -3 Analytic ( β =0)Analytic ( β =10 -3 ) Figure 4.
Growth rate, σ , vs. oscillation frequency, ω , for VSImodes with wavenumber k x = ˆ k/H = 10 /H in the fiducial diskmodel with rapid cooling, β = 10 − . We investigate the effect ofvertical box size, from z max = 3 H ( blue triangles ) to 7 H ( green x’s) with a free surface boundary (∆ P = 0). For z max = 5 H (black dots) we also compare free surface ( black dots ) and rigid( red plusses ), finding little difference. The fundamental mode ismarked by ‘f’. Lines show the analytic predictions for isothermalperturbations ( β = 0, dashed ) and for β = 10 − ( dotted ). See textfor discussion. Rapid thermal relaxation
We first calculate the VSI in a disk with rapid thermalrelaxation by setting β = 10 − . Since β ≪ β crit , this casegives similar results to the well studied case of isothermalperturbations (e.g. N13; MP14; BL15). This case allowsus to explore and test our finite cooling time model in afamiliar context.
Case study: ˆ k = 10 Fig. 4 maps the growth rates σ and wave frequencies ω of unstable modes for ˆ k = 10. Each symbol correspondsto a different order mode, i.e. a different number of ver-tical node crossings. Lower order modes have smaller | ω | for fixed radial wavenumber, as expected for inertial-gravity waves (see § C.2.1).We compare our numerical results to the analytic dis-persion relations for β = 0 and 10 − , from Eqs. C16 and44, respectively. The analytic models also have discretemodes, which lie on the continuous curves that are plot-ted for clarity in Fig. 4. While the analytic models havemany simplifications, they crucially lack an artificial ver-tical surface (as they assume an infinite vertical domain). We confirmed that smaller β values gave similar results. M-K. Lin, A. N. Youdin ω = -2.04h Ω K , σ =0.45h Ω K , β =0.001, k x H=10.0 -1.0-0.50.00.51.0 δ P / ρ K K x -4 -2 0 2 4z/H-0.4-0.20.00.20.40.60.81.0 δ v z Figure 5.
Pseudo-enthalpy perturbation W (top) and verticalvelocity perturbation δv z of the fundamental VSI with ˆ k = 10 inthe fiducial disk with rapid thermal relaxation β = 10 − . The real(imaginary) parts of W and δv z are plotted as solid (dotted) lines. k x H=10, β =10 -3 , ρ ( δ v x , δ v z ), δ P -0.3 -0.2 -0.1 0.0 0.1 0.2 0.3x/H-4-2024 z / H -1.00-0.67-0.33 0.00 0.33 0.67 1.00 Figure 6.
Visualization of the fundamental VSI mode in Fig. 5.For this plot the mode has been transformed back to real space.Thus, the color scale corresponds to the real pressure perturbationscaled by its maximum value. Arrows show the real meridionalflow multiplied by √ ρ . The horizontal axis has been stretched forclarity. For our fiducial case (black dots), the lower order andlower frequency modes closely follow the analytic pre-diction with growth rates increasing with | ω | . However,for | ω | & h Ω K the growth rate declines with | ω | . Thisbreak from the analytic prediction is not due the choiceof boundary condition, as we demonstrate by considering ω /(h Ω K )0.00.51.01.52.0 σ / ( h Ω K ) f ss sss z max =5H, ∆ P=0z max =5H, δ v z =0z max =3H, ∆ P=0z max =7H, ∆ P=0k x H=30, β =10 -3 Analytic ( β =0)Analytic ( β =10 -3 ) Figure 7.
Same as Fig. 4 but for ˆ k = 30. Examples of surfacemodes are marked by ‘s’. rigid boundaries (red plusses).Rather, the decline in growth rate for high order mod-els is due to the vertical box size, as seen by compar-ing the z max /H = 3 , z max give better agreement with the analytic theory,which have z max → ∞ . Larger boxes include regions oflarger vertical shear (see Eqs. 2 and 20), which higherorder modes can tap to give larger growth rates. Forthe z max = 7 H case, we begin to see another branch ofmodes with the highest growth rates at ω ∼ − h Ω K .This branch contains the ‘surface modes’ described be-low.The need for vertically extended domains to capturethe largest VSI growth rates is problematic, especiallyfor hydrodynamic simulations. This complication is mit-igated by at least two factors. First, we show below thatwith slower cooling, the growth of higher order modesis preferentially damped, consistent with the analysis in § H . Quantifying which modeswill contribute most to non-linear transport is an impor-tant issue, but beyond the scope of this work.Fig. 5 shows eigenfunctions W = δP/ρ and δv z forthe lowest order fundamental mode. (In this and similarplots below, we normalize eigenfunctions such that itsmaximum amplitude is unity with vanishing imaginarypart at the lower boundary.) Fig. 6 maps the pressureperturbation and meridional flow, scaled by √ ρ to reflectthe contribution to kinetic energy. Notice the stretched x axis. Radial velocities are in fact typically much smallerthan vertical velocities, as expected for a vertically elon-gated, anelastic mode (N13). Most of the kinetic energyis contained within ∼ H of the midplane due to thedensity stratification. Surface modes of the VSI
The surface modes mentioned above are more promi-nent for larger wavenumbers. Fig. 7 maps the ˆ k = 30eigenvalues, with the surface modes labeled. Surfacemodes strongly depend on the location of the imposedvertical boundary, and disagree with the analytic mod-els, which lack an imposed surface. We thus discount thephysical significance of surface modes for our model, aswe explain further below.Surface modes are a well known feature of the VSI infinite vertical domains (N13; MP14). BL15 show thatsurface modes arise whether the surface is imposed (asertical shear instability 9 ω = -1.58h Ω K , σ =1.48h Ω K , β =0.001, k x H= 30 -4 -2 0 2 4z/H-1.0-0.50.00.51.0 δ v z Figure 8.
Example of a ‘surface mode’ in the fiducial disk model. in our models) or natural (for disk models with a zero-density surface). BL15 mention the interface between adisk’s interior and corona as a physical boundary thatcan trigger surface modes. Since vertically isothermaldisk models lack a surface, modes which depend on theadopted value of z max are not physically meaningful.The possibility of artificial surface modes seedinggrowth in numerical simulations of the VSI merits fur-ther study. Indeed N13 find that the initial growth ofperturbations primarily occurs near the vertical bound-aries. The motion in surface modes is indeed concen-trated near the surface, as shown in Fig. 8, where thedensity is low. Thus, their contribution to transportmight (by themselves) be weak, following the argumentsin § § Slower thermal relaxation
We now consider the effect of longer cooling timescalesby gradually increasing β . We expect VSI growth for β < β crit = 0 . β values, as argued in § k = 10 and β from 0 .
01 to 0 .
1. Theanalytic results from Eq. 44 are now plotted as differentcurves for each mode order M , with β varying along eachcurve. We see the standard increase in frequency, | ω | ,with mode order.As expected, the higher order modes are preferentiallydamped. For β ≥ .
03 the fundamental ( M = 0) mode isthe fastest growing. For β = 0 .
1, only the M = 0 modegrows. This growth is slow since β is near β crit .Fig. 9 also shows how cooling times affect the de-pendence on z max , the size of the vertical domain. Forall cooling times, the fundamental mode depends onlyweakly on z max , and there is good agreement with an-alytic values. This convergence is reassuring given theimportance of the fundamental mode at longer coolingtimes. For higher order modes and short cooling times,however, the eigenvalues vary strongly with z max . Thedisagreement with analytic theory is strongest for thesmallest domain, with z max = 3 H . For β ≥ .
03, the z max /H = 5 , σ / ( h Ω K ) β =0.01 β =0.03 β =0.05 β =0.10k x H=10, z max =3HM=0, analyticM=1, analyticM=2, analyticM=3, analytic0.00.20.40.60.8 σ / ( h Ω K ) β =0.01 β =0.03 β =0.05 β =0.10k x H=10, z max =5HM=0, analyticM=1, analyticM=2, analyticM=3, analytic0 2 4 6 8 10 12- ω /(h Ω K )0.00.20.40.60.8 σ / ( h Ω K ) β =0.01 β =0.03 β =0.05 β =0.10k x H=10, z max =7HM=0, analyticM=1, analyticM=2, analyticM=3, analytic
Figure 9.
Unstable modes with ˆ k = 10 and thermal relaxationtimescales β = 0 .
01 (black dots), β = 0 .
03 (red crosses), β = 0 . β = 0 . z max = 3 H (top), z max = 5 H (middle, the fiducialsetup) and z max = 7 H (bottom). Lines are computed from Eq.44 for M = 0 (solid, fundamental mode) and M = 1 , , β increases con-tinuously from 0 .
01 to 0 . β values with numerical results. (aided by the fact that only M < k = 10 and β = 0 .
1. Com-pared with the small β case in Figs. 5 and 6, the eigen-fuctions show a more complex variation of phase withheight. This dependence yields the ‘tilted’ appearance ofpressure field in Fig. 11. Physically, the increased roleof buoyancy, which increases with height, explains whylarger β values produce more complex vertical structure.Fig. 12 shows how cooling affects higher wavenumbermodes, specifically ˆ k = 30. We only show the largerdomains with z max = 5 H, H . These cases again showgood convergence for β ≥ .
03, which the smaller domainwith z max = 3 H lacks.With slower cooling, the wave frequency shows a morecomplex dependence on mode order, both analyticallyand numerically. For β = 0 .
1, the fundamental mode(which lies on the M = 0 curve) no longer has the small-est | ω | value.Moreover, the fundamental mode is no longer the0 M-K. Lin, A. N. Youdin ω = -1.97h Ω k , σ =0.10h Ω k , β =0.100, k x H=10.0 -1.0-0.50.00.51.0 δ P / ρ k k x -4 -2 0 2 4z/H-1.0-0.50.00.51.0 δ v z Figure 10.
Same as Fig. 5 but with a thermal relaxationtimescale β = 0 . k x H=10, β =10 -1 , ρ ( δ v x , δ v z ), δ P -0.3 -0.2 -0.1 0.0 0.1 0.2 0.3x/H-4-2024 z / H -1.00-0.67-0.33 0.00 0.33 0.67 1.00 Figure 11.
Same as Fig. 6 but with a thermal relaxationtimescale β = 0 . fastest growing mode for β ≥ .
03, or even for β = 0 . | hq ˆ k | =1 . § β ≥ .
03. We are thus confident that ar-tificial surface modes should not affect our determinationof the critical cooling time. σ / ( h Ω K ) s β =0.01 β =0.03 β =0.05 β =0.10k x H=30, z max =5HM=0, analyticM=1, analyticM=2, analytic0 1 2 3 4- ω /(h Ω K )0.00.20.40.60.81.01.2 σ / ( h Ω K ) s s β =0.01 β =0.03 β =0.05 β =0.10k x H=30, z max =7HM=0, analyticM=1, analyticM=2, analytic
Figure 12.
Same as Fig. 9 but for ˆ k = 30 and domain sizes z max = 5 H ( top ) and 7 H ( bottom ). Examples of surface modesare marked with ‘s’. -3 -2 -1 M a x . g r o w t h r a t e i n h Ω K β crit k x H=1k x H=5k x H=10k x H=30k x H=50k x H=100z max =5H10 -3 -2 -1 β M a x . g r o w t h r a t e i n h Ω K β crit k x H=1k x H=5k x H=10k x H=30k x H=50k x H=100z max =7H
Figure 13.
Maximum VSI growth rates in the fiducial disk modelas a function of the thermal relaxation timescale β , for z max = 5 H (top) and z max = 7 H (bottom). The vertical line is the criticalthermal timescale β crit obtained from Eq. 46. Critical thermal relaxation timescale
Having explored the behavior of VSI modes with cool-ing, we turn to the numerical validity of the analytic cool-ing criterion for vertically isothermal disks, β < β crit = h | q | / ( γ − β crit = 0 . k = 5 ,
10, the growth rate drops to zero at theexpected β = β crit . For longer wavelength modes, withˆ k = 1, growth persists to slightly longer cooling times.This difference is not surprising since our analytic deriva-tion assumed ˆ k ≫
1, but the change is quantitativelyminor.For shorter wavelengths, with ˆ k ≥
30, growth persistsfor β significantly larger than β crit . This tail of growth ispartly explained by the breaking of the ˆ k < / | qh | = 20approximation, used in the analytic derivation. Despitethe lack of a clear stability boundary at high ˆ k , the β crit threshold remains useful, since growth rates dropto .
10% of their maximum value at β crit and continueto fall for larger β . Moreover, we expect that longerwavelength modes with ˆ k .
20 are more significant fordisk transport, see § z max /H = 5 , β → β crit = 0 .
125 is con-sistent with the simulations shown in their Fig. 12. N13found nonlinear VSI growth for β = 0 . < β crit but nogrowth for β = 0 . > β crit . Fig. 14 confirms that the numerically determined β crit shows the expected scaling with disk parameters h , q , and γ . For this test we fix ˆ k = 10 (an appropriate value forall the reasons discussed above) and measure the small-est β value where growth vanishes. The agreement withthe analytic scaling of Eq. 46 is quite good, confirmingthe applicability of our critical cooling time in verticallyisothermal disks. Vertically non-isothermal disk model
We now generalize our critical cooling time by con-sidering disks which are not vertically isothermal. Werely on physical arguments for this generalization. Indisks with weaker vertical stratification, i.e. Γ > γ , we expect β crit to be larger. We thus transform1 / ( γ − → / ( γ − Γ).We also expect β crit to scale with the vertical shear.In general the vertical shear rate ∂ z Ω ∝ s , the radialentropy gradient, see Eq. 20. We thus transform q → s = q + p (1 − Γ).Our estimate for the generalized cooling time crite-rion — valid for both vertically isothermal and non-isothermal disks (cf. Eq. 46) — is thus β crit → β crit , gen = h | s | γ − Γ . (48)We test β crit , gen using a disk model with ( p, q, h ) =( − . , − , .
05) and ( γ, Γ) = (1 . , . Since the dimensionless cooling time T relax in N13 is normalizedto the orbital period, we convert β = Ω K P orb T relax = 2 πT relax . β c r i t NumericalAnalytical1 2 3 4 5( γ -1) -1 β c r i t NumericalAnalytical0.02 0.04 0.06 0.08 0.10h0.000.050.100.150.200.250.30 β c r i t NumericalAnalytical
Figure 14.
Dependence of the upper limit to the thermal relax-ation timescale β crit for the fundamental VSI mode on disk pa-rameters. The fiducial setup is ( γ, Γ = (1 . , . p, q, h ) =( − . , − , . q ∈ [ − . , − . γ ∈ [1 . , . h ∈ [0 . , . k = 10. s = − .
85 and β crit , gen = 0 . ∼ β crit = 0 . z max = 0 . H s ≃ . H , see Eq.16. At this surface, the vertical buoyancy frequency, N z ,diverges as | z | → H s . While we might expect patholog-ical behavior for this model, we fortunately do not findit.In Fig. 15 we plot the mode diagram for ˆ k = 30 andseveral values of β . This plot is qualitatively similar tothe vertically isothermal case, Fig. 9. As before, larger β values rapidly stabilize higher order modes. For suf-ficiently large β only the fundamental mode is unsta-ble, with a growth rate that decreases as β approaches β crit , gen .Fig. 16 shows how maximum VSI growth rates varywith β in the vertically non-isothermal model. The qual-itative behavior is similar to Fig. 13 for the verticallyisothermal disk. For ˆ k = O (10), growth rates drop to2 M-K. Lin, A. N. Youdin ω /(h Ω K )0.000.050.100.150.200.25 σ / ( h Ω K ) β =0.01 β =0.05 β =0.1 β =0.2 β =0.4k x H=30
Figure 15.
Unstable modes in a vertically non-isothermal diskwith ˆ k = 30 and a range of thermal relaxation timescales. -3 -2 -1 β M a x . g r o w t h r a t e i n h Ω K β crit,gen β crit k x H=5k x H=10k x H=50k x H=100k x H=200k x H=300
Figure 16.
Maximum VSI growth rates in the vertically non-isothermal disk model as a function of the thermal relaxationtimescale t c = β Ω − . The vertical dashed line is the critical ther-mal timescale given by Eq. 46. A conjectured critical thermaltimescale accounting for a non-isothermal background, given by β crit , gen = 0 .
425 (see Eq. 48), is also shown as the vertical dottedline. zero at the predicted β crit , gen . As in the vertically isother-mal case, modes with both lower and higher wavenum-bers exhibit slow growth rates beyond the critical coolingtime. As argued in § β crit , gen should be agood approximation. This expectation remains to betested. THE VSI WITH REALISTIC PROTOPLANETARY DISKCOOLING
To this point we have treated the cooling time, β , as afree parameter that, moreover, is independent of height.To understand the applicability of the VSI to PPDs, wenow consider β values that are consistent with PPD mod-els. Section 7.1 develops our PPD cooling law, which de-pends on disk radius, disk height and perturbation wave-length. We compare these PPD β values to β crit in § β crit was derived assuming a vertically constant β .Thus in § β values. Cooling model
Radiative diffusion and Newtonian cooling
The thermal relaxation of a perturbation depends onthe relative sizes of the perturbation’s lengthscale, l , andthe photon mean-free-path, l rad ≡ κ d ρ (49)where κ d is the opacity, specifically the dust opacity ap-propriate for cold PPDs.In the optically thick limit, l ≫ l rad , radiative diffusionsmooths out thermal perturbations. In this regime, thelinearized cooling function is δ Λ diff = PρT C v ∇ · ( k rad ∇ δT ) , (50)where k rad is the radiative conduction coefficient definedbelow, and C v is the specific heat capacity at constantvolume.Since the VSI is characterized by vertically-elongated,radially-narrow disturbances, we retain only the radialderivatives of the perturbations in Eq. 50. Thus, in theradially local approximation, we have δ Λ diff ≃ − ηk x P δTT , (51)where the radiative diffusion coefficient is η ≡ k rad ρC v = 16 σ s T κ d ρ C v , (52)and σ s is the Stefan-Boltzmann constant. The thermalrelaxation timescale (defined in Eq. 10) for radiativediffusion is thus t diff = 1 ηk x . (53)In the optically thin regime, l ≪ l rad , thermal re-laxation operates by ‘Newtonian cooling’. The coolingtime is independent of l and inversely proportional tothe opacity, t thin ∝ /κ d . Specifically t thin = l η , (54)and t thin does not depend on ρ (because our adopted κ d depends on T only, see below).Our general cooling time for all perturbations, t c ≡ t thin + t diff , (55)is a simple prescription to smoothly match the opticallythick and thin limits. PPD cooling times
In a vertically isothermal disk with surface density Σ = ρ √ πH , the dimensionless thermal relaxation time β becomes β ( r, z ; ˆ k ) ≡ t c Ω K = Σ Ω K ηρ (cid:20) κ d Σ + ˆ ρ ( z )2 π ˆ k (cid:21) . (56)where ˆ ρ = ρ/ρ . The first and second terms in squarebrackets represent the optically thin and thick coolingregimes, respectively.ertical shear instability 13 -4 -2 0 2 4z/H0.010.10 β = t c Ω K β crit β crit r=50AUr=5AU Figure 17.
Thermal relaxation timescales in the fiducial MMSNat r = 50AU and r = 5AU for ˆ k = 30 (solid lines). The corre-sponding horizontal dashed lines are the critical thermal relaxationtimescales derived in linear theory. We adopt the Minimum Mass Solar Nebula (MMSN)disk model of Chiang & Youdin (2010) which specifiesΣ = 2200 r − / g cm − , (57a) T = 120 r − / K , (57b)with r AU ≡ r/ AU. This model has q = − / p = − /
14 and h = 0 . r / , using µ = 2 .
33 and the aboverelation between Σ and ρ ∝ r p .For the dust opacity we adopt the Bell & Lin (1994)law with κ d ∝ T , giving κ d = 2 . κ d r − / cm g − , (58)where the opacity normalization factor, ˆ κ d , scales withthe ratio of small dust to gas; ˆ κ d = 1 is our fiducial case.For this disk and opacity model, the cooling time be-comes β ( r, z ; ˆ k ) =3 . × − ˆ κ − d r / × h . × ˆ κ d ˆ ρ ( z ) r − / ˆ k − i . (59)Perturbations are in the optically thin regime forˆ k & . × ˆ κ d ˆ ρ ( z ) r − / . (60)In the inner disk, e.g. r AU ∼
1, only extremely small-scale perturbations are in the optically-thin regime nearthe midplane. However, in the outer disk, e.g. r AU ∼
10, more moderate wavenumbers ˆ k &
10 experienceoptically-thin cooling. β vs. β crit in PPDs The simplest way to estimate whether the VSI can op-erate at a given radius in a PPD is to compare the coolingtime (Eq. 59) to its critical value (Eq. 46), evaluated forthe MMSN as β crit = 0 . r / . (61)This comparison is strictly valid only for optically thincooling which is independent of height, as assumed in thederivation of β crit . This condition typically holds in theouter disk. Fig. 17 shows that a ˆ k = 30 perturbationat 50AU has vertically constant β . Furthermore, since β < β crit , the disturbance would be unstable. For optically thick cooling, β varies with height, com-plicating the comparison with β crit . At 5AU, Fig. 17shows that cooling times are long in the midplane with β > β crit . However, since β < β crit away from the mid-plane, we require the analysis of § β < β crit at all heights.In Fig. 18, we compare β to β crit across a range of diskradii for different heights, opacity values and wavenum-bers. For optically thin cooling in the outer disk, curvesfor different wavenumbers overlap, as expected from Eq.59. Since the (optically thin) slope of β is steeper thanfor β crit , VSI growth can be suppressed at large radii (forthe chosen opacity law). This effect is seen for ˆ κ d = 1in the central panels of Fig. 18, where VSI is dampedoutside ∼
150 AU.We find that the most important factor for VSI growthis the opacity. With a smaller opacity, ˆ κ d = 0 .
1, growthis suppressed at all radii, as shown in the top panels ofFig. 18. Since optically thin cooling is too slow in thiscase, optically thick cooling — above the floor set byoptically thin cooling — is also too slow.Larger opacities make optically thin cooling muchfaster than β crit , as shown in the top panels of Fig. 18.However with larger opacities, optically thick cooling af-fects larger disk radii, slowing the cooling. Remarkably,the adopted MMSN value of opacity hits the a sweet spotwhere optically thin cooling is fast enough, yet slower op-tically thick can be restricted to inner disk radii.At smaller disk radii, Fig. 18 shows the hallmark of dif-fusive cooling, that larger wavenumbers can cool faster,but not below the floor set by optically thin cooling.Optically thick cooling times also rise sharply towardsmaller radii (as β ∝ r − / ) due to high densities andshort orbital times. This effect suppresses VSI growth atsmall radii, but with a strong wavenumber dependence.Fig. 19 highlights this wavenumber dependence byplotting the wavenumbers for which β = β crit . Abovethe solid curves (i.e. for larger wavenumbers), β < β crit at all disk heights, implying linear VSI growth. Below thedashed curves, VSI growth is strongly suppressed since β > β crit for all | z | < H . Between the solid and dashedcurves some growth is possible, but only fairly close tothe solid curve (according to § k x & /H . Asargued elsewhere, such small-scale disturbances may notdrive significant turbulence or transport.We thus doubt that the VSI is significant at 1 AU inMMSN-like PPDs, for any opacity. At higher opacities,the required wavenumbers become shorter and even moreproblematic, as shown by the red curves in Fig. 19. Atlower opacities, even optically thin cooling (the fastestpossible) is too slow, as discussed above. VSI growth in the MMSN
We now consider the growth of the VSI in a MMSNdisk model. This numerical calculation is similar to § -3 -2 -1 β = t c Ω K × MMSN opacity, z=0 -3-2-1 × MMSN opacity, z=3H -3 -2 -1 β = t c Ω K β crit β , k x H=1 β , k x H=10 β , k x H=100MMSN opacity, z=0 -3-2-1 r/AU10 -3 -2 -1 β = t c Ω K × MMSN opacity, z=0 r/AU -3-2-1 × MMSN opacity, z=3H
Figure 18.
Dimensionless thermal relaxation timescales β , evaluated at the midplane (left) and at z = 3 H (right) in the fiducial PPD.Eq. 59 is plotted for three values of the perturbation radial wavenumber: ˆ k = 1 (dotted), ˆ k = 10 (dashed) and ˆ k = 100 (dashed-dot), forthree values of the opacity relative to the MMSN: ˆ κ d = 0 . κ d = 1 (middle) and ˆ κ d = 10 (bottom). The solid line is the criticalthermal relaxation timescale β crit . r/AU10 -4 -2 m i n ( k x H ) f o r β < β c r i t MMSN opacity, z=0MMSN opacity, z=3H10 × MMSN opacity, z=010 × MMSN opacity, z=3H
Figure 19.
The minimum perturbation wavenumber ˆ k in theMMSN such that the associated dimensionless thermal relaxationtime β at z = 0 (solid) and z = 3 H (dashed) is less than the criti-cal timescale β crit . Asterisks correspond to the inner cut-off radiishown in Fig. 20. mode for a range of wavenumbers and disk radii. Wefocus on the fundamental, i.e. lowest order vertical, modebecause it is the fastest growing mode except for somesurface modes at high wavenumbers. We neglect thesesurface modes for reasons discussed in § § r AU ∼ r AU ∼
50 with growth timescales ∼ r increases towards100AU in Fig. 20. This trend is expected as the outerradius cutoff at ∼ CAVEATS AND NEGLECTED EFFECTS
Viscosity
Our neglect of viscosity is valid in a laminar disk, be-cause molecular viscosity is so small. However, turbu-lence could act as an effective viscosity which preferen-tially damps small scale modes. Since the goal of VSI isto drive turbulent transport, the most relevant modes forsustaining the VSI should be able to overcome turbulentertical shear instability 15 t g r o w / P o r b k x H=1k x H=5k x H=10k x H=20k x H=30k x H=50
Figure 20.
VSI growth times ( t grow = 1 /σ ) in orbital units( P orb = 2 π/ Ω K ) for the reference MMSN disk model with self-consistent dust cooling. damping.To roughly estimate the wavelengths which are proneto damping, we consider the standard prescription forthe kinematic viscosity ν = αc s H (Shakura & Sunyaev1973), in terms of the dimensionless α parameter. Theviscous timescale for a perturbation lengthscale l ∼ /k x = H/ ˆ k is t visc = 1 /α ˆ k Ω. For significant growth, t visc should be longer than the characteristic VSI growthtimescale 1 / ( h | q | Ω), i.e. ˆ k . | q | h/α . With h ∼ . q = − α from simulations (N13), we estimate thatgrowth requires ˆ k .
20 for α ∼ − or ˆ k .
70 for α ∼ − .This argument justifies our focus on moderatewavenumbers, ˆ k = O (10). Future work should considerthe viscous effects in more detail, to better understandhow the VSI operates in nature and in simulations. Convective overstability
The VSI bears some similarity to the ‘convective over-stability’ (CONO, Klahr & Hubbard 2014; Lyra 2014).Both instabilities rely on thermal relaxation to overcomethe stability of disks to adiabatic perturbations, i.e. theSolberg-Hoiland criteria (see § N r ≡ − ρC P ∂P∂r ∂S∂r < . (62)Since vertical buoyancy is an important stabilizing influ-ence for the VSI, it should be considered in future studiesof CONO.The sign of N r depends on disk parameters. Todemonstrate that the parameters needed for N r < ∝ r n so that n = p + q/ /
2. Using Eqs. 22a and 62 with γ = 1 . N r < n > / q ). For stan-dard disk temperature laws, this requirement implies aflat or rising surface density profile, e.g. n = 3 /
14 for q = − / n = 0 for q = − /
2. Since standard Σ pro-files decline with radius, CONO is most like to operateat special locations, like the outside edges of disk gaps orholes and/or shadowed regions where q < − / k . β ∼ Radiative transfer
Our relatively simple treatment of cooling with an ide-alized dust opacity could certainly be generalized in fu-ture works. In hotter disks, gas phase opacities must beconsidered. In cold disks, there are many choices for thedust opacity, which varies with grain sizes and composi-tions. Changing dust properties would alter the viabilityof the VSI for better or worse. The radiative transferitself could be calculated with higher levels of sophisti-cation, as is already being done in numerical simulationsof the VSI (Stoll & Kley 2014). SUMMARY AND DISCUSSION
In this paper we study the vertical shear instability(VSI) with a focus the role of radiative cooling. In turnwe assess the viability of the VSI as an angular mo-mentum transport mechanism in protoplanetary disks(PPDs).Our linear, axisymmetric analysis of the VSI considers(uniquely to our knowledge) finite cooling times in a ver-tically global model. In order for vertical shear to drivethe VSI, short cooling times are needed to weaken thestabilizing effects of vertical buoyancy. Our main analyt-ical finding, which we confirm numerically, is the criticalcooling timescale above which VSI growth is suppressed,Eqs. 5 and 46.Our main focus is irradiated, vertically isothermaldisks which have strong vertical buoyancy. The criti-cal cooling time is thus short, shorter than the orbitaltime by a factor of the disk aspect ratio. This findingis consistent with, and helps explain, the results of re-cent numerical simulations (N13). We briefly considervertically non-isothermal disks in § . ∼ < APPENDIX A. DERIVATION OF THE APPROXIMATE EQUATIONS
Here we detail the derivation of Eq. 29 used in the analytical discussion of §
5. Starting with Eqs. 27a — 27e, weset Γ = 1 and ζ = 0 for the vertically isothermal limit and the fully-radially-local approximation, respectively. Weeliminate the horizontal velocity perturbations ( δv x , δv y ) to obtaini υδv z = dWdz + ∂ ln ρ∂z ( W − Q ) , (A1a) υ c s Q + υ k x D W = i υ (cid:18) i k x rD ∂ Ω ∂z − ∂ ln ρ∂z (cid:19) δv z − i υ dδv z dz , (A1b) υ W − γυ Q + i υt c ( W − Q ) = i υc s ∂ ln ρ∂z ( γ − δv z , (A1c)where D ≡ κ − υ . (A2)Reduction to a single ODE requires ∂ z D = ∂ z κ . At this point we could apply the low frequency and Keplerianapproximations to set D → κ → Ω , then D is vertically constant, and we can obtain Eq. 29 more directly. However,to demonstrate that the order of approximation is irrelevant, we will retain ∂ z κ initially. Using Eq. 19—21, ∂κ ∂z = 4 ∂ Ω ∂z + r ∂∂r ∂ Ω ∂z = − d ln ρdz qc s r (cid:18) z r + z − (cid:19) ≡ − d ln ρdz qc s r F ( r, z ) . (A3)The function F increases monotonically from F = − z = 0 to F → | z | → ∞ , so | F | = O (1).We eliminate W, Q from Eqs. A1, using Eqs. 31 and A3 to obtain0 = d δv z dz + (cid:20) k x c s qDr − k x c s ( k x c s + χD ) qc s FDr (cid:21) d ln ρdz dδv z dz + ( υ (cid:18) k x D + χc s (cid:19) + (cid:18) χ + i k x c s qDr (cid:19) d ln ρdz − c s D (cid:18) d ln ρdz (cid:19) (cid:18) k x − i k x qr (cid:19) (cid:20) (1 − χ ) + χ ( k x c s + χD ) qc s Fr (cid:21)) δv z , (A4)where we have replaced ∂/∂z by d/dz for the fully-radially-local treatment. Now we make the low frequency andKeplerian approximations, setting D → Ω , to give δv ′′ z + hq ˆ k − ˆ k (cid:16) ˆ k + χ (cid:17) qh F ln ρ ′ δv ′ z + (cid:16) χ + i hq ˆ k (cid:17) ln ρ ′′ − ln ρ ′ (cid:16) ˆ k − i hq ˆ k (cid:17) − χ + χ (cid:16) ˆ k + χ (cid:17) qh F δv z = − ˆ υ (cid:16) ˆ k + χ (cid:17) δv z , (A5)in terms of dimensionless variables of Eq. 28. Retaining terms to first order in the disk aspect ratio, h ≪
1, gives0 = δv ′′ z + (cid:16) hq ˆ k (cid:17) ln ρ ′ δv ′ z + h ˆ υ (cid:16) ˆ k + χ (cid:17) + (cid:16) χ + i hq ˆ k (cid:17) ln ρ ′′ − ln ρ ′ (cid:16) ˆ k − i hq ˆ k (cid:17) (1 − χ ) i δv z . (A6)Approximating the density field by Eq. 17b then gives Eq. 29.We can establish a correspondence between our Eq. 29 and Eq. 41 in Lubow & Pringle (1993), which describesadiabatic axisymmetric waves in a vertically isothermal disk without vertical shear. Accounting for the required changeof variables, the correspondence is exact after setting q = 0 (no vertical shear) and χ = 1 /γ (adiabatic flow) in ourEq. 29, and applying the approximations in our § B. APPLICABILITY OF THE FULLY-RADIALLY-LOCAL APPROXIMATION
In the fully-radially-local approximation, background radial gradients are ignored except where it appears implicitlyin the expression for the vertical shear rate ∂ z Ω (via Eq. 21). This is done by neglecting the terms proportional to ζ in Eq. 27a, 27b and 27e, i.e. setting ζ = 0. (Nominally ζ = 1, which is used in our numerical study.)For a power-law disk, the neglected radial gradients are O ( r − ), and they appear in comparison with terms of O ( k x ).The neglected terms therefore have a relative magnitude of O ( h/ ˆ k ), which is small for thin disks ( h ≪
1) and/or smallradial wavelengths (ˆ k ≫ B.1.
Spurious growth of adiabatic perturbations when ζ = 0A limitation of the reduced model described in Appendix A, § §
5, is that it cannot be employed foradiabatic flow when there is vertical shear, even if h/ ˆ k ≪
1. We explain this by setting β → ∞ and hence χ = 1 /γ inEq. A6 to give0 = δv ′′ z + (cid:16) hq ˆ k (cid:17) (ln ρ ′ δv z ) ′ + (cid:26) ˆ υ (cid:18) γ + ˆ k (cid:19) − (cid:18) γ − γ (cid:19) (cid:20) ln ρ ′′ + ˆ k (cid:18) − i hq ˆ k (cid:19) ln ρ ′ (cid:21)(cid:27) δv z . (B1)We multiply Eq. B1 by ρδv ∗ z and integrate vertically, assume boundary terms vanish when integrating by parts, toobtainˆ υ (cid:18) γ + ˆ k (cid:19) Z ˆ z ˆ z ρ | δv z | d ˆ z = (cid:18) γ − γ (cid:19) Z ˆ z ˆ z ρ | δv ′ z | d ˆ z + 1 γ Z ˆ z ˆ z ρ | ( ρδv z ) ′ | d ˆ z + (cid:18) γ − γ (cid:19) ˆ k (cid:18) − i hq ˆ k (cid:19) Z ˆ z ˆ z ρ ln ρ ′ | δv z | d ˆ z + i hq ˆ k Z ˆ z ˆ z ln ρ ′ ( ρδv ∗ z ) ′ δv z d ˆ z. (B2)In the presence of vertical shear q = 0, Eq. B2 shows that ˆ υ is complex for real ˆ k , which indicates instability forany value of γ >
1. This stability condition is inconsistent with the second Solberg-Hoiland criterion (Eq. 25), whichstates that instability requires the disk to be close to neutral stratification (i.e. | γ − | ≪ ζ = 0 model arises because we have retained the global radial temperature gradientto obtain vertical shear, but have ignored it elsewhere in the linear problem (as well as the background radial densitygradient). Nevertheless, we demonstrate below that this inconsistency is unimportant for the VSI, which occurs for β ≪
1, not for adiabatic perturbations.
B.2.
Effect of global radial gradients
In Fig. 21, we plot the effect of the global radial gradient terms proportional to ζ in Eq. 27a—27e by calculatingthe fundamental VSI growth rates using three approaches. We compute growth rates from the dispersion relation Eq.44, which assumes ζ = 0; from Eq. 27a—27e with ζ = 0; and from Eq. 27a—27e with ζ = 1.All three methods give similar behavior, and growth rates are in close agreement for β .
1. Differences arise for β &
1, and as β → ∞ the fully-radially-local approximation, where ζ = 0, gives a (spurious) growth rate as expectedfrom the discussion above. Inclusion of the global radial gradient terms results in the expected behavior ( σ → β → ∞ ). Despite this caveat, Fig. 21 shows that provided we consider β . O (1), then setting ζ = 0 does not affectthe VSI significantly. C. LINEAR PROBLEM IN THE ISOTHERMAL LIMIT
Here we summarize selected results for isothermal perturbations ( β ≡
0) in vertically isothermal disks (Γ = 1) inthe fully-radially local approximation ( ζ = 0). In this case Eq. A1c becomes Q = W (C1)(i.e. δP = c s δρ ). For isothermal perturbations it is simpler to work with an equation for W by substituting Q = W into Eq. A1b and using Eq. A1a to eliminate δv z . In this case, we shall not yet make the low frequency approximation,but first make the Keplerian approximation. We obtain, in terms of dimensionless variables,0 = W ′′ + " ln ρ ′ − i hq ˆ kf (ˆ z )1 − ˆ υ W ′ + ˆ υ k − ˆ υ ! W, (C2)where for discussion purposes we have defined f (ˆ z ) such that d Ω d ˆ z = h qf (ˆ z )Ω . (C3)8 M-K. Lin, A. N. Youdin -3 -2 -1 β -0.6-0.4-0.20.00.20.4 σ / ( h Ω K ) Analytical ( ζ =0)Numerical ( ζ =0)Numerical ( ζ =1) Figure 21.
Growth rate of the fundamental VSI mode as a function of the thermal relaxation time β . The ‘Analytical ( ζ = 0)’ curveis calculated from the dispersion relation Eq. 44. The ‘Numerical ( ζ = 0)’ curve is obtained from the numerical solution to Eq. 27a—27e with ζ = 0, which neglects explicit dependencies on the radial disk structure. The ‘Numerical ( ζ = 1)’ curve is obtained from Eq.27a—27e with ζ = 1, which accounts for the radial disk structure self-consistently. The disk parameters are ( γ, Γ) = (1 . , . p, q, h ) = ( − . , − , . k = 30. By comparing Eq. C3 with Eq. 20, we see that f (ˆ z ) = ˆ z (cid:0) h ˆ z (cid:1) − / . More generally, though, f may be regarded asa representation of the vertical shear profile. Physically, we expect there is a maximum value of | d Ω /dz | , the existenceof which should limit the growth rate of the VSI, as remarked in § C.1.
Maximum growth rate in the low-frequency limit
Here we consider the low-frequency limit | ˆ υ | ≪ W ′′ + h ln ρ ′ − i hq ˆ kf (ˆ z ) i W ′ + ˆ υ (cid:16) k (cid:17) W. (C4)We multiply Eq. C4 by ρW ∗ and integrate vertically from ˆ z = ˆ z to z = ˆ z . We neglect boundary terms whenintegrating by parts, by assuming W or W ′ vanishes at the boundaries, or that the boundaries are sufficiently far awayso that the boundary terms are negligible because of the decaying background density with increasing height. Then,ˆ υ (cid:16) k (cid:17) Z ˆ z ˆ z ρ | W | d ˆ z = Z ˆ z ˆ z ρ | W ′ | d ˆ z + i hq ˆ k Z ˆ z ˆ z ρf (ˆ z ) W ∗ W ′ d ˆ z. (C5)It follows that for instability (Im ˆ υ > q = 0 or more generally d Ω /dz = 0, i.e. there must bevertical shear.The real and imaginary parts of Eq. C5 are (cid:0) ˆ ω − ˆ σ (cid:1) (cid:16) k (cid:17) Z ˆ z ˆ z ρ | W | d ˆ z − Z ˆ z ˆ z ρ | W ′ | d ˆ z = Re " i h ˆ kq Z ˆ z ˆ z ρf (ˆ z ) W ∗ W ′ d ˆ z , ω ˆ σ (cid:16) k (cid:17) Z ˆ z ˆ z ρ | W | d ˆ z = Im " i h ˆ kq Z ˆ z ˆ z ρf (ˆ z ) W ∗ W ′ d ˆ z , (C6)ertical shear instability 19where we recall ˆ ω = Re ˆ υ and ˆ σ = Im ˆ υ . Adding the square of these equations give " | ˆ υ | (cid:16) k (cid:17) Z ˆ z ˆ z ρ | W | d ˆ z − Z ˆ z ˆ z ρ | W ′ | d ˆ z + 4ˆ σ (cid:16) k (cid:17) Z ˆ z ˆ z ρ | W | d ˆ z Z ˆ z ˆ z ρ | W ′ | d ˆ z = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) i h ˆ kq Z ˆ z ˆ z ρf (ˆ z ) W ∗ W ′ d ˆ z (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) . (C7)It is clear that 4ˆ σ (cid:16) k (cid:17) Z ˆ z ˆ z ρ | W | d ˆ z Z ˆ z ˆ z ρ | W ′ | d ˆ z ≤ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) h ˆ kq Z ˆ z ˆ z ρf (ˆ z ) W ∗ W ′ d ˆ z (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) . (C8)On the left hand side of this inequality, we apply the Cauchy-Schwarz inequality to obtain Z ˆ z ˆ z ρ | W | | W ′ | d ˆ z ! ≤ Z ˆ z ˆ z ρ | W | d ˆ z Z ˆ z ˆ z ρ | W ′ | d ˆ z. (C9)On the right hand side of Eq. C8 we have (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)Z ˆ z ˆ z ρf (ˆ z ) W ∗ W ′ d ˆ z (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ≤ Z ˆ z ˆ z ρ | f (ˆ z ) W ∗ W ′ | d ˆ z ≤ max ( | f | ) Z ˆ z ˆ z ρ | W | | W ′ | d ˆ z, (C10)where max( | f | ) is the maximum value of | f | in ˆ z ∈ [ˆ z , ˆ z ]. Inserting these inequalities into Eq. C8 gives | ˆ σ | ≤ h | ˆ kq | p k max( | f | ) < h | q | | f | ) . (C11)It follows that the maximum possible growth rate of unstable modes, satisfying the above boundary conditions, islimited by the maximum vertical shear rate in the domain considered, | σ | < max (cid:12)(cid:12)(cid:12)(cid:12) r d Ω dz (cid:12)(cid:12)(cid:12)(cid:12) , (C12)as expected on physical grounds. Note that if the thin-disk approximation is used in an infinite domain, then themaximum growth rate is unbounded (since in that case max | f | → ∞ ). However, large growth rates invalidate thelow-frequency approximation and the above analysis breaks down.In practice, one might consider a vertical domain of a few scale heights in a thin disk with | q | = O (1). Then f ≃ ˆ z ,so that max | f | = O (1), implying a maximum growth rate O ( h Ω K ), consistent with numerical results. C.2.
Explicit solutions in the thin-disk limit
Here we assume a thin disk ( h ≪
1) so that ln ρ ≃ − ˆ z / f (ˆ z ) ≃ ˆ z . However, we do not assume low frequencyfrom the onset. Then Eq. C2 becomes W ′′ − qh ˆ k − ˆ υ ! ˆ zW ′ + ˆ υ k − ˆ υ ! W = 0 . (C13)We remark that taking the low frequency limit and considering ˆ k ≫
1, Eq. C13 becomes equivalent to Eq. 39 in N13or Eq. 28 in BL15, although we have taken a different route.We seek power series solutions to Eq. C13, W (ˆ z ) = ∞ X l =0 a l ˆ z l . (C14)Then the coefficients a l are given by the recurrence relation( l + 2)( l + 1) a l +2 + " ˆ υ k − ˆ υ ! − l hq ˆ k − ˆ υ ! a l = 0 . (C15)We demand the series to terminate at l = L , i.e. a polynomial of order L , so that the vertical kinetic energy densityremain bounded as | ˆ z | → ∞ . Then the eigenfrequency ˆ υ is given viaˆ υ − (cid:16) L + 1 + ˆ k (cid:17) ˆ υ + L (cid:16) hq ˆ k (cid:17) = 0 . (C16)Note that we have applied a regularity condition at infinity, since the vertically isothermal disk has no surface. Ifvertical boundaries are imposed at finite height, as done in numerical calculations, then the above solution needs to be0 M-K. Lin, A. N. Youdinmodified to match the desired boundary conditions. This enables the ‘surface modes’ seen in numerical calculations(BL15). C.2.1.
Stability without vertical shear
In the absence of vertical shear q = 0, Eq. C16 can be written asˆ υ ˆ k = (cid:0) L − ˆ υ (cid:1) (cid:0) − ˆ υ (cid:1) , (C17)which is just the dispersion relation for axisymmetric isothermal waves in a vertically isothermal disk (e.g.Okazaki et al. 1987; Takeuchi & Miyama 1998; Tanaka et al. 2002; Zhang & Lai 2006; Ogilvie & Latter 2013;Barker & Ogilvie 2014; BL15). In this case the solutions are Hermite polynomials, W ∝ He L (ˆ z ). The eigenfre-quency ˆ υ = ˆ ω is real and the disk is stable. The low frequency branch of Eq. C17 are inertial waves (Balbus 2003).For | ˆ ω | ≪ L ≥ ω ˆ k = L , or ˆ ω ∝ ˆ k − for fixed L . This inverse relation has beenqualitatively observed in numerical simulations of Stoll & Kley (2014). Note that L = M +1, where M ≥ §
5) based on the reduced equation for δv z , rather than for W as consideredhere. C.2.2.
Instability with vertical shear
The VSI corresponds to unstable inertial waves. This is readily obtained for large ˆ k by balancing the last two termsof Eq. C16 to give the low frequency branch. Thenˆ υ ≃ L qh ˆ kL + 1 + ˆ k ! , (C18)which is equivalent to Eq. 34 of BL15 in the limit ˆ k ≫ L,
1. This signifies instability for L ≥ υ >
0. These are the low frequency unstable modes seen in Fig. 4 and Fig. 7,for which σ ∝ | ω | . D. ANALYTIC DISPERSION RELATION WITH THERMAL RELAXATION
D.1.
Coefficients
The coefficients of the dispersion relation Eq. 44 are: c = M ( M + 1) A , (D1a) c = i β n (1 − γ ) h k (1 + 2 M ) − hq ˆ kM ( M + 1) i − A γM ( M + 1) o , (D1b) c = (cid:16) ˆ k + 1 (cid:17) A + β n (1 − γ ) h γ ˆ k (1 + 2 M ) − hq ˆ kγM ( M + 1) i − γ A M ( M + 1) o , (D1c) c = β n hq ˆ k + γ h i + hq ˆ k (cid:16) k (cid:17)i − − k o , (D1d) c = β (cid:16) γ ˆ k (cid:17) h γ (cid:16) − i hq ˆ k (cid:17) − i − (cid:16) k (cid:17) , (D1e) c = 2i β (cid:16) k (cid:17) (cid:16) γ ˆ k (cid:17) , (D1f) c = β (cid:16) γ ˆ k (cid:17) . (D1g) D.2.
Finding marginal stability
To investigate marginal stability we set ˆ σ = 0 so the frequency, ˆ υ = ˆ ω c , is real; and β = β c , the cooling time formarginal stability. We take the short radial wavelength limit, ˆ k ≫
1, of the coefficients in Eq. D1. We consider M . O (1) and assume β c ≪ ˆ k . The real and imaginary parts of Eq. 44 then give relations for ˆ ω c and β c as0 = M ( M + 1)(1 − h q ˆ k ) + 4 β c hqM ( M + 1) (cid:16) ˆ ω c ˆ k (cid:17) + (cid:2) γβ c (1 − γ ) (1 + 2 M ) (cid:3) (cid:16) ˆ ω c ˆ k (cid:17) (D2)+ 2 hqγβ c (cid:16) ˆ ω c ˆ k (cid:17) − (cid:16) ˆ ω c ˆ k (cid:17) + β c γ ˆ ω c (cid:16) ˆ ω c ˆ k (cid:17) , hqM ( M + 1)ˆ k + β c (1 − γ )(1 + 2 M ) ˆ k (cid:16) ˆ ω c ˆ k (cid:17) + hq ˆ k (cid:16) ˆ ω c ˆ k (cid:17) − β c ˆ ω c (cid:16) ˆ ω c ˆ k (cid:17) − hqγ β c ˆ ω c (cid:16) ˆ ω c ˆ k (cid:17) (D3)+ 2 β c γ ˆ ω c (cid:16) ˆ ω c ˆ k (cid:17) . ertical shear instability 21We recall that in the low frequency and thin-disk approximations, h, | ˆ ω c | ≪
1. We note that | q | = O (1) and γ > O (1). We assume | ˆ ω c ˆ k | = O (1), since for inertial waves | ˆ ω ˆ k | ≃ √ M . Finally, we further assume β c ≪
1, to bejustified a posteriori , to give Eq. 45a—45b.
D.3.
Maximum critical cooling time
Here we show that for sufficiently thin disks the longest critical cooling time is that for the M = 0 or fundamentalmode. This allows us to focus on the fundamental mode to obtain an overall cooling requirement for the VSI.Consider the simplified dispersion relations for marginal stability, Eq. 45a—45b. We write X = ˆ ω ˆ k , θ = ( hq ˆ k ) andtreat M as a continuous variable. We find from Eq. 45a that2 X dXdM = X (1 − θ )(2 M + 1) X + 2 M ( M + 1)(1 − θ ) , (D4)and from Eq. 45b that2 X dXdM = (cid:2) X + 2 M ( M + 1) (cid:3) (cid:18) d ln β c dM + 42 M + 1 + d ln XdM (cid:19) − M + 1) . (D5)We eliminate dX/dM , making use of Eq. 45a in the process, to obtain(2 M + 1) d ln β c dM = − (1 − θ ) M ( M + 1) (cid:2) M + 1) + 12 X (cid:3) + (3 + θ ) X X + 2 M ( M + 1)(1 − θ )] [ X + 2 M ( M + 1)] . (D6)Hence, dβ c dM < , (D7)for all M ≥ θ ≤
1, which imply max( β c ) occurs at M = 0. This conclusion may also be reached by explicitlysolving Eq. 45a—45b with θ as a small parameter. For fixed ˆ k the condition θ ≤ h . All such modes are stabilized if β > β c ( M = 0) ≡ β crit .This result highlights the importance of the fundamental mode — it is the most difficult mode to stabilize with finitecooling. Furthermore, for M = 0 we may obtain the expression for β crit from the dispersion relations Eq. D2—D3without assuming it is much less than unity at the outset or place restrictions on θ . REFERENCESArmitage, P. J. 2011, ARA&A, 49, 195Bai, X.-N. 2015, ApJ, 798, 84Balbus, S. A. 2003, ARA&A, 41, 555Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214Balbus, S. A., Hawley, J. F., & Stone, J. M. 1996, ApJ, 467, 76Barker, A. J., & Latter, H. N. 2015, ArXiv e-prints, arXiv:1503.06953Barker, A. J., & Ogilvie, G. I. 2014, MNRAS, 445, 2637Bell, K. R., & Lin, D. N. C. 1994, ApJ, 427, 987Bitsch, B., Johansen, A., Lambrechts, M., & Morbidelli, A. 2015, A&A, 575, A28Blaes, O. M., & Balbus, S. A. 1994, ApJ, 421, 163Chandrasekhar, S. 1961, Hydrodynamic and hydromagnetic stabilityChiang, E., & Youdin, A. N. 2010, Annual Review of Earth and Planetary Sciences, 38, 493Chiang, E. I., & Goldreich, P. 1997, ApJ, 490, 368Fricke, K. 1968, ZAp, 68, 317Gammie, C. F. 1996, ApJ, 457, 355Goldreich, P., & Lynden-Bell, D. 1965, MNRAS, 130, 125Goldreich, P., & Schubert, G. 1967, ApJ, 150, 571Gressel, O., Turner, N. J., Nelson, R. P., & McNally, C. P. 2015, ArXiv e-prints, arXiv:1501.05431Klahr, H., & Hubbard, A. 2014, ApJ, 788, 21Lee, A. T., Chiang, E., Asay-Davis, X., & Barranco, J. 2010, ApJ, 718, 1367Lesniak, M. V., & Desch, S. J. 2011, ApJ, 740, 118Lesur, G., Kunz, M. W., & Fromang, S. 2014, A&A, 566, A56Lubow, S. H., & Pringle, J. E. 1993, ApJ, 409, 360Lyra, W. 2014, ApJ, 789, 77McNally, C. P., & Pessah, M. E. 2014, ArXiv e-prints, arXiv:1406.4864Nelson, R. P., Gressel, O., & Umurhan, O. M. 2013, MNRAS, 435, 2610Ogilvie, G. I., & Latter, H. N. 2013, MNRAS, 433, 2420Okazaki, A. T., Kato, S., & Fukue, J. 1987, PASJ, 39, 457Salmeron, R., & Wardle, M. 2003, MNRAS, 345, 992Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337Simon, J. B., Bai, X.-N., Stone, J. M., Armitage, P. J., & Beckwith, K. 2013, ApJ, 764, 66Stoll, M. H. R., & Kley, W. 2014, A&A, 572, A77Takeuchi, T., & Miyama, S. M. 1998, PASJ, 50, 141Tanaka, H., Takeuchi, T., & Ward, W. R. 2002, ApJ, 565, 12572 M-K. Lin, A. N. Youdin