On the Vertical Shear Instability in Magnetized Protoplanetary Disks

The vertical shear instability (VSI) is a robust phenomenon in irradiated protoplanetary disks (PPDs). While there is extensive literature on the VSI in the hydrodynamic limit, PPDs are expected to be magnetized and their extremely low ionization fractions imply that non-ideal magneto-hydrodynamic (MHD) effects should be properly considered. To this end, we present linear analyses of the VSI in magnetized disks with Ohmic resistivity. We primarily consider toroidal magnetic fields, which are likely to dominate the field geometry in PPDs. We perform vertically global and radially local analyses to capture characteristic VSI modes with extended vertical structures. To focus on the effect of magnetism, we use a locally isothermal equation of state. We find that magnetism provides a stabilizing effect to dampen the VSI, with surface modes, rather than body modes, being the first to vanish with increasing magnetization. Subdued VSI modes can be revived by Ohmic resistivity, where sufficient magnetic diffusion overcome magnetic stabilization, and hydrodynamic results are recovered. We also briefly consider poloidal fields to account for the magnetorotational instability (MRI), which may develop towards surface layers in the outer parts of PPDs. The MRI grows efficiently at small radial wavenumbers, in contrast to the VSI. When resistivity is considered, we find the VSI dominates over the MRI for Ohmic Els\"{a}sser numbers $\lesssim 0.09$ at plasma beta parameter $\beta_Z \sim 10^4$.


INTRODUCTION
It has been postulated over decades that the turbulence and angular momentum transport in most astrophysical accretion disks are mediated by the magnetorotational instability (MRI; Balbus & Hawley 1991). However, protoplanetary disks (PPDs) are distinguished by their extremely weakly ionized gas (Gammie 1996;Armitage 2011), where gas and magnetic fields are poorly coupled, and the MRI turbulence is either quenched or dampened in the bulk of the disk (Perez-Becker & Chiang 2011;Bai & Stone 2011;Simon et al. 2013a,b). Instead, angular momentum transport is dominated by magnetized disk winds, leaving the main disk mostly laminar (Bai & Stone 2013;Gressel et al. 2015;Bai et al. 2016;Bai 2017;Béthune et al. 2017;Gressel et al. 2020).
Nevertheless, some level of turbulence is expected in PPDs to account for the recent ALMA observations of molecular line emissions (Teague et al. 2016;Flaherty et al. 2017Flaherty et al. , 2018Flaherty et al. , 2020. Furthermore, turbulence may serve as an essential ingredient in many stages of planet formation. Turbu-⋆ E-mail: cc795@cam.ac.uk † E-mail: mklin@asiaa.sinica.edu.tw lence affects the gravitational sedimentation (Dubrulle et al. 1995;Johansen & Klahr 2005;Youdin & Lithwick 2007), radial diffusion (Clarke & Pringle 1988), and collisional growth (Ormel & Cuzzi 2007;Birnstiel et al. 2010) of dust particles. Long-lived vortices induced by turbulence (e.g. Raettig et al. 2015;Manger & Klahr 2018) can concentrate dust particles (Barge & Sommeria 1995;Klahr & Henning 1997;Cuzzi et al. 2008) and seed planetesimal formation through streaming instability or gravitational instability (Youdin & Goodman 2005a;Johansen et al. 2007;Chiang & Youdin 2010), whereas the growth of streaming instability can be substantially diminished by a moderate level of turbulent viscosity (Chen & Lin 2020;Umurhan et al. 2020). Turbulence also influences the radial migration of planets, as well as the flow morphology and gap formation around them Papaloizou et al. 2004). Therefore, understanding the origin and characteristics of turbulence in PPDs is essential to many aspects of planet formation and evolution.
The VSI is inherited from the Goldreich-Schubert-Fricke instability (Goldreich & Schubert 1967;Fricke 1968) and is initially discovered in the context of differentially rotating stars. Its importance to accretion disks was later explored by Urpin & Brandenburg (1998), Urpin (2003) and Arlt & Urpin (2004). The applicability of the VSI to PPDs has been demonstrated only recently in Nelson et al. 2013. It quickly drew intensive interests (e.g. Stoll & Kley 2014;Barker & Latter 2015;Umurhan et al. 2016a;Latter & Papaloizou 2018;Lin 2019;Cui & Bai 2020;Schäfer et al. 2020) and is considered to be a promising hydrodynamic mechanism in driving turbulence in PPDs.
A differentially rotating disk with Keplerian profile is stable according to the Rayleigh criterion (Chandrasekhar 1961). The presence of vertical shear can destabilize inertial waves in a vertically global disk model (Barker & Latter 2015). Nevertheless, a fluid element also experiences stabilizing effects from vertical buoyancy, which impedes the VSI growth. This can be overcome by sufficiently rapid cooling that brings the perturbed fluid element to reach local thermal equilibrium with its surroundings, hence diminishing the buoyancy. Local linear analyses demonstrate that the unstable modes are characterized by short radial wavelengths and maximum growth rates much smaller than the orbital frequency (Urpin 2003;N13). Two classes of VSI modes have been identified: rapidly-growing surface modes concentrated near the disk surface and more vertically extended body modes (N13; Barker & Latter 2015;Lin & Youdin 2015). The body modes can be further categorized into breathing and corrugation modes, depending on the symmetry about the midplane (N13).
The non-linear evolution of the VSI has been examined by hydrodynamic simulations. In accordance with the linear theory, the VSI is triggered when thermal relaxation timescales are less than 0.01 − 0.1 times the local dynamical timescales, and the wave modes exhibit elongated vertical wavelengths (N13). While the surface modes possess the fastest growth rate, the body (corrugation) modes eventually take over, dominating the non-linear evolution. Fully developed VSI turbulence yield a Shakura-Sunyaev (Shakura & Sunyaev 1973) α value on the order of 10 −4 − 10 −3 (N13, Stoll & Kley 2014). Non-axisymmetric 3D simulations show the development of vortices (Richard et al. 2016;Manger & Klahr 2018;Flock et al. 2020;Manger et al. 2020;Pfeil & Klahr 2020). Incorporating dust particles, numerical simulations show that VSI can stir up dust grains against vertical settling, but may also concentrate them through inducing dusttrapping vortices (Stoll & Kley 2016;Flock et al. 2017;Lin 2019;Schäfer et al. 2020).
Most studies of the VSI to date have neglected magnetic fields despite their importance in the evolution of PPDs. Two recent works extend analyses of the VSI to MHD regimes. Local linear stability analyses in the ideal MHD limit show that weak magnetic fields are favoured to excite the VSI, specifically when plasma beta β 400 for thin disks (LP18), where β is the ratio of gas to magnetic pressure. They also find the MRI growth rates exceed that for VSI modes, but the wave vectors of the two are perpendicular to each other. For resistive disks, LP18 estimates a critical Ohmic Elsässer number of ∼ 1/h, where h is the disk aspect ratio, for the VSI to operate, by requiring the magnetic diffusion timescale to be shorter than the Alfvén wave propagation timescales. Non-linear MHD simulations, applicable to outer regions of the disk, demonstrate that the VSI can initiate and sustain turbulence in magnetized disks (Cui & Bai 2020). Weak ambipolar diffusion strength, or the enhanced coupling between gas and magnetic fields, works as stabilizing effects to dampen the VSI growth.
Previous hydrodynamic models suggest that effective VSI growths span over ∼ 5 − 100 AU in PPDs (Lin & Youdin 2015;Pfeil & Klahr 2019). These regions are susceptible to all three non-ideal MHD effects -Ohmic resistivity, Hall effect, and ambipolar diffusion (Wardle 2007;Bai 2011). However, a quantitative analysis in vertically-global, magnetized disks with non-ideal effects is still lacking. Such analyses can be useful for understanding the VSI mode properties in real PPDs and interpreting nonlinear simulations. Hence, in this work, we extend the linear stability analysis of the VSI to weakly ionized gas in a vertically global disk model. We remark that a vertically-global analysis is necessary for a proper description of elongated body modes of the VSI, which have been found to dominate in numerical simulations (N13; Stoll & Kley 2014;Cui & Bai 2020). We consider ideal MHD and further include Ohmic resistivity as a proxy for non-ideal effects. We focus primarily on the effect of toroidal magnetic fields. However, in a local model, we also investigate the dominance between MRI and VSI by considering purely poloidal magnetic fields.
The plan of the paper is as follows: in §2, we introduce the basic formulation and establish the equilibrium state of the problem. In §3, we derive and discuss the Solberg-Hoiland stability criteria for magetized disks. In §4, we present the linearized equations and detail the analytical and numerical methods used, with results shown in §5, for a purely toroidal background magnetic field. We also conduct a brief analysis for a purely poloidal background magnetic field and compare the MRI growth rates with the VSI in §6. Finally, we discuss the results in §7 and summarize our main findings in §8.

