On The Vertical Shear Instability In Magnetized Protoplanetary Discs | Monthly Notices Of The Royal Astronomical Society

Can Cui, Min-Kai Lin, On the vertical shear instability in magnetized protoplanetary discs, Monthly Notices of the Royal Astronomical Society, Volume 505, Issue 2, August 2021, Pages 2983–2998, https://doi.org/10.1093/mnras/stab1511

The vertical shear instability (VSI) is a robust phenomenon in irradiated protoplanetary discs (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 magnetohydrodynamic effects should be properly considered. To this end, we present linear analyses of the VSI in magnetized discs 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, a locally isothermal equation of state is employed. 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 overcomes magnetic stabilization, and hydrodynamic results are recovered. We also briefly consider poloidal magnetic fields to account for the magnetorotational instability (MRI), which may develop towards surface layers in the outer regions of PPDs. The MRI grows efficiently at small radial wavenumbers, in contrast to the VSI. When resistivity is considered, the VSI dominates over the MRI for Ohmic Elsässer numbers ≲0.09 at plasma beta parameter βZ ∼ 104.

It has been postulated over decades that the turbulence and angular momentum transport in most astrophysical accretion discs are mediated by the magnetorotational instability (MRI; Chandrasekhar 1961; Balbus & Hawley 1991). However, protoplanetary discs (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 disc (Bai & Stone 2011; Perez-Becker & Chiang 2011; Simon et al. 2013a, b). Instead, angular momentum transport is dominated by magnetized disc winds, leaving the bulk disc mostly laminar (Bai & Stone 2013; Gressel et al. 2015, 2020; Bai et al. 2016; Bai 2017; Béthune, Lesur & Ferreira 2017).

Nevertheless, some level of turbulence is expected in PPDs to account for the recent Atacama Large Millimeter/submillimeter Array (ALMA) observations of molecular line emissions (Teague et al. 2016; Flaherty et al. 2017, 2018, 2020). Furthermore, turbulence may serve as an essential ingredient in many stages of planet formation. Turbulence affects the gravitational sedimentation (Dubrulle, Morfill & Sterzik 1995; Johansen & Klahr 2005; Youdin & Lithwick 2007), radial diffusion (Clarke & Pringle 1988), and collisional growth (Ormel & Cuzzi 2007; Birnstiel, Dullemond & Brauer 2010) of dust particles. Long-lived vortices induced by turbulence (e.g. Raettig, Klahr & Lyra 2015; Manger & Klahr 2018) can concentrate dust particles (Barge & Sommeria 1995; Klahr & Henning 1997; Cuzzi, Hogan & Shariff 2008) and seed planetesimal formation through streaming instability or gravitational instability (Youdin & Goodman 2005; 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, Estrada & Cuzzi 2020). Turbulence also influences the radial migration of planets, as well as the flow morphology and gap formation around them (Nelson & Papaloizou 2004; Papaloizou, Nelson & Snellgrove 2004). Therefore, understanding the origin and characteristics of turbulence in PPDs is essential to many aspects of planet formation and evolution.

The lack of magnetohydrodynamic (MHD) turbulence in a PPD led to a surge in the interest of purely hydrodynamic instabilities (Lyra & Umurhan 2019; Weiss, Bai & Fu 2021). Among the most explored are the vertical shear instability (VSI; Nelson, Gressel & Umurhan 2013, hereafter N13; Lin & Youdin 2015; Latter & Papaloizou 2018, hereafter LP18; Cui & Bai 2020), the convective overstability in its linear (Klahr & Hubbard 2014; Lyra 2014; Latter 2016) and non-linear (subcritical baroclinic instability; Klahr & Bodenheimer 2003; Petersen, Julien & Stewart 2007a; Petersen, Stewart & Julien 2007b; Lesur & Papaloizou 2010) phases, and the zombie vortex instability (Marcus et al. 2013, 2015; Lesur & Latter 2016; Umurhan, Shariff & Cuzzi 2016b). These instabilities set in under certain thermodynamic and structural conditions, thereby operating at distinct regions of PPDs (Malygin et al. 2017; Lyra & Umurhan 2019; Pfeil & Klahr 2019). The VSI is of particular interest as it extends a large portion of the disc (e.g. Lin & Youdin 2015; Lyra & Umurhan 2019).

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 discs 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 N13. It quickly drew intensive interests (e.g. Stoll & Kley 2014; Barker & Latter 2015; Umurhan, Nelson & Gressel 2016a; LP18; Lin 2019; Cui & Bai 2020; Flock et al. 2020; Schäfer, Johansen & Banerjee 2020; Flores-Rivera et. al. 2020) and is considered to be a promising hydrodynamic mechanism in driving turbulence in PPDs.

A differentially rotating disc 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 disc 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 disc 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 mid-plane ( 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 time-scales are less than 0.01–0.1 times the local dynamical time-scales, 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 yields a Shakura–Sunyaev (Shakura & Sunyaev 2009) α value on the order of 10−4–10−3 ( N13; Stoll & Kley 2014; Cui & Bai 2020). Non-axisymmetric 3D simulations show the development of vortices (Richard, Nelson & Umurhan 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 dust-trapping 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 discs ( 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 discs, LP18 estimate a critical Ohmic Elsässer number of ∼1/h, where h is the disc aspect ratio, for the VSI to operate, by requiring the magnetic diffusion time-scale to be shorter than the Alfvén wave propagation time-scales. Non-linear MHD simulations, applicable to outer regions of the disc, demonstrate that the VSI can initiate and sustain turbulence in magnetized discs (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 discs with non-ideal effects is still lacking. Such analyses can be useful for understanding the VSI mode properties in real PPDs and interpreting non-linear simulations. Hence, in this work, we extend the linear stability analysis of the VSI to weakly ionized gas in a vertically global disc 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 Section 2, we introduce the basic formulation and establish the equilibrium state of the problem. In Section 3, we derive and discuss the Solberg–Hoiland stability criteria for magnetized discs. In Section 4, we present the linearized equations and detail the analytical and numerical methods used, with results shown in Section 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 Section 6. Finally, we discuss the results in Section 7 and summarize our main findings in Section 8.

The Solberg–Hoiland criteria describe the linear hydrodynamic stability of ideal fluids against axisymmetric and adiabatic perturbations (Tassoul 1978). Now, our axisymmetric, magnetized discs 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 disc models satisfy the same form of equations as in adiabatic hydrodynamics. This motivates us to formulate an equivalent Solberg–Hoiland stability criterion for magnetized discs as follows.

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 disc (equations 1– 2, and 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

a strictly isothermal disc (constant cs in R, Z);

weak magnetic fields (β ≫ 1) so that curvature terms associated with magnetic fields can be neglected (Pessah & Psaltis 2005).

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).

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.

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).

Comparison of growth rates s and oscillation frequencies ω of N13 and this work (βϕ → ∞). Line denotes analytical solutions obtained in Section 4.2. Diamonds denote numerical solutions obtained in N13. Crosses denote numerical solutions obtained in Section 4.3. Labels B, C, and S represent breathing, corrugation, and surface modes, respectively. Numbers represent fundamental and first overtone modes. Modes reside in the lower right are high-order body modes.

In this section, we present example solutions in the hydrodynamic limit (Section 5.1), the ideal MHD limit (Section 5.2), and the non-ideal MHD limit (Section 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 disc 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 |$\beta _{\phi } =2\mu _0 P/B_{\phi }^2$| and |$\beta _Z =2\mu _0 P/B_Z^2$|⁠ . We denote mid-plane plasma beta parameters by βϕ0 and βZ0, and mid-plane Elsässer number Λ0. For reference, values of βZ ∼ 104 and βϕ ∼ 102 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 that have discrete modes are plotted as continuous curves.

We first compare our numerical and analytical solutions with N13, who considered purely hydrodynamic discs. 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 (⁠|$\partial$| ρ/|$\partial$| t = 0), while we account for full compressibility in numerical solutions. The results are shown in Fig. 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 Fig. 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 disc surface in the former (Barker & Latter 2015). In practice, the growth rate is limited by the maximum vertical shear rate within the domain (Lin & Youdin 2015), as reflected in the numerical solutions. Overall, the comparison is satisfactory.

In the ideal MHD limit, we examine the behaviour of the VSI modes as functions of disc magnetizations βϕ (Section 5.2.1), radial wavenumbers K (Section 5.2.2), and disc aspect ratios h (Section 5.2.3).

We now examine the strengths of toroidal magnetic fields on the VSI. In Fig. 2, we show example growth rates and frequencies for constant-β discs in the left-hand panel and constant-Bϕ discs in the right-hand panel. Similarly, Fig. 3 shows how growth rates of various modes vary with plasma beta. The curves without labels are high-order body modes.

Growth rates s and oscillation frequencies ω of unstable modes at discrete plasma βϕ or βϕ0. Left: a constant-β disc. Right: a constant-Bϕ disc. Curves denote analytic solutions (Section 4.2) and crosses denote numerical solutions (Section 4.3).

Growth rates s of all unstable modes as a function of plasma βϕ or βϕ0. Left: a constant-β disc. Right: a constant-Bϕ disc.

We highlight three major findings. First, 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 Z2 for both models as seen in equations ( 37) and ( 38), hence the gas is subject to stronger stabilization at the disc surface. Finally, the critical βϕ to recover hydrodynamic results for a constant-β disc, βϕ ≳ 105, is smaller than that for the mid-plane value in a constant-Bϕ disc, βϕ0 ≳ 109. This is because in a constant-Bϕ disc, equation ( 29) shows that βϕ decreases with height, so the vertically averaged β is smaller than its mid-plane value.

Fig. 4 shows the flow structure in a constant-β disc. The radial domain of K(R − R0) = 2|$\pi$| corresponds to an interval of 0.18H. The left-hand panels show the fundamental corrugation modes in discs with βϕ = 105 (top) and βϕ = 102 (bottom). The perturbed vertical velocities show even symmetry about the mid-plane. The right-hand panels are corresponding fundamental breathing modes, where the perturbed vertical velocities have odd symmetry. The contours show the magnetic field perturbations |$\Re \lbrace {B}_\phi ^{\prime } \, \exp [\mathrm{i}K(R-R_0)]\rbrace /B_\phi$| 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 mid-plane since the stabilizing vertical buoyancy increases with height. The same arguments also apply to a constant-Bϕ disc.

Fundamental corrugation modes C1 (left) and fundamental breathing modes B1 (right) at βϕ = 105 (top) and βϕ = 102 (bottom) in constant-β discs. Contours show the magnetic field perturbations, |$\Re \lbrace {B}_\phi ^{\prime } \, \exp [\mathrm{i}K(R-R_0)]\rbrace /B_\phi$|⁠ , normalized by its maximum. Arrows denote perturbed velocity vectors (⁠|${v}_R^{\prime }, {v}_Z^{\prime }$|⁠ ). The radial interval of K(R − R0) = 2|$\pi$| corresponds to 0.18H.

The left-hand panel of Fig. 5 depicts contours of maximum growth rates as a function of βϕ and radial wavenumber K for constant-β discs. A critical βc ∼ 103 can be defined to separate the discs into two regimes. In constant-Bϕ discs it is βc ∼ 105. 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 disc magnetization below which surface modes are quenched. In Fig. 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.

Contours of maximum growth rates in logarithmic scale as functions of βϕ and radial wavenumber K (left) or disc aspect ratio h (right) in a constant-β disc. Dashed line represents βϕ = h−2.

In the right-hand panel of Fig. 5, we show contours of maximum growth rates as functions of βϕ and disc aspect ratio h, again for constant-β discs. The maximum growth rates increase with the disc aspect ratio for a given βϕ. Requiring modes to fit into the vertical height of the disc, a lower limit can be placed on βϕ for the VSI to operate, βmin ≳ (−R|$\partial$| ln Ω/|$\partial$| 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 disc. Nevertheless, we find the local condition βϕ = h−2, shown in Fig. 5 as the dashed line, successfully predicts the quenching of the VSI in our disc model.

In this section, 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ϕ disc model is employed (Section 2.3.2). With purely toroidal magnetic fields in axisymmetric discs, 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 (Section 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-hand panel of Fig. 6, we show the growth rates of all unstable modes as a function Λ0 at βϕ0 = 102. The labels correspond to surface and body modes in the hydrodynamic limit (Fig. 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 |$\Lambda _0\, \gtrsim 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 ∼ 103, then they re-emerge and converge to hydrodynamic results. The growth rates of all modes converge to hydrodynamic results for Λ0 ≲ 10, with the transition starting at Λ0 ∼ 102. Local analyses demonstrate that the stabilizing effect by magnetic fields will be overcome by magnetic diffusion when Λ ≲ h−1 (=20 in our fiducial disc) 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 effect of Ohmic resistivity on the VSI growth rate. Left: growth rate s versus Ohmic Elsässer number Λ0 at βϕ0 = 102. Right: contour of maximum growth rates in logarithmic scales as functions of plasma βϕ0 and Ohmic Elsässer number Λ0.

The right-hand panel of Fig. 6 shows the maximum growth rates as functions of βϕ0 and Λ0. For βϕ0 > 103, the maximum growth rate is a monotonically decreasing function with increasing Λ0, whereas for βϕ0 < 103, the maximum growth rate has its minimum resides at some intermediate Λ0, corresponding to the left-hand panel of Fig. 6. The hydrodynamic result is recovered for sufficiently weak fields (βϕ0 ≳ 105) or sufficiently strong resistivity (λ0 ≲ 10).

The above analyses focus on discs 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 disc 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 far-ultraviolet (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 disc with purely poloidal magnetic fields, |$\boldsymbol {B}=(B_R,0,B_Z)$|⁠ . In Section 6.1, we study the effects of poloidal magnetic fields on the VSI in a vertically global disc model. In Section 6.2, we compare the MRI with the VSI in a local disc model. Although toroidal fields are absent in the background, it is allowed in the perturbed state.

We make several simplifying assumptions to establish disc equilibria with a purely poloidal magnetic field.

The Lorentz force in the momentum equation is ignored because the disc is weakly magnetized. We therefore use equilibrium solutions for Ω and ρ in the hydrodynamic limit (equations 26 and 27 with β → ∞, see also N13).

We assume thin discs and consider h ≪ 1.

A constant background vertical magnetic field BZ is assumed, and we seek the required equilibrium radial magnetic field BR, as follows.

The quantities σ, Z, and K reported below are non-dimensionlized 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. ℑ(σ) ≡ ω = 0 (Sano & Miyama 1999), while it can be seen in Fig. 2 that VSI modes generally have 0.01 < ω < 1, hence the MRI modes do not contaminate the numerical results of the VSI.

The left-hand panel of Fig. 7 shows growth rates of all unstable modes as a function of disc 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 ≳ 1010, because the total magnetic field strength increases over height as radial magnetic field develops away from the mid-plane in equation ( 65).

Left: growth rates s of all unstable modes as a function of mid-plane βZ0 in the ideal MHD limit. Middle: growth rates s versus mid-plane Ohmic Elsässer number Λ0 at βZ0 = 104. Right: contour of maximum growth rates as functions of βZ0 and Ohmic Elsässer number Λ0.

Growth rates including Ohmic resistivity are shown in the middle and right-hand panels of Fig. 7. The middle panel shows the growth rates of all unstable modes as a function of Ohmic Elsässer number Λ0 at βZ0 = 104. Nearly all the VSI modes are suppressed when Ohmic resistivity is weak at Λ0 ≳ 10. The VSI modes start to grow when |$\Lambda _0 \lesssim 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-hand panel of Fig. 7 shows the maximum growth rate as a function of βZ0 and Λ0. Large Λ0 ∼ 102 corresponds to the ideal MHD limit, where maximum growth rate drops to ∼5 × 10−3 for βZ0 = 104. 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 104 to 1010.

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.

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 (Section 6.1) at some fiducial height. We take the vertical shear rate R|$\partial$| Ω2/|$\partial$| Z as an input parameter. We consider axisymmetric perturbations that are proportional to exp (σt + ikRR + ikZZ). The wavenumber vector is denoted by |$\boldsymbol {k}=(k_R,0,k_Z)$|⁠ .

The left-hand panel of Fig. 8 shows local growth rates of MRI modes in the ideal MHD limit (η = 0) without vertical shear (R|$\partial$| Ω2/|$\partial$| 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-hand panel of Fig. 8, we show growth rates of VSI modes by setting R|$\partial$| Ω2/|$\partial$| 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 = 109 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-hand panel of Fig. 3.

Contours of growth rates as functions of βZ and radial wavenumber |K| by equation ( 84) in the ideal MHD limit. The vertical wavenumber is set to be kZ = 2|$\pi$| /H, and |$\epsilon ={k}_{Z}^2/{k}^2$|⁠ . Left: no vertical shear R|$\partial$| Ω2/|$\partial$| Z = 0. Right: with vertical shear R|$\partial$| Ω2/|$\partial$| Z = −0.1.

We next consider resistive discs. In the left-hand panel of Fig. 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-hand panel of Fig. 9, we show growth rates of mostly VSI modes by setting K = 150. For βZ ≳ 105, we obtain VSI growth rates in the hydrodynamic limit. For βZ ≲ 105, the VSI is dampened for Λ ≳ 10, but revived with small Λ or strong resistivity.

Contours of growth rates as functions of β and Ohmic Elsässer number Λ by equation ( 84). The vertical wavenumber and vertical shear are set to be kZ = 2|$\pi$| /H and R|$\partial$| Ω2/|$\partial$| Z = −0.1. Left: MRI modes (|K| = 0). Right: VSI modes (|K| = 150).

Finally, in Fig. 10, we show the ratio of MRI growth rates obtained from the left-hand panel of Fig. 9 to VSI growth rates obtained from the right-hand panel of Fig. 9. For βZ ≳ 105 the VSI dominates over the MRI. For βZ ≲ 105, which includes PPDs with typical βZ ∼ 104, the VSI dominates if Λ ≲ 0.09.

The ratio of MRI growth rates to VSI growth rates computed from local dispersion relation equation ( 84).

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 large-scale poloidal magnetic field suggest a field configuration dominated by the toroidal component once the disc 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 disc 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 disc mid-plane because of the Keplerian shear.

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 Section 2.1.1). Hence, our results for Ohmic resistivity are also applicable to ambipolar diffusion. A value of Λ0 = 1 and βϕ0 = 102 gives a maximum growth rate of s = 0.087 (Fig. 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 contributes. At 2 au, the mid-plane Λ0 = 5 × 10−4 (Bai 2017) gives a maximum growth rate of s = 0.094 for a wide range from βϕ0 = 10 to 109 (Fig. 6), as a small Λ0 enables the recovery of hydrodynamic results.

Our locally isothermal disc models, which correspond to instantaneous cooling, favour the VSI because there is no stabilizing effect from vertical gas buoyancy ( N13; Lin & Youdin 2015). However, in a realistic PPD cooling time-scales are finite and is sensitive to stellar irradiation and dust properties (Malygin et al. 2017; Flock et al. 2020; Pfeil & Klahr 2020). In magnetized discs, 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 time-scales in magnetized discs.

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 need to be resolved when incorporating these two non-ideal MHD effects. First, 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 disc models in the future.

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_\mathrm{d} \simeq \sqrt{\delta _Z/\mathrm{St}}H_\mathrm{g}$| (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 stronger magnetization and to increase with stronger 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.

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 disc 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 discs, surface modes show maximum growth rates at large radial wavenumbers, while in strongly magnetized discs, surface modes are dampened, while body modes are dominant and prefer intermediate radial wavenumbers. Large disc aspect ratios or vertical shear rates lead 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 disc 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 = 104 in PPDs results in a critical Λ ≲ 0.09 for the dominance of the VSI.

We are pleased to thank the anonymous referee for helpful comments and Xue-Ning Bai, Henrik Latter, and Gordon Ogilvie for fruitful discussions. CC acknowledges the support from Department of Applied Mathematics and Theoretical Physics at University of Cambridge. M-KL is supported by the Ministry of Science and Technology of Taiwan under grant 107-2112-M-001-043-MY3 and an Academia Sinica Career Development Award (AS-CDA-110-M06). Software used: dedalus (Burns et al. 2020) and eigentools package.

The data underlying this paper will be shared on reasonable request to the corresponding author.

We remark that this analogy can be generalized to locally polytropic discs where |$P=K\rho ^\Gamma$|⁠ , where K is a prescribed function of position and Γ is the polytropic index. In this case the effective adiabatic index becomes γB = (Γβ + 2)/(β + 1).

https://github.com/DedalusProject/eigentools

Equation (64) of LP18 contains a typographical error, the corrected expression is |$E_\eta \lesssim 1/q$| (Latter, private communication).

For the hydrodynamic VSI, ζ ∼ Ωtcool, where tcool is the cooling time-scale (Lin & Youdin 2015).

Arlt R. , Urpin V. , 2004, A&A, 426, 755 10.1051/0004-6361:20035896

Armitage P. J. , 2011, ARA&A, 49, 195 10.1146/annurev-astro-081710-102521

Bai X.-N. , 2011, ApJ, 739, 51 10.1088/0004-637X/739/1/51

Bai X.-N. , 2013, ApJ, 772, 96 10.1088/0004-637X/772/2/96

Bai X.-N. , 2015, ApJ, 798, 84 10.1088/0004-637X/798/2/84

Bai X.-N. , 2017, ApJ, 845, 75 10.3847/1538-4357/aa7dda

Bai X.-N. , Stone J. M. , 2011, ApJ, 736, 144 10.1088/0004-637X/736/2/144

Bai X.-N. , Stone J. M. , 2013, ApJ, 769, 76 10.1088/0004-637X/769/1/76

Bai X.-N. , Ye J. , Goodman J. , Yuan F. , 2016, ApJ, 818, 152 10.3847/0004-637X/818/2/152

Balbus S. A. , 2009, preprint (arXiv:0906.0854)

Balbus S. A. , Hawley J. F. , 1991, ApJ, 376, 214 10.1086/170270

Balbus S. A. , Terquem C. , 2001, ApJ, 552, 235 10.1086/320452

Barge P. , Sommeria J. , 1995, A&A, 295, L1

Barker A. J. , Latter H. N. , 2015, MNRAS, 450, 21 10.1093/mnras/stv640

Béthune W. , Lesur G. , Ferreira J. , 2017, A&A, 600, A75 10.1051/0004-6361/201630056

Birnstiel T. , Dullemond C. P. , Brauer F. , 2010, A&A, 513, A79 10.1051/0004-6361/200913731

Blaes O. M. , Balbus S. A. , 1994, ApJ, 421, 163 10.1086/173634

Burns K. J. , Vasil G. M. , Oishi J. S. , Lecoanet D. , Brown B. P. , 2020, Phys. Rev. Res., 2, 023068 10.1103/PhysRevResearch.2.023068

Chandrasekhar S. , 1961, Hydrodynamic and Hydromagnetic Stability. Clarendon Press, Oxford

Chen K. , Lin M.-K. , 2020, ApJ, 891, 132 10.3847/1538-4357/ab76ca

Chiang E. , Youdin A. N. , 2010, Annu. Rev. Earth Planet. Sci., 38, 493 10.1146/annurev-earth-040809-152513

Chiang E. I. , Goldreich P. , 1997, ApJ, 490, 368 10.1086/304869

Clarke C. J. , Pringle J. E. , 1988, MNRAS, 235, 365 10.1093/mnras/235.2.365

Cui C. , Bai X.-N. , 2020, ApJ, 891, 30 10.3847/1538-4357/ab7194

Cuzzi J. N. , Hogan R. C. , Shariff K. , 2008, ApJ, 687, 1432 10.1086/591239

Dubrulle B. , Morfill G. , Sterzik M. , 1995, Icarus, 114, 237 10.1006/icar.1995.1058

Flaherty K. et al. , 2020, ApJ, 895, 109 10.3847/1538-4357/ab8cc5

Flaherty K. M. , Hughes A. M. , Teague R. , Simon J. B. , Andrews S. M. , Wilner D. J. , 2018, ApJ, 856, 117 10.3847/1538-4357/aab615

Flaherty K. M. et al. , 2017, ApJ, 843, 150 10.3847/1538-4357/aa79f9

Flock M. , Nelson R. P. , Turner N. J. , Bertrang G. H. M. , Carrasco-González C. , Henning T. , Lyra W. , Teague R. , 2017, ApJ, 850, 131 10.3847/1538-4357/aa943f

Flock M. , Turner N. J. , Nelson R. P. , Lyra W. , Manger N. , Klahr H. , 2020, ApJ, 897, 155 10.3847/1538-4357/ab9641

Lizxandra F.-R. , Mario F. , Riouhei N. , 2020, A&A, 644, 8 10.1051/0004-6361/202039294

Gammie C. F. , Balbus S. A. , 1994, MNRAS, 270, 138 10.1093/mnras/270.1.138

Goldreich P. , Schubert G. , 1967, ApJ, 150, 571 10.1086/149360

Gressel O. , Ramsey J. P. , Brinch C. , Nelson R. P. , Turner N. J. , Bruderer S. , 2020, ApJ, 896, 126 10.3847/1538-4357/ab91b7

Gressel O. , Turner N. J. , Nelson R. P. , McNally C. P. , 2015, ApJ, 801, 84 10.1088/0004-637X/801/2/84

Johansen A. , Klahr H. , 2005, ApJ, 634, 1353 10.1086/497118

Johansen A. , Oishi J. S. , Mac Low M.-M. , Klahr H. , Henning T. , Youdin A. , 2007, Nature, 448, 1022 10.1038/nature06086

Johansen A. , Youdin A. , Mac Low M.-M. , 2009, ApJ, 704, L75 10.1088/0004-637X/704/2/L75

Klahr H. , Hubbard A. , 2014, ApJ, 788, 21 10.1088/0004-637X/788/1/21

Klahr H. H. , Bodenheimer P. , 2003, ApJ, 582, 869 10.1086/344743

Klahr H. H. , Henning T. , 1997, Icarus, 128, 213 10.1006/icar.1997.5720

Kunz M. W. , 2008, MNRAS, 385, 1494 10.1111/j.1365-2966.2008.12928.x

Kunz M. W. , Balbus S. A. , 2004, MNRAS, 348, 355 10.1111/j.1365-2966.2004.07383.x

Latter H. N. , 2016, MNRAS, 455, 2608 10.1093/mnras/stv2449

Latter H. N. , Papaloizou J. , 2018, MNRAS, 474, 3110(LP18)10.1093/mnras/stx3031

Lesur G. , Papaloizou J. C. B. , 2010, A&A, 513, A60 10.1051/0004-6361/200913594

Lesur G. R. J. , Latter H. , 2016, MNRAS, 462, 4549 10.1093/mnras/stw2172

Lin M.-K. , 2019, MNRAS, 485, 5221 10.1093/mnras/stz701

Lin M.-K. , Youdin A. N. , 2015, ApJ, 811, 17 10.1088/0004-637X/811/1/17

Lubow S. H. , Pringle J. E. , 1993, ApJ, 409, 360 10.1086/172669

Lyra W. , 2014, ApJ, 789, 77 10.1088/0004-637X/789/1/77

Lyra W. , Umurhan O. M. , 2019, PASP, 131, 072001 10.1088/1538-3873/aaf5ff

Malygin M. G. , Klahr H. , Semenov D. , Henning T. , Dullemond C. P. , 2017, A&A, 605, A30 10.1051/0004-6361/201629933

Manger N. , Klahr H. , 2018, MNRAS, 480, 2125 10.1093/mnras/sty1909

Manger N. , Klahr H. , Kley W. , Flock M. , 2020, MNRAS, 499, 1841

Marcus P. S. , Pei S. , Jiang C.-H. , Barranco J. A. , Hassanzadeh P. , Lecoanet D. , 2015, ApJ, 808, 87 10.1088/0004-637X/808/1/87

Marcus P. S. , Pei S. , Jiang C.-H. , Hassanzadeh P. , 2013, Phys. Rev. Lett., 111, 084501 10.1103/PhysRevLett.111.084501

Nelson R. P. , Gressel O. , Umurhan O. M. , 2013, MNRAS, 435, 2610(N13)10.1093/mnras/stt1475

Nelson R. P. , Papaloizou J. C. B. , 2004, MNRAS, 350, 849 10.1111/j.1365-2966.2004.07406.x

Ogilvie G. I. , 1997, MNRAS, 288, 63 10.1093/mnras/288.1.63

Ormel C. W. , Cuzzi J. N. , 2007, A&A, 466, 413 10.1051/0004-6361:20066899

Papaloizou J. C. B. , Nelson R. P. , Snellgrove M. D. , 2004, MNRAS, 350, 829 10.1111/j.1365-2966.2004.07566.x

Perez-Becker D. , Chiang E. , 2011, ApJ, 727, 2 10.1088/0004-637X/727/1/2

Pessah M. E. , Psaltis D. , 2005, ApJ, 628, 879 10.1086/430940

Petersen M. R. , Julien K. , Stewart G. R. , 2007a, ApJ, 658, 1236 10.1086/511513

Petersen M. R. , Stewart G. R. , Julien K. , 2007b, ApJ, 658, 1252 10.1086/511523

Pfeil T. , Klahr H. , 2019, ApJ, 871, 150 10.3847/1538-4357/aaf962

Pfeil T. , Klahr H. , 2020, preprint (arXiv:2008.11195)

Raettig N. , Klahr H. , Lyra W. , 2015, ApJ, 804, 35 10.1088/0004-637X/804/1/35

Richard S. , Nelson R. P. , Umurhan O. M. , 2016, MNRAS, 456, 3571 10.1093/mnras/stv2898

Sano T. , Miyama S. M. , 1999, ApJ, 515, 776 10.1086/307063

Schäfer U. , Johansen A. , Banerjee R. , 2020, A&A, 635, A190 10.1051/0004-6361/201937371

Shakura N. I. , Sunyaev R. A. , 2009, A&A, 500, 33 (special issue 500/01: reprint of 1973, A&A, 24, 337)

Simon J. B. , Bai X.-N. , Armitage P. J. , Stone J. M. , Beckwith K. , 2013b, ApJ, 775, 73 10.1088/0004-637X/775/1/73

Simon J. B. , Bai X.-N. , Stone J. M. , Armitage P. J. , Beckwith K. , 2013a, ApJ, 764, 66 10.1088/0004-637X/764/1/66

Stoll M. H. R. , Kley W. , 2014, A&A, 572, A77 10.1051/0004-6361/201424114

Stoll M. H. R. , Kley W. , 2016, A&A, 594, A57 10.1051/0004-6361/201527716

Tassoul J.-L. , 1978, Theory of Rotating Stars. Princeton Univ. Press, Princeton, NJ

Teague R. et al. , 2016, A&A, 592, A49 10.1051/0004-6361/201628550

Umurhan O. M. , Estrada P. R. , Cuzzi J. N. , 2020, ApJ, 895, 4 10.3847/1538-4357/ab899d

Umurhan O. M. , Nelson R. P. , Gressel O. , 2016a, A&A, 586, A33 10.1051/0004-6361/201526494

Umurhan O. M. , Shariff K. , Cuzzi J. N. , 2016b, ApJ, 830, 95 10.3847/0004-637X/830/2/95

Urpin V. , 2003, A&A, 404, 397 10.1051/0004-6361:20030513

Urpin V. , Brandenburg A. , 1998, MNRAS, 294, 399 10.1046/j.1365-8711.1998.01118.x

Wardle M. , 2007, Ap&SS, 311, 35 10.1007/s10509-007-9575-8

Weiss B. P. , Bai X.-N. , Fu R. R. , 2021, Sci. Adv., 7, eaba5967 10.1126/sciadv.aba5967

Youdin A. N. , Goodman J. , 2005, ApJ, 620, 459 10.1086/426895

Youdin A. N. , Lithwick Y. , 2007, Icarus, 192, 588 10.1016/j.icarus.2007.07.012

We explore the properties of the local dispersion relation equation ( 84) in the limit of |kRBR| ≪ |kZBZ|. When there is no vertical shear, R|$\partial$| Ω2/|$\partial$| 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.

Maximum growth rates and most unstable wavenumber are shown as functions of Ohmic Elsässer number Λ by equation ( A4) at fixed Alfvén velocity and R|$\partial$| Ω2/|$\partial$| Z = −0.1. Dashed lines are asymptotic solutions in the limit of η → 0 (equations A5 and A6) and η → ∞ (equations A7 and A8).

Oxford University Press is a department of the University of Oxford. It furthers the University's objective of excellence in research, scholarship, and education by publishing worldwide

This PDF is available to Subscribers Only

For full access to this pdf, sign in to an existing account, or purchase an annual subscription.

Free Design & Artwork Within 24 Hours

High Quality Is Our Obligation!

Have An Excellent Elite Sales Team

On the vertical shear instability in magnetized protoplanetary discs | Monthly Notices of the Royal Astronomical Society | Oxford Academic

Featured Products

News & Blog

China calls on rare earths companies to bring prices back to "reasonable" level | Reuters

Hydrogen fuels might have just got a huge leg-up | Top Gear

Vivid Kaya S12 loudspeaker |Stereophile.com

Vivid Kaya S12 loudspeaker |Stereophile.com

Rare-earth metals: the key to the renewable energy economy.

Consider the Mighty Bass Speaker - Premier Guitar

Build Your Own Speaker - Scientific American

Audio 101: How Do Headphones Work?

PM_Logo