BASIC EQUATIONS AND EQUILIBRIA
Consider a gaseous, inviscid, magnetized PPD. The gas density, velocity, and magnetic field are denoted by ρ, v, B, respectively. Our formulations are presented in cylindrical (R, φ, Z) coordinates centered on the protostar, although the spherical radius r is also used to simplify expressions. The basic dynamical equations written in SI units are where d/dt ≡ ∂/∂t + v · ∇ is the material derivative and µ0 is the magnetic permeability. The gravitational potential is given by where G is gravitational constant, and M is the mass of the central star. The total pressure is the sum of thermal pressure P and magnetic pressure where B = |B|. The strength of the magnetic field is parametrized by the ratio of gas pressure to magnetic pressure,

Induction Equation
The evolution of magnetic fields is governed by the induction equation, in which three non-ideal MHD effects manifest, Here, σB is resistivity, e is electric charge, and ne is electron number density. The drag coefficient γi represents the momentum transfer between ion-neutral collisions, and the ion density is denoted by ρi. On the right-hand side, the four terms each corresponds to the standard inductive term, Ohmic resistivity, the Hall effect, and ambipolar diffusion. Ohmic resistivity is relevant to electron-neutral collisions. Ambipolar diffusion corresponds to ion-neutral drift. The Hall effect is associated with ion-electron drift, and it differs fundamentally from Ohmic resistivity and ambipolar diffusion as being a non-dissipative process. The relation between the current density J and the magnetic field is completed by the Maxwell equation, The displacement current in Equation (8) is self-consistently neglected for non-relativistic MHD (Balbus 2009). Finally, the magnetic field also satisfies the solenoidal condition, To characterize the non-ideal MHD effects, it is convenient to define the diffusion coefficients for Ohmic resistivity, the Hall effect, and ambipolar diffusion. At a given ionization fraction, the Ohmic diffusivity ηO is independent of field strength and density, whereas the Hall diffusivity ηH ∝ B/ρ, and the ambipolar diffusivity ηA ∝ B 2 /ρ 2 . Hence, ambipolar diffusion dominates in regions of strong fields or low densities, Ohmic resistivity dominates in regions of weak fields or high densities, and the Hall effect governs in between (Lesur 2020). We further introduce dimensionless Elsässer numbers, where the Alfvén velocity vA is given by 2.1.1 Toroidal fields in axisymmetric disks Toroidal magnetic fields have been shown to dominate in global non-ideal MHD simulations, with the saturation of VSI turbulence (e.g Béthune et al. 2017;Bai 2017;Cui & Bai 2020). Furthermore, the VSI is an axisymmetirc instability (N13). These results motivate us to primarily consider purely toroidal fields in axisymmetric disks. In this case, the Hall effect vanishes, and we show here that ambipolar diffusion behaves the same way as Ohmic resistivity. In the limit of B = (0, B φ , 0), we expand the numerator of the last term on the right-hand side of Equation (7), corresponding to ambipolar diffusion, This implies that for electric currents being everywhere perpendicular to magnetic fields, as the case for purely toroidal field geometry, ambipolar diffusion can be treated as an effective resistivity with field dependency (Balbus & Terquem 2001). Consequently, the induction equation reduces to where η = ηO + ηA, and we assume a constant η throughout. The second and final terms on the right-hand side result from the curvilinear geometry.

Effective Energy Equation
The ideal gas law is given by where R is the gas constant, µ is the mean molecular weight, and T is the gas temperature. We define the isothermal sound speed cs, pressure scale-height H, and disk aspect-ratio h through where ΩK = GM/R 3 is the Keplerian angular velocity. In the hydrodynamic limit, isothermality is most favourable for the VSI because the stabilizing effect stabilizing from gas buoyancy is absent (N13; Lin & Youdin 2015). To focus on the effect of magnetic fields, we adopt a locally isothermal thermodynamic response, which is applicable to the outer parts of PPDs wherein the temperature is regulated by stellar irradiation (Chiang & Goldreich 1997). Thereby, the energy equation of gas pressure becomes Consider the evolution of the magnetic pressure associated with a purely azimuthal field, Combining Equations (14), (17), and (18), we formulate an effective energy equation where the effective adiabatic index is defined by Equation (19) resembles the energy equation for a fluid subject to heating or cooling in the hydrodynamic limit (Lin & Youdin 2015). It shows that an axisymmetric, locally isothermal, magnetized gas with a purely azimuthal field behaves like an unmagnetized gas with adiabatic index γB 1 , while it is also subject to non-adiabatic effects on the righthand side: the locally isothermal thermodynamic response, magnetic diffusion, and curvature effects. Note that γB(β) here is not necessarily constant because β(R, Z, t) can be non-uniform and can evolve in time. For β → ∞, we recover a locally isothermal gas with unit adiabatic index γB → 1. On the other hand, a strongly magnetized disk with β → 0 is equivalent to a fluid with an effective adiabatic index γB → 2.

Equilibrium State
We consider axisymmetric steady states with a purely azimuthal velocity and magnetic field. The midplane gas density and temperature are prescribed as where R0 is a reference radius, and qT and qD are constant power-law indices. The equilibrium solutions satisfy the equation of motion, where Ω ≡ v φ /R is the angular frequency. These equations may be solved explicitly for a given magnetic field configuration. We consider two special cases in the following, depending on whether resistivity is included.

Constant-β disks
We consider constant-β disks for ideal MHD (η = 0). The effective energy equation (19) is then satisfied identically so the solutions below represent exact equilibria. In this case, the vertical gradient of equilibrium angular velocity is The density and rotation profiles are As β → ∞, we recover the hydrodynamic limit (N13). The thin disk approximation (h ≪ 1) leads to Ω ≈ ΩK. As β → 0, the density gradient profile shows that a strongly magnetized disk becomes vertically unstratified. A small β 1 leads to strong deviations from the Keplerian rotation, but a minimum value of β is required to ensure Ω 2 > 0 at the disk midplane, For typical PPD parameters, h ∼ 0.05 and qD, qT are of order unity, so the above requirement becomes β O(10 −3 ). This is easily satisfied for the weakly magnetized disks with β 1 that we consider.

Constant-B φ disks
We consider constant-B φ disks for ideal MHD and resistive disks (η = 0). Resistive disks present approximate equilibrium solutions because there is a slow diffusion of the magnetic field due to the global curvature term in Equation (19). However, this is not expected to significantly affect radially localized dynamics such as the VSI. In this case, β is no longer a constant but declines with height, where β0 = β(R, 0). The disk becomes more strongly magnetized with increasing height. For constant-B φ disks, magnetic pressure does not contribute to the vertical equilibrium, and we recover the hydrodynamic limit for density profile The vertical shear gradient differs from constant-β disks because the curvature term in Equation (23) now depends on Z through ρ. We find and the rotation profile is Accordingly, even a strictly isothermal disk (qT = 0) can exhibit vertical shear due to the curvature term. Thus, we may expect a corresponding VSI growth, but this should be examined further in radially global models. For PPDs where magnetic fields are weak, the angular velocity profiles in Equations (27) and (32) converge to the hydrodynamic limit.

SOLBERG-HOILAND CRITERIA
The Solberg-Hoiland criteria describe the linear hydrodynamic stability of ideal fluids against axisymmetric and adiabatic perturbations (Tassoul 1978). Now, our axisymmetric, magnetized disks with purely azimuthal fields obey a similar energy Equation (19) as in hydrodynamics. If, in addition, curvature terms and magnetic tension forces can be neglected, then our disk models satisfy the same form of equations as in adiabatic hydrodynamics. This motivates us to formulate an equivalent Solberg-Hoiland stability criteria for magnetized disks as follows.
The standard hydrodynamic Solberg-Hoiland criteria are expressed in terms of the gradient of pressure P , entropy S ∝ ln P 1/γ /ρ , and angular frequency Ω, where γ is the adiabatic index. By setting P → Π and γ → γB , we find that the stability is ensured in a magnetized disk if where ∂R is the square of epicyclic frequency, and defines the gradient of effective entropy SB of a locally isothermal fluid with an azimuthal field. Correspondingly, we define as the square of buoyancy frequencies in the radial and vertical directions, respectively. To examine whether magnetized disk models are stable subject to the modified Solberg-Hoiland criteria, we first derive the vertical buoyancy frequency for constant-β and constant-B φ thin disks, As β → ∞, N 2 Z ∝ 1/β → 0 for both models, and we recover a locally isothermal hydrodynamic disk with vanishing vertical buoyancy. As β → 0, for constant-β disks, N 2 Z ∝ β and the buoyancy effect vanishes since there is no vertical density stratification in this limit. In a constant-B φ disk, β → 0 occurs for large |Z| as in Equation (29), so that N 2 Z → Ω 2 K Z 2 /H 2 , and the disk expected to be strongly stabilized by magnetic buoyancy.
We now examine whether the first Solberg-Hoiland criterion is satisfied. Equation (33) can be cast into A stable stratification, N 2 R,Z > 0, requires the effective entropy to increase with R or Z because total pressure gradients are usually negative. In thin weakly magnetized disks, (37) and (38)), thus the first criterion is generally satisfied.
Next, we examine the second Solberg-Hoiland criterion. Recall that the Solberg-Hoiland criteria apply to adiabatic flows. Thus, Equations (33) and (34) should only be applied if the governing equations of the magnetized disk (Equations (1)-(2), (19)) can map exactly to adiabatic hydrodynamics. That is, the right-hand side of Equation (19) and the magnetic tension force in Equation (2) should be negligible. This requires (i) a strictly isothermal disk (constant cs in R, Z); (ii) ideal MHD (η = 0); (iii) weak magnetic fields (β ≫ 1) so that curvature terms associated with magnetic fields can be neglected (Pessah & Psaltis 2005).
The first restriction implies ∂Ω/∂Z = 0. The second criterion then becomes which is satisfied in both of our disk models. To study stability with vertical shear, which requires a radially varying disk temperature and non-ideal MHD effects, we must solve the linearized equations explicitly as conducted in the following section.

Perturbation Equations
We consider axisymmetric Eulerian perturbations for v ′ , ρ ′ , and Π ′ of the form The complex frequency is denoted by σ = s + iω, and a real radial wavenumber k is taken. We assume radially localized disturbances, kRR ≫ 1. The background disk variables are evaluated at the reference radius R0, but their vertical dependence is retained. Curvature terms resulting from the cylindrical geometry are neglected, which restricts the analyses to weak magnetic fields (Pessah & Psaltis 2005). By assuming B ′ = (0, B ′ φ , 0) for simplicity, the linearized perturbation equations read Note that the magnetic tension force vanishes for a purely toroidal field and neglecting curvature terms. The perturbed azimuthal magnetic field B ′ φ can be expressed as When β → ∞ and hence Π ′ → P ′ , the set of linearized equations recover the hydrodynamic limit (Lin & Youdin 2015). For finite magnetic field strengths and η = 0, these equations describe the ideal MHD limit. For η > 0, the set of equations describe non-ideal MHD.

Analytical Solutions
Analytic solutions can be obtained in the limit of ideal MHD (η = 0) with a large and constant β, by solving the linearized equations with polynomial solutions. To do so, we make the following simplifying assumptions (Lin & Youdin 2015): (i) A fully radially local approximation (∂/∂R = 0) for background disks. However, the vertical shear that originates from the radial temperature profile (∂T /∂R) is retained.
(ii) We set Ω = ΩK and κ = ΩK where they appear explicitly and without vertical derivatives. As seen in Equation (27), the Keplerian approximation is only valid for large β. A strong magnetic field will lead to substantial deviation from ΩK.
(iii) In the hydrodynamic limit, the VSI is an overstability due to the destabilization of inertial waves (Barker & Latter 2015). We expect a similar result for weak magnetizations and consider low frequency modes with |σ| ≪ ΩK to filter out acoustic waves (Lubow & Pringle 1993).
(iv) We consider thin disks and set To non-dimensionalize the perturbation equations, we choose the appropriate scalings for timescales, lengthscales, and velocities to be the Keplerian orbital time Ω −1 K , background sound speed cs, and pressure scale-height H. We thus write where scaled variables are denoted by an asterisk and are omitted below. Equations (42) can be combined into a single second-order ordinary differential equation, where a1, a2 and a3 are constants in Z and are defined by Equation (50) is equivalent to the hydrodynamic result derived by Lin & Youdin (2015, see their Equation (29)) without cooling and assuming K ≫ 1. This can be seen by . However, they also show that in the hydrodynamic limit, neglecting the radial structure of disk while retaining vertical shear is only valid for gas that is nearly isothermal in its thermodynamic response. We thus expect Equation (50) to only apply for weak magnetizations where γB → 1. In the limit β → ∞ and hence a3 → 0, we recover the equation for the VSI in locally isothermal disks (Lin & Youdin 2015;Barker & Latter 2015). When there is no vertical shear (qT = 0), Equation (50) is analogous to that obtained by Lubow & Pringle (1993) for axisymmetric waves in adiabatic disks. We bring Equation (50) into a form that can be solved analytically (Lubow & Pringle 1993). Define a variable Y (Z) where ξ is a constant to be determined. Plugging this into Equation (50) yields an ordinary differential equation in Y , Then choose ξ to eliminate the Z 2 term, which brings Equation (53) into We proceed to represent solutions of Y by Yn(Z), an nth order polynomial in Z, where Bm are constant coefficients. Substituting this into Equation (55), we arrive at the recurrence relation between Bm and Bm+2, where we have relabled ξ → ξn. For physical solutions to Equation (55), we demand the vertical kinetic energy density of perturbations to remain bound, i.e. ρ|v ′ Z | 2 approaches zero for large |Z|. Then Yn(Z) should be a polynomial, as assumed. Thus Bn+2 = 0 when m = n, or (a1 − 2ξn)n = a2 + ξn.
This relation can be also quickly obtained by recognizing Equation (55) as the Hermite differential equation (Barker & Latter 2015). We choose the negative root in Equation (54) for physical solutions, giving For a given set of disk parameters (q, β, h), Equations (51), (58), and (59) can be readily solved to obtain the dispersion relation for the complex frequency σn = σn(K; n) of the nth mode. The corresponding eigenfunction v ′ Z is then given by Equations (52), (56), and (57). We note that the vertical shear rate increases without bound with height in thin disks (Equation (48)), leading to unbound growth rates (as seen in Figure 1 discussed below), which violates the low frequency approximation.

Numerical Solutions
We also solve the full linearized Equations (42) -(46) numerically. The equations can be written in a form of standard matrix eigenvalue problem, where σ is the eigenvalue, L is a 5 × 5 matrix of linear operators, and T is a vector of eigenfunctions. We use dedalus 2 (Burns et al. 2020), a general purpose spectral code for differential equations, to solve the linear eigenvalue problem. We employ a Chebyshev collocation grid of N = 150 points. Spurious solutions are filtered out and numerical convergence is verified against double-resolution calculations with N = 300 by using the eigentools package 3 . Unlike the analytic model in §4.2, where the disk surface extends to infinity due to the vertically isothermal background state, for numerical solutions we consider a finite vertical domain with Z ∈ [−5, 5]H. The boundary conditions imposed at the disk surfaces are v ′ Z = 0, The former condition is adopted in the limit of ideal MHD, and the latter being additional conditions when resistivity is included. We have experimented with boundary conditions and found that our main findings are insensitive to them.

RESULTS
In this section, we present example solutions in the hydrodynamic limit ( §5.1), the ideal MHD limit ( §5.2), and the non-ideal MHD limit ( §5.3). The fiducial parameter values are h = 0.05, qT = 1, qD = 1.5, and K = 35. Note that qT = 1 gives a constant disk aspect ratio h. We present the non-dimensionlized perturbed quantities as in Equation (49). Plasma beta parameters associated with azimuthal and vertical fields are denoted by β φ = 2µ0P/B 2 φ and βZ = 2µ0P/B 2 Z . We denote midplane plasma beta parameters by β φ0 and βZ0, and midplane Elsässer number Λ0. For reference, values of 2 https://dedalus-project.org/ 3 https://github.com/DedalusProject/eigentools βZ ∼ 10 4 and β φ ∼ 10 2 are found to account for the accretion rate in PPDs (Simon et al. 2013a;Bai 2015). Note, however, as in the discussion above in this section we only consider azimuthal fields, so that β = β φ . Poloidal fields will be explored in Section 6. We follow N13 to denote breathing, corrugation, and surface modes as B, C, and S, respectively, with numbers 1 and 2 representing the fundamental and first overtone modes. For clarity, analytic solutions which have discrete modes are plotted as continuous curves.

The Hydrodynamic Limit
We first compare our numerical and analytical solutions with N13, who considered purely hydrodynamic disks. To this end, we compute the numerical solutions to Equations (42) -(46) and analytical solutions given via Equation (58) in the hydrodynamic limit (β φ → ∞), and compare with numerical solutions to Equation (39) in N13. Note that N13 employed the anelastic approximation (∂ρ/∂t = 0), while we account for full compressibility in numerical solutions. The results are shown in Figure 1.
For the fundamental and first overtone breathing and corrugation modes (B1, B2, C1, and C2), all three methods yield consistent results. For surface modes and higher-order body modes in the lower right of Figure 1, the two numerical solutions also show consistency, especially at low oscillation frequencies. The analytic solutions for these higher-order body modes do not match with numerical solutions due to the lack of a disk surface in the former (Barker & Latter 2015). In practice, the growth rate is limited by the maximum verti-  cal shear rate within the domain (Lin & Youdin 2015), as reflected in the numerical solutions. Overall, the comparison is satisfactory.

The Ideal MHD Limit
In the ideal MHD limit, we examine the behaviour of the VSI modes as functions of disk magnetizations β φ ( §5.2.1), radial wavenumbers K ( §5.2.2), and disk aspect ratios h ( §5.2.3).

Disk magnetization
We now examine the strengths of toroidal magnetic fields on the VSI. In Figure 2, we show example growth rates and fre-quencies for constant-β disks in the left panel and constant-B φ disks in the right panel. Similarly, Figure 3 shows how growth rates of various modes vary with plasma beta. The curves without labels are high-order body modes. We highlight three major findings. Firstly, strong magnetization reduces the VSI growth. Physically, this is because the gas and the magnetic fields are perfectly coupled in the limit of ideal MHD, so that magnetic fields impede the free movement of the perturbed gas. Furthermore, surface modes are the first to vanish with strong magnetization. This can be understood by the fact that the stabilizing vertical buoyancy scales as Z 2 for both models as seen in Equations (37) and (38), hence the gas is subject to stronger stabilization at the disk surface. Finally, the critical β φ to recover hydrodynamic results for a constant-β disk, β φ 10 5 , is smaller than that for the midplane value in a constant-B φ disk, β φ0 10 9 . This is because in a constant-B φ disk, Equation (29) shows that β φ decreases with height, so the vertically averaged β is smaller than its midplane value. Figure 4 shows the flow structure in a constant-β disk. The radial domain of K(R − R0) = 2π corresponds to an interval of 0.18H. The left panels show the fundamental corrugation modes in disks with β φ = 10 5 (top) and β φ = 10 2 (bottom). The perturbed vertical velocities show even symmetry about the midplane. The right panels are corresponding fundamental breathing modes, where the perturbed vertical velocities have odd symmetry. The contours show the magnetic field perturbations Re{B ′ φ exp[iK(R − R0)]}/B φ and is normalized by its maximum value. The perturbed magnetic fields possess opposite symmetry to perturbed vertical velocities. Importantly, we find that strong magnetization confines VSI activity towards the midplane since the stabilizing vertical buoyancy increases with height. The same arguments also apply to a constant-B φ disk.

Radial wavenumber
The left panel of Figure 5 depicts contours of maximum growth rates as a function of β φ and radial wavenumber K for constant-β disks. A critical βc ∼ 10 3 can be defined to separate the disks into two regimes. In constant-B φ disks it is βc ∼ 10 5 . For β φ βc, the maximum growth rate is a monotonically increasing function of K, whereas for β φ βc, the maximum growth rate peaks at some intermediate K. We explain below that βc is in fact the critical disk magnetization below which surface modes are quenched. In Figure 3, we see that at β φ βc, surface modes, which prefer very small radial wavelengths, dominate the maximum growth rates resulting in fast growth rates at large K. At β φ βc, surface modes are suppressed, while body modes that prefer longer radial wavelengths persist, and thus maximum growth rates appear at intermediate radial wavenumbers.

Disk aspect ratio
In the right panel of Figure 5, we show contours of maximum growth rates as functions of β φ and disk aspect ratio h, again for constant-β disks. The maximum growth rates increase with the disk aspect ratio for a given β φ . Requiring modes to fit into the vertical height of the disk, a lower limit can be placed on β φ for the VSI to operate, βmin (−R∂ ln Ω/∂Z) −2 (LP18). This is set by the vertical shear rate, and can be simplified to βmin h −2 using Equation (48). Note, however, this criterion was derived for purely poloidal background fields in a local approximation, while we consider purely toroidal magnetic fields in a vertically global disk. Nevertheless, we find the local condition β φ = h −2 , shown in Figure 5 as the dashed line, successfully predicts the quenching of the VSI in our disk model.

The Non-ideal MHD Limit
In this subsection, we show that the VSI can be revived when non-ideal MHD effects are included. To assure the existence of equilibrium solutions, a constant-B φ disk model is employed ( §2.3.2). With purely toroidal magnetic fields in axisymmetric disks, the three non-ideal MHD effects reduce to only Ohmic resistivity, because the Hall effect vanishes, and ambipolar diffusion acts as an effective resistivity with field dependency ( §2.1.1). Therefore, we only explore the dependency of Ohmic Elsässer number Λ, while we expect the same results apply to ambipolar diffusion. We take the diffusivity η to be constant so that the Elsässer number increases with |Z|, as shown in Equation (11).
In the left panel of Figure 6, we show the growth rates of all unstable modes as a function Λ0 at β φ0 = 10 2 . The labels correspond to surface and body modes in the hydrodynamic limit (Figure 1), which is recovered for small Λ0 or strong Ohmic resistivity. On the other hand, Λ0 → ∞ tends to the ideal MHD limit. We find for Λ0 10 3 , the surface modes vanish and the growth rate of body modes is significantly reduced due to strong magnetization. As Λ0 declines from larger values, the growth rates of these body modes drop at around Λ0 ∼ 10 3 , then they re-emerge and converge to hy-drodynamic results. The growth rates of all modes converge to hydrodynamic results for Λ0 10, with the transition starting at Λ0 ∼ 10 2 . Local analyses demonstrate that the stabilizing effect by magnetic fields will be overcome by magnetic diffusion when Λ h −1 (= 20 in our fiducial disk) for a mode with growth rate ∼ hΩ (LP18) 4 . This is in agreement with our results, though the growth rates from our solutions are only reduced rather than completely suppressed.
The right panel of Figure 6 shows the maximum growth rates as functions of β φ0 and Λ0. For β φ0 > 10 3 , the maximum growth rate is a monotonically decreasing function with increasing Λ0, whereas for β φ0 < 10 3 , the maximum growth rate has its minimum resides at some intermediate Λ0, corresponding to the left panel of Figure 6. The hydrodynamic result is recovered for sufficiently weak fields (β φ0 10 5 ) or sufficiently strong resistivity (λ0 10).

PURELY POLOIDAL BACKGROUND MAGNETIC FIELDS
The above analyses focus on disks threaded by a toroidal magnetic field, which is expected to dominate over poloidal field strengths in PPDs (e.g. Bai 2017; Béthune et al. 2017;Cui & Bai 2020). However, the presence of a poloidal field, even weak, can lead to new effects such as MHD disk winds and the MRI (Bai 2013;Simon et al. 2013b;Gressel et al. 2020). Specifically, the surface layers in outer regions of PPDs are likely sufficiently ionized by stellar FUV radiation to trigger the MRI (Perez-Becker & Chiang 2011;Simon et al. 2013a,b;Bai 2015). These regions are also prone to the VSI since the vertical shear rate increases with height. In this section, we investigate the VSI modes in a disk with purely poloidal magnetic fields, B = (BR, 0, BZ ). In §6.1, we study the effects of poloidal magnetic fields on the VSI in a vertically global disk model. In §6.2, we compare the MRI with the VSI in a local disk model. Although toroidal fields are absent in the background, it is allowed in the perturbed state.

Vertically Global Model
We make several simplifying assumptions to establish disk equilibria with a purely poloidal magnetic field: (i) The Lorentz force in the momentum equation is ignored because the disk is weakly magnetized. We therefore use equilibrium solutions for Ω and ρ in the hydrodynamic limit (Equations (26)-(27) with β → ∞, see also N13).
(ii) We assume thin disks and consider h ≪ 1.
(iii) A constant background vertical magnetic field BZ is assumed, and we seek the required equilibrium radial magnetic field BR, as follows.
The equilibrium magnetic fields must satisfy solenoidal condition and induction equation. Considering only Ohmic resistivity with a constant diffusivity, the equilibrium induction equation is where v φ = RΩφ. In the thin disk approximation, the gradi- The induction equation can be satisfied by the radial field, Notice that this field configuration is not subject to Ohmic diffusion and satisfies the solenoidal condition, Therefore, an approximate equilibrium magnetic field configuration is obtained. The equilibrium solution, Equation (65), resembles Equation (50) in LP18. Since |BR| ∼ O(h)|BZ|, the strength of the magnetic field is dominated by the vertical field, which is assumed a constant so that this is similar to the constant-B φ disks considered in §2.3. With a poloidal field it is not possible to map the problem to adiabatic hydrodynamics. We therefore work with the MHD equations directly. The radial derivatives of background quantities are omitted. The set of linearized equations are Note that there is now a magnetic tension force in the momentum equations. We solve the linear eigenvalue problem, in its dimensionless form, numerically using the spectral method described in §4.3. A resolution of N = 100 is used. Spurious solutions are filtered out with double resolution calculations. Boundary conditions imposed at upper and lower disk surfaces for ideal MHD limit are (Gammie & Balbus 1994;Sano & Miyama 1999) Additional conditions are imposed with Ohmic resistivity, The quantities σ, Z, K reported below are nondimensionlized as in Equation (49). With the inclusion of vertical magnetic fields, MRI modes emerge in the numerical solutions. Unlike the VSI, however, MRI modes are not overstable even with resistivity, i.e. Im(σ) ≡ ω = 0 (Sano & Miyama 1999); while it can be seen in Figure 2 that VSI modes generally have 0.01 < ω < 1, hence the MRI modes do not contaminate the numerical results of the VSI.
The left panel of Figure 7 shows growth rates of all unstable modes as a function of disk magnetization βZ0 in the ideal MHD limit. Consistent with purely the toroidal field model, strong magnetization suppresses VSI growth and surface modes are the first to vanish with increasing field strengths. The critical βZ0 to recover hydrodynamic results is even larger, βZ0 10 10 , because the total magnetic field strength increases over height as radial magnetic field develops away from the midplane in Equation (65).
Growth rates including Ohmic resistivity are shown in the middle and right panels of Figure 7. The middle panel shows the growth rates of all unstable modes as a function of Ohmic Elsässer number Λ0 at βZ0 = 10 4 . Nearly all the VSI modes are suppressed when Ohmic resistivity is weak at Λ0 10. The VSI modes start to grow when Λ0 10. Fundamental body modes B1 and C1 converge to the hydrodynamic growth rates at Λ0 ∼ 1. On the other hand, surface modes and high order body modes show slower transitions to their hydrodynamic growth rates, requiring Λ0 ∼ 0.1. The right panel of Figure 7 shows the maximum growth rate as a function of βZ0 and Λ0. Large Λ0 ∼ 10 2 corresponds to the ideal MHD limit, where maximum growth rate drops to ∼ 5 × 10 −3 for βZ0 = 10 4 . It can be seen that Λ0 < 1 is required for the fastest growing modes to recover hydrodynamic results in a wide range of βZ0 from 10 4 to 10 10 .
The above results are similar to that for toroidal fields, which indicates that the qualitative effect of a magnetic field and resistivity on the VSI does not depend on the background field geometry.

VSI VS MRI
To better compare the VSI and the MRI, we perform a vertically local linear analysis in the incompressible limit with Ohmic resistivity. In this local model, background quantities are assumed to be uniform and its values set to that in the vertically global model ( §6.1) at some fiducial height. We take the vertical shear rate R∂Ω 2 /∂Z as an input parameter. We consider axisymmetric perturbations that are proportional to exp(σt + ikRR + ikZZ). The wavenumber vector is denoted by k = (kR, 0, kZ ).
The above equations give a dispersion relation where A = kR kZ R∂Ω 2 ∂Z .
Equation (84) generalizes that of LP18 to include Ohmic re- sistivity. Appendix A explores this dispersion relation in more detail in the limit |kRBR| ≪ |kZBZ|. Here, we solve Equation (84) in full to investigate the dominance of MRI and VSI as functions of disk magnetizations βZ , radial wavenumbers K, and Elsässer numbers Λ. A vertical wavenumber kZ = 2π/H is fixed in numerical calculations. VSI modes require kR and kZ with opposite signs if ∂Ω/∂Z < 0 (Nelson et al. 2013;Latter & Papaloizou 2018), as we will consider below, hence we use |K| to denote the absolute value of K. Other background quantities are evaluated at Z = 2H. The left panel of Figure 8 shows local growth rates of MRI modes in the ideal MHD limit (η = 0) without vertical shear (R∂Ω 2 /∂Z = 0). In the white regions, the MRI is quenched by strong magnetizations. It can be seen that the MRI modes prefer small |K|, with a maximum growth rate of s = 0.75 at βZ = 137. In the right panel of Figure 8, we show growth rates of VSI modes by setting R∂Ω 2 /∂Z = −0.1. The fastest growing VSI modes prefer weak magnetization and large |K|, in contrast to fast growing MRI modes. Local VSI modes with K = 150 and βZ = 10 9 have a growth rate s = 0.05, which is less than that of the fastest growing surface mode s = 0.094 obtained from vertically global analysis shown in the right panel of Figure 3.
We next consider resistive disks. In the left panel of Figure  9, we show growth rates as functions of βZ and Λ for MRI modes by setting K = 0. In contrast to the VSI, MRI growth rates decline towards small Λ as the MRI is dampened. In the right panel of Figure 9, we show growth rates of mostly VSI modes by setting K = 150. For βZ 10 5 , we obtain VSI growth rates in the hydrodynamic limit. For βZ 10 5 , the VSI is dampened for Λ 10, but revived with small Λ or strong resistivity.
Finally, in Figure 10, we show the ratio of MRI growth rates obtained from the left panel of Figure 9, to VSI growth rates obtained from the right panel of Figure 9. For βZ 10 5 the VSI dominates over the MRI. For βZ 10 5 , which includes PPDs with typical βZ ∼ 10 4 , the VSI dominates if Λ 0.09.

DISCUSSION
7.1 Comparison with previous works LP18 carried out local linear analyses of the VSI in the ideal MHD limit with an exact background equilibrium solution for a purely poloidal field. In this work, we set up global equilibria for a purely azimuthal field. Comparing our numerical solutions to local analytical results, we find that the overall mode behaviour that magnetization can stabilize the VSI is in agreement with each other. Global numerical simulations of magnetized PPDs carried out by Cui & Bai (2020) indeed show that magnetism tends to suppress the VSI growth, whereas ambipolar diffusion acts to revive VSI modes. Their simulations also show the absence of surface modes, which is consistent with our findings that surface modes are the first to be dampened with increasing field strengths.
Global simulations of PPDs initially threaded by a largescale poloidal magnetic field suggest a field configuration dominated by the toroidal component once the disk reaches a quasi-steady state (e.g. Bai 2017; Béthune et al. 2017;Cui & Bai 2020). This suggests that, as far as the VSI is concerned, it is the disk model with an azimuthal field that is more relevant, as employed in most of this paper. Furthermore, the magnetization in our toroidal field model is parametrized by β, which does not depend on the orientation of the magnetic field. Our results are thus applicable to the aforementioned simulations wherein the toroidal field reverses polarity across the disk midplane because of the Keplerian shear.

Application to PPDs
In the outer part of the PPDs ( 30 AU), ambipolar diffusion is the dominant non-ideal MHD effect, with Elsässer numbers approximately unity. For a purely azimuthal field, ambipolar diffusion acts as an effective resistivity with field dependency (see §2.1.1). Hence, our results for Ohmic resistivity are also applicable to ambipolar diffusion. A value of Λ0 = 1 and β φ0 = 10 2 gives a maximum growth rate of s = 0.087 ( Figure  6), which is close to the hydrodynamic result, s = 0.094ΩK. In the inner part of the PPDs, Ohmic resistivity becomes the dominant non-ideal MHD effect, though the Hall effect also contribute. At 2 AU, the midplane Λ0 = 5 × 10 −4 (Bai 2017) gives a maximum growth rate of s = 0.094 for a wide range from β φ0 = 10 to β φ0 = 10 9 (Figure 6), as a small Λ0 enables the recovery of hydrodynamic results.
Our locally isothermal disk models, which correspond to instantaneous cooling, favor the VSI because there is no stabilizing effect from vertical gas buoyancy (N13, Lin & Youdin 2015). However, in a realistic PPD cooling timescales are finite and is sensitive to stellar irradiation and dust properties (Malygin et al. 2017;Pfeil & Klahr 2020;Flock et al. 2020). In magnetized disks, the magnetic field will provide extra stabilization via magnetic buoyancy in addition to gas buoyancy. Thus, we expect that in PPDs the required cooling time may be shorter than that estimated based on purely hydrodynamic models (Lin & Youdin 2015), unless Λ0 is small enough to diminish the stabilizing effect from magnetic fields. Detailed analyses should be conducted to give new critical cooling timescales in magnetized disks.
The analysis we present with Ohmic resistivity points to future directions in including ambipolar diffusion and the Hall effect, the latter of which requires a poloidal field. A few obstacles and complications needs to be resolved when incorporating these two non-ideal MHD effects. Firstly, it is difficult to find appropriate background equilibria in a global model because of the vertical shear, especially for in presence of poloidal magnetic fields (Ogilvie 1997). Furthermore, ambipolar diffusion gives rise to anisotropic damping and introduces the ambipolar shear instability (Blaes & Balbus 1994;Kunz & Balbus 2004;Kunz 2008). The Hall effect will further introduce the Hall shear instability, and its effect is polarity dependent (Balbus & Terquem 2001;Kunz 2008). All of these effects will complicate the problem and deserve step-by-step analyses in local and global disk models in the future.

Implications to dust dynamics
Small dust grains tend to settle towards the disk midplane (Dubrulle et al. 1995). However, the VSI drives turbulence that can vertically mix up dust particles (Stoll & Kley 2016;Flock et al. 2017Flock et al. , 2020. On the other hand, our results show that strong magnetization can stabilize the VSI away from the disk midplane, implying a limit on the vertical extent of the ensuing VSI turbulence and therefore a maximum dust layer thickness, H d,max . For definiteness, consider a purely toroidal field with constant β and neglect non-ideal MHD effects. The destabilizing vertical shear, R∂Ω 2 /∂Z, competes against the stabilizing magnetic buoyancy, N 2 Z . We therefore expect the VSI to be suppressed where R∂Ω 2 /∂Z /N 2 Z < ζ, where ζ is some critical ratio 5 . Using Equations (25) and (37), we estimate assuming β ≫ 1 and a thin disk. The example in Figure 4 with β = 10 2 , h = 0.05, and qT = 1 show that gas motions are negligible for |Z| 2H, which suggest ζ ≃ 2.5. One may ask if H d,max can be made sufficiently small to constrain particles to a dense midplane layer that can undergo, for example, the streaming instability and hence facilitate planetesimal formation (Youdin & Goodman 2005b;Johansen et al. 2009). For dynamical growth, the streaming instability requires a local dust-to-gas mass ratio ν 1. The metallicity is Σ d /Σg = νH d /H, where Σ d and Σg are the dust and gas surface densities, respectively, and H d is the characteristic dust scale height. Thus, if we take H d = H d,max /2, then Inserting typical PPD values (qT , h, Σ d /Σg) ≃ (1, 0.05, 0.01), we find an equipartition field strength, β = 1, would be required to confine dust into a thin layer such that ν ∼ 1 by quenching the VSI elsewhere. Such a strong field is unrealistic for PPDs (e.g. Simon et al. 2013a;Bai 2015). This suggests that magnetic fields do not affect the vertical dust structure in PPDs through its geometric effect on the VSI.
Instead, magnetic effects likely manifest through weakening the VSI and hence the ensuing turbulence, as found in this work and Cui & Bai (2020). This determines H d ≃ δZ/StHg (Dubrulle et al. 1995), where St is the particle Stokes number and δZ is the dimensionless vertical diffusion coefficient associated with VSI turbulence. We expect δZ to drop with larger β and δZ to increase with non-ideal MHD effects. This relation should be calibrated with future simulations, which can then be used to estimate magnetic field strengths from the vertical distribution of dust in PPDs.

CONCLUSIONS
In this work, we perform linear analyses of the VSI under the ideal MHD limit and with Ohmic resistivity. A vertically global and radially local disk model is employed to properly accommodate the characteristic VSI modes of elongated vertical wavelengths. A locally isothermal equation of state is assumed to better focus on the effect of magnetism. Our main findings are summarized as follows.
• In the ideal MHD limit, magnetic fields operate as a stabilizing effect to suppress the growth of VSI modes. Surface modes are the first to vanish rather than body modes with increasing magnetic field strengths. Ohmic resistivity acts as destabilizing effect to assist the VSI growth.
• In weakly magnetized disks, surface modes show maximum growth rates at large radial wavenumbers, while in strongly magnetized disks, surface modes are dampened, while body modes are dominant and prefer intermediate radial wavenumbers. Large disk aspect ratios or vertical shear rates leads to fast VSI growth.
• The MRI modes appear when a poloidal magnetic field is present. In the local analysis, we find that MRI and VSI modes dominate at different βZ and Λ in the ideal MHD limit. MRI prefers relatively strong disk magnetizations and small radial wavenumbers. The VSI modes are most effective at weak magnetizations and large radial wavenumbers. With Ohmic resistivity, a typical value of βZ = 10 4 in PPDs results in a critical Λ 0.09 for the dominance of the VSI.

DATA AVAILABILITY
The data underlying this article will be shared on reasonable request to the corresponding author.

APPENDIX A: PROPERTIES OF THE LOCAL DISPERSION RELATION
We explore the properties of the local dispersion relation Equation (84) in the limit of |kRBR| ≪ |kZ BZ|. When there is no vertical shear, R∂Ω 2 /∂Z = 0, Equation (84) is identical to the dispersion relation for the MRI in Equation (22) of Sano & Miyama (1999). When there is no magnetic field, vAZ = 0 and η = 0, our dispersion relation recovers Equation (34) of Goldreich & Schubert (1967) in the case when Brunt-Väisälä frequency vanishes. When there is a magnetic field but η = 0, the dispersion relation recovers Equation (58) of LP18.
The full dispersion relation (84) is implicit in σ and k. To obtain the maximum growth rate, we take the derivative ∂/∂k for each term on the left-hand side and assume ǫ to be a constant. This yields a relation between the maximum growth rate σmax and the most unstable wavenumber k, The maximum growth rate is computed numerically by substituting k in Equation (84), and the results are shown in Figure A1. In the limit of weak resistivity η → 0, the maximum growth rate and the most unstable wavenumber are When there is no vertical shear A = 0, these expressions recover the MRI channel modes. When η → ∞, and we recover Sano & Miyama (1999) when A = 0. This paper has been typeset from a T E X/L A T E X file prepared by the author.