{'fraction_non_alphanumeric': 0.03295063686764775, 'fraction_numerical': 0.02046457570950668, 'mean_word_length': 4.465876899304662, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Global Simulations of the Vertical Shear Instability with Non-ideal Magnetohydrodynamical Effects

The mechanisms of angular momentum transport and level of turbulence in protoplanetary disks (PPDs) are crucial for understanding many aspects of planet formation. In the recent years, it has been realized that the magneto-rotational instability (MRI) tends to be suppressed in PPDs due to non-ideal MHD effects, and the disk is largely laminar with accretion driven by magnetized disk winds. In parallel, several hydrodynamical mechanisms have been identified that likely also generate vigorous turbulence and drive disk accretion. We study the interplay between MHD winds in PPDs with the vertical shear instability (VSI), one of the most promising hydrodynamical mechanisms, through 2D global non-ideal MHD simulations with ambipolar diffusion and Ohmic resistivity. We find that for typical disk parameters, MHD winds can coexist with the VSI with accretion primarily wind-driven accompanied by vigorous VSI turbulence. The properties of the VSI remain similar to unmagnetized case, and the wind and overall field configuration are not strongly affected by VSI turbulence, showing modest level of variability and corrugation of midplane current sheet. Enhanced coupling between gas and magnetic field weakens the VSI. The VSI is also weakened with increasing magnetization, and we find that corrugation motions characteristic of the VSI transitions to low-amplitude breathing mode oscillations.


INTRODUCTION
Understanding the gas dynamics of protoplanetary disks (PPDs) is crucial for many aspects of planet formation. The key element to PPD gas dynamics lies in the mechanisms for disk angular momentum transport. This is firstly becuase it determines the overall disk structure and controls long-term disk evolution and secondly mechanisms responsible for angular momentum transport generally give rise to complex internal flow structures, particularly turbulence. Most stages of planet formation take place in PPDs, which sensitively depend on both of these factors (e.g., Armitage 2011). While angular momentum transport is not yet directly observable, theoretical studies are mainly guided by the fact that PPDs are actively accreting onto protostars at a typical rate of ∼ 10 −8 M yr −1 (e.g., Hartmann et al. 1998;Herczeg & Hillenbrand 2008), which places strong constraints on the possible mechanisms involved.
Angular momentum transport in PPDs has been conventionally attributed to the turbulence driven by the magneto-rotational instability (MRI, Balbus & Hawley 1991). However, this scenario is complicated by the fact that PPDs are extremely weakly ionized (e.g., Wardle 2007;Bai 2011a), which substantially weakens the coupling between gas and magnetic fields. This coupling is described by three non-ideal magnetohydrodynamic (MHD) effects: Ohmic resistivity, Hall effect and ambipolar diffusion (AD). They affect the MRI in different ways (e.g., Wardle 1999;Balbus & Terquem 2001;Blaes & Balbus 1994;Jin 1996;Kunz & Balbus 2004;Desch 2004), and particularly dissipations by resistivity and AD lead to damping or even complete quenching of the MRI (Gammie 1996;Bai & Stone 2011;Simon et al. 2013a,b). Via detailed simulations with increasing realistic disk microphysics, it has been shown that PPD accretion is likely primarily driven by magnetized disk winds Bai 2013;Gressel et al. 2015;Bai 2017) in a largely laminar disk: a major shift of paradigm.
Parallel to the development of the MHD theory, pure hydrodynamic mechanisms for driving angular momentum transport have also been discussed extensively in recent years. These include the vertical sheear instability (VSI; Nelson et al. 2013;hereafter N13), convective overstability in its linear (Klahr & Hubbard 2014;Lyra 2014;Latter 2016) and non-linear state (subcritical baroclinic instability; Klahr & Bodenheimer 2003;Petersen et al. 2007a,b;Lesur & Papaloizou 2010), and zombie vortex instability (Marcus et al. 2013(Marcus et al. , 2015Umurhan et al. 2016a). These instabilities generally operate under certain specific thermodynamic conditions, and can potentially be activated in different regions of PPDs (Fromang & Lesur 2017;Malygin et al. 2017;Pfeil & Klahr 2019). It is generally found in numerical simulations that the level of turbulence resulting from these hydrodynamic instabilities is not sufficiently strong to account for the typical disk accretion rates; they yield Shakura-Sunyaev (Shakura & Sunyaev 1973) α values up to ∼ 10 −4 −10 −3 , while α ∼ 10 −2 is likely required for the outer disk (e.g., Hartmann et al. 1998). But as potential sources of disk turbulence, they have important implications to planet formation. They can efficiently stir up dust particles (e.g. Stoll & Kley 2016;Flock et al. 2017a;Lin 2019), trap them in vortices and facilitate planetesimal formation (e.g. Raettig et al. 2015;Manger & Klahr 2018). Such turbulence are also needed so that sub-micron dusts are suspended to account for the disk SEDs (D'Alessio et al. 2001) and as seen in scattered light images.
So far, studies of wind-driven accretion generally adopted simplistic treatment of thermodynamics, partly to avoid the development of these hydrodynamic instabilities. On the other hand, studies of these hydrodynamic instabilities mostly ignored magnetic fields. A natural question arises: are hydrodynamic and MHD mechanisms compatible with each other? In other words, can these hydrodynamic instabilities operate in disks that launch magnetized disk winds?
In this paper, we focus on one of the most promising hydrodynamic instabilities, the VSI, and its interplay with magnetic fields. The VSI is an application of the Goldreic-Schubert-Fricke instability (Goldreich & Schubert 1967;Fricke 1968), originally derived in the context of differentially rotating stars, to accretion disks (Urpin & Brandenburg 1998;Urpin 2003;Arlt & Urpin 2004). Applications of the VSI to PPDs by N13 have drawn significant attention, leading to intense followup studies since then (e.g., Barker & Latter 2015;Umurhan et al. 2016a). Essentially, a disk that is Rayleigh-stable from radial shear becomes unstable in the presence of vertical shear where the rotational velocity varies over height, provided that cooling timescale is much shorter than orbital timescale to overcome the stabilizing effect of vertical buoyancy. In the limit of instant cooling, the criterion becomes (N13) where j is the specific angular momentum and is a function of cylindrical coordinates (R, Z), and k R , k Z are the radial and vertical wave numbers. The VSI taps free energy from the vertical shear in disk rotational velocities, which destablizes inertial-gravity waves (e.g., Lin & Youdin 2015). Such vertical shear is present when the fluid is baroclinic where contours of constant density and pressure are not aligned, and this is generally unavoidable in passively heated disks. The requirement for rapid cooling or thermal relaxation is more demanding. For standard models of PPDs, calculations suggest that the VSI can be triggered in the outer regions beyond ∼ 5 AU as well as the very inner regions (Lin & Youdin 2015;Malygin et al. 2017;Pfeil & Klahr 2019). Unstable VSI modes typically have large ratio of k R /k Z to take the advantage of the vertical shear, and non-linear simulations by N13 in idealized setup assuming vertically isothermal disks have found that the VSI develops into vigorous turbulence with prominent vertical oscillations. Cooling time well within ∼ 0.1 orbital time is needed to trigger the VSI, giving α values up to ∼ 10 −3 . More realistic simulations by Stoll & Kley (2014) incorporating radiative transfer show agreement with N13 results, though a reduced efficiency of angular momentum transport is reported. Recent 3D simulations further suggest the development of vortices (Richard et al. 2016;Manger & Klahr 2018), which is likely due to secondary Kelvin-Helmholtz instabilities (Latter & Papaloizou 2018).
In this paper, we conduct 2D global non-ideal MHD simulations of PPDs and study the possible development of the VSI, focusing on outer regions dominated by AD. While the VSI may operate in the innermost disk regions that are highly optically thick, the dynamics of these regions are much more complex and are likely to be dominated by the MRI turbulence due to the thermal ionization (e.g., Fromang et al. 2002;Desch & Turner 2015;Flock et al. 2017b). We restrict ourselves to cases where the disk is otherwise laminar (without VSI) to avoid complications from external driving sources. Note that even the system is laminar, there is still finite coupling between gas and magnetic fields, and hence a pure hydrodynamic understanding is incomplete. Recently, Latter & Papaloizou (2018) analyzed the linear properties of the VSI in the presence of magnetic fields. They found that in ideal MHD limit, while the MRI and the VSI modes lie in the same branch of the dispersion relation, MRI modes always dominate given its much higher growth rate. It is speculated that the VSI modes might survive when the MRI is suppressed, though they only mentioned Ohmic resistivity. We will directly examine whether the VSI can survive in the AD-dominated regime more relevant to realistic PPDs. If so, we will further address how the VSI is affected by magnetic fields, and in turn, whether the development of the VSI affects wind properties. As the first numerical simulations on the interplay between the VSI and the non-ideal MHD, we aim to make the problem as clean as possible. Therefore, we do not employ sophisticated implementations of realistic disk physics as in Bai (2017), but rather conduct numerical experiments similar to those in Bai & Stone (2017) with prescribed AD coefficients to mimic outer PPD conditions. We anticipate our results to serve as a benchmark for future studies that incorporate more realistic thermodynamics together with non-ideal MHD physics. This paper is organized as follows. In §2, we provide detailed descriptions of the numerical methods and the simulation setup. In §3, we list diagnostic quantities to furnish analyses of simulation results. Discussion on fiducial models of instability features are detailed in §4. The flow structure, mass accretion, and disk winds are discussed in §5. We conduct parameter study on magnetic field strengths, cooling timescales, and ambipolar diffusion strengths in §6. Finally, we summarize and discuss the main findings in §7.

METHODS AND SIMULATION SETUP
In this section, we describe the method and setup of our simulations, and in the meantime demonstrate the basic considerations that enter these simulations.

Dynamical equations
We use the grid-based high-order Godunov MHD code Athena++ (Stone et al., submitted) to conduct global simulations of the VSI in the context of PPDs. Athena++ is the successor of the Athena MHD code (Gardiner & Stone 2005, 2008Stone et al. 2008), but is rewritten in C++ which has much more flexible coordinate and grid options with significantly improved performance, scalability, and source code modularity. We solve standard MHD equations in spherical polar coordinates in conservative form, including non-ideal MHD effects as implemented in Bai & Stone (2017) where ρ, v, P are gas density, velocity, and pressure. The sum of the thermal and magnetic pressure is denoted by P * = P + B 2 /8π, where P * = P * I, with I being the identity tensor, and B = |B|. The gravitational potential of the central star with mass M * has the form Φ = −GM * /r. The total energy density is given by E = P/(γ − 1) + ρv 2 /2 + B 2 /8π, where v = |v|. The ideal gas law is adopted with adiabatic index γ = 7/5 for molecular gas in the bulk disk. 1 The term Λ c in the last equality represents the cooling rate to be detailed in §2.4. The electric field involves components from non-ideal MHD effects. In the local fluid rest frame, it reads where the Ohmic and ambipolar diffusivities are denoted by η O and η A , and we have ignored Hall term (see Section §2.3). The current density is J = c∇ × B/4π, and we express J ⊥ = −(J ×b) ×b as the component of J that is perpendicular to the magnetic field. The magnetic field B has its unit vector denoted byb = B/B. Poynting flux associated with non-ideal MHD is given by S = cE × B/4π. We use Gaussian units in the above equations, whereas in code units a factor of 4π is absorbed so that magnetic permeability is µ = 1. We perform global 2D simulations in spherical polar coordinates (r, θ, φ). We also use R = r sin θ and z = r cos θ to denote radial and vertical components in cylindrical coordinates. The unit system in the code has G = M = R 0 = 1, where R 0 is the reference radius fixed at the location of the inner boundary. The simulation domain covers from r = 1 to 100, and from θ = 0 to π (i.e., including the polar region) so that the domain contains sufficient dynamical range and is fully extended to accommodate wind launching. We use the van Leer time integrator, the HLLD Riemann solver, with piecewise linear reconstruction. Super time-stepping is used to accelerate the calculations of non-ideal MHD as in Bai & Stone (2017).

Disk model and initial conditions
As mentioned earlier, we aim to conduct relatively clean numerical experiments, hence our disk model is 1 The VSI primarily occurs in the disk where gas is molecular with γ = 7/5. By a transition above the disk surface, the gas is mainly atomic in the wind zone where γ = 5/3 (e.g., Wang et al. 2019). The exact value of γ affects the critical thermal relaxation timescale to trigger the VSI (Lin & Youdin 2015). Nevertheless, as we have experimented, we find effectively no difference in the outcome of the simulations between these two values.
set to be scale-free. In particular, disk temperature is set by the aspect ratio = H/r, where H = c s /Ω is the disk scale height, and c 2 s = P/ρ is the isothermal sound speed. Being scale-free (self-similar) requires to be independent of r, so that disk temperature is given by where T 0 = v 2 K0 = GM/R 0 . We consider disk temperature profile in θ to be In the above, the aspect ratio morphs from a disk value d to a surface value w . We employ a modestly thin disk with d = H d /r = 0.1 and set w = 0.5. The transition happens at θ trans = 3.5H d , with δθ the angle from θ trans . The transition width is controlled by n, where in our prescription n = 2, corresponding to a transition occurs within 0.05H d . Physically, the temperature transition is motivated by the external UV and X-rays heating from the protostar in the wind zone (e.g., Glassgold et al. 2004;Walsh etl. al. 2010).
We consider a density profile that is a power-law in radius where we set q D = 2. The function f (θ) is obtained through hydrostatic equilibrium solution. Solving for the Euler's equation in spherical polar coordinates in r and θ directions yields (Bai & Stone 2017) where we define F (θ) ≡ f (θ) 2 (θ). Note that f (θ) is an implicit function of (θ), hence numerical integration is applied to calculate F (θ) and then to compute f (θ). Radial and θ− velocities are set to zero for equilibrium, but we further add random noises to these velocity components with amplitude of ±5 percent of the local sound speed.
The poloidal magnetic fields are initialized by specifying an azimuthal vector potential (Zanni et al. 2007;Bai & Stone 2017), (12) Here, B z0 is the midplane field strength at R 0 . The parameter m determines the initial magnetic field configuration, quantifying the degree for which the magnetic field lines bend. A value approaches infinity m → ∞ gives a pure vertical field. In this paper, we choose m = 0.5. The poloidal field is computed by B = ∇ × A φ , so that one can arrive at B mid = B z0 (R/R 0 ) −(q D +1)/2ẑ for magnetic field at the midplane. Using vector potential guarantees that the resulting field is divergence free. Midplane poloidal field B is set so that the plasma β, defined as the ratio of gas to magnetic pressure, is constant. With midplane β 0 = 8πp 0 /B 2 mid , p 0 being midplane pressure, we choose β 0 = 10 4 in the fiducial model and further explore stronger and weaker fields.

Elsässer numbers
The strength of the non-ideal MHD effects are characterized by the dimensionless Elsässer numbers. For resistivity and ambipolar diffusion, they are given by where v A = B 2 /4πρ is the Alfvén speed. At a given ionization fraction, η O is independent of field strength or density, whereas η A ∝ B 2 /ρ 2 , and the Hall diffusivity η H ∝ B/ρ. Note that by definition, Am is generally independent of field strength. The Ohmic and Ambipolar Elsässer numbers are known to largely control operation of the MRI, with threshold about unity (e.g., Turner et al. 2007;Bai & Stone 2011). Here, we primarily consider the outer regions of PPDs where low densities make AD the dominant non-ideal MHD effect. While the Hall effect likely plays an important role at intermediate disk radii (e.g., Bai 2017), we do not include it to avoid complications, and that it would break the scale-free nature of our numerical experiments. We do include resistivity for purely numerical reasons (see below). In our fiducial simulations, we set Am = 0.3 for the bulk disk. Note that the value of Am is found to be on the order of unity towards PPD outer regions (Bai 2011a,b), and we choose it to be on the lower end which helps suppress the MRI. The value of Am then smoothly increases from the disk zone to the wind zone where Am is set to 100, attributed to the stellar irradiation of FUV and X-rays that substantially elevate the ionization level above the disk surface (Perez-Becker & Chiang 2011), recovering the ideal MHD regime. The functional form of the transition is similar to that of temperature, since heating and ionization in the disk atmosphere/wind zone are both due to UV/X-rays.
For simulations with pure AD, the system unavoidably develops a current sheet in the midplane region where B φ flips sign (e.g., Bai & Stone 2017;Bai 2017;Suriano et al. 2018). This current sheet tends to be unstable, which then corrugates and develops into more complex structures (see Appendix B for more discussion). While this is interesting on its own right and will be investigated in a future publication, it imposes difficulties to assess the development and saturation of the VSI and causes problems near the inner radial boundary. Given the experimental nature of our simulations, we thus apply resistivity near the inner boundary and a thin layer (∼ ±H d ) in the midplane region throughout to stabilize the system (see Figures 5). It allows us to investigate the development of the VSI on top of a laminar disk that launches MHD disk winds. The value of midplane resistivity is fixed to η O = 0.05c s H d , and it rapidly declines towards the surface being negligible within 0.1H d . The fiducial run with β 0 = 10 4 has an initial midplane Ohmic Elsässer number of Λ = 4 × 10 −3 .

Thermal relaxation
Cooling of the system is achieved through thermal relaxation, which is associated with the Λ c term on the right hand side of Equation (5). We adjust disk temperature as where T is the temperature at time t, and the target temperature T eq is set to the equilibrium temperature by initial condition. A relaxation time τ is prescribed to be a fraction of the local Keplerian orbital period, i.e. τ (R) ∝ P orb = 2π/Ω K (R). More precisely, we adjust disk temperature at each timestep by where ∆T = T (t + ∆t) − T (t). The terminology utilized for thermodynamics throughout this paper is i) locally isothermal, which for ∆t τ the disk temperature is kept to its equilibrium value at each position; ii) thermal relaxation, for which disk temperature is recovered to its equilibrium value on some timescale determined by τ . Following Bai (2017) and Bai & Stone (2017), the wind zone is set to be locally isothermal for all runs with a smooth transition from disk zone. 2

Boundary conditions
The conditions of inner radial boundary deserve special attention. We fix the hydrodynamic variables to initial equilibrium state, with temperature and density computed from Equations (7) and (9). The angular velocity is set to the minimum between initial v φ through Equation (11) and Ω K (r 0 )R, along with the other two velocity components v r = v θ = 0. Such a fixed state boundary condition can provide more stable flow structure and minimize the influence from inner boundary to the main simulation domain. Furthermore, we add a buffer zone between r = r 0 and r = 1.5r 0 by including a constant resistivity from the midplane all the way up to near the pole with values at each θ equal to midplane resistivity. It helps stabilize the magnetic flux originating from the inner boundary though not completely over long times.
The rest of the boundary conditions are standard. The hydrodynamic variables in the outer radial boundary are extrapolated from the last grid zone, assuming ρ ∝ r −q D , T ∝ r −q T , v φ ∝ r −1/2 . Radial and meridional velocities v r and v θ are copied directly from the last grid zone, except setting v r = 0 when v r < 0. Magnetic field variables in the inner and outer ghost zones are determined via B r ∝ r −2 , B θ ∝ const., B φ ∝ r −1 . With θ domain reaching to the pole, we employ polar boundary conditions, where θ and φ components are reversed across the polar boundary.

Simulation runs
Our simulation runs are listed in Table 1, with three main physical parameters: the disk magnetization β 0 , the thermal relaxation timescale τ , and the disk ambipolar diffusion Elsässer number Am. For the fiducial run, we choose β 0 = 10 4 , τ = 0 and Am = 0.3.

TABLE 1 List of Simulation Models and Parameters
Note. -FidH and Fid are hydrodynamic and MHD fiducial simulations, respectively. P orb and P 0 denote local orbital period and innermost orbital period. Run FidH has runtime of 600P 0 and B3 has runtime of 800P 0 , and the rest models have runtime of 500P 0 .
The simulation models listed in Table 1 are all carried out with the same resolution, 1536×480 grid cells in (r, θ), except model ResH. The radial domain spans from r = 1 to 100 and has logarithmic spacing. Note that our radial domain is much wider than previous hydrodynamic simulations of the VSI. This is because to properly accommodate the MHD winds, domain size should generally be much larger than wind launching radius. Grid spacing in θ direction increases by a constant factor of 1.006 per grid cell from the midplane toward the two poles so that the resolution at the midplane is four times finer than that at the pole. This enables us to achieve a resolution of about 32 cells per H d in r and 32 cells per H d in θ at the disk midplane. Run ResH is a hydrodynamic run with limited θ domain and no temperature transition for comparison with earlier works in literature. Its limited spatial domain spans r ∈ [1, 4] and θ ∈ [−5H d /r, 5H d /r]. The resolution achieves 96 cells per H d in r and 108 cells per H d in θ with uniform grid spacing. In Appendix A, we demonstrate our adopted resolution is sufficient for numerical convergence.
For comparison, pure hydrodynamic simulations with the identical setup as in MHD runs are also conducted by the same set of thermal relaxation timescales. Note that these hydrodynamic simulations differ from existing simulations (e.g., those in N13) in a way that we cover much more extended θ−domain with a temperature transition. We will see that most VSI activity is bounded by this temperature transition region even the disk atmosphere is set to have instantaneous cooling.

DIAGNOSTICS
We list relevant diagnostic quantities in this section to facilitate the analyses of simulation results.

Kinetic energy
One quantity of great interest is the perturbed energy, which is defined as the sum of volume-integrated meridional and radial kinetic energies nomarlized by the initial equilibrium angular velocity v φ0 (N13): We compute the kinetic energies in a box of r ∈ [2, 4] and θ ∈ [−3H d /r, 3H d /r]. The evolution of kinetic energies will be frequently invoked in result sections to facilitate comparisons of different simulation models.

Stress tensors
In steady state, the equation of angular momentum conservation in cylindrical coordinates readṡ where the angle brackets stand for temporal and spatial average. The net mass accretion rate on the left hand side of Equation (17) is given bẏ and two terms on the right associate to transport due to the radial and vertical components of the angular momentum flux, corresponding to the Rφ and zφ components of the stress tensor, are defined as (Balbus & Hawley 1998) The hydrodynamic and magnetic parts of the stress tensors are known as Reynolds (T Rey ) and Maxwell (T Max ) stress, respectively. For vertical transport, we retain only the Maxwell stress, because the Reynolds component, ρv z v φ does not carry excess angular momentum from the disk. The vertical height of the wind base is marked by |z wb | = 3.5H d which separates the disk and atmosphere. Empirically, it is the location where non-ideal MHD dominated disk zone transitions to ideal MHD dominated wind zone. For better statistics, we do it over a time interval from 200 to 300P 0 as well as over a cylindrical radial extent R ± 0.1 around the radius of interest. Radial transport is generally mediated by turbulence, and is characterized by the the classic dimensionless α parameter (Shakura & Sunyaev 1973), and at each location it is defined as The role of radial (viscous) transport on disk accretion is reflected in the vertically integrated α, defined as and by assuming steady-state viscously driven accretion, it yields an accretion rate given bẏ where Σ = z wb −z wb ρ dz is the surface density. In the second relation, we replace ∂/∂R by R −1 for an order of magnitude estimate.
Vertical transport of angular momentum is mediated by magnetized disk winds. Because of long lever arm, it is more efficient than radial transport by a factor ∼ R/H, given B R and B z are on the same order of magnitude. Inside the disk, wind-driven accretion is achieved by a torque exerted by the Lorentz force, and in the thin disk limit, B z is approximately constant, the local accretion velocity given by (Wardle 2007 which states that accreting velocity is proportional to the vertical gradient of B φ .

Wind kinematics
The local mass loss rate per logarithmic radius of wind is formulated by whereṀ wind (R) is the cumulative wind mass loss rate within radius R. An important diagnostic of the magnetized winds is the Alfvén radius, which is the point along a filed line where the poloidal gas velocity v p equals the local Alfvén velocity v Ap = B p / √ 4πρ. The wind mass loss rate is related to the mass accretion rate in the bulk disk by (Ferreira & Pelletier 1995), where ξ is called the ejection index, and the ratio R A /R wb is often referred to the magnetic lever arm. Thus, the location of the Alfvén point can provide a convenient measure of the wind mass loss rate.
3.4. VSI linear modes General linear analysis on pure hydrodynamic VSI can be found in N13. With locally isothermal equation of state, and in the short wavelength limit, the VSI growth rate σ is given by where Ω 0 is the angular velocity of disk at center of the shearing box, and κ 0 denotes epicyclic frequency. This expression directly relates the vertical shear ∂v φ /∂Z to the growth rate of the wave modes. The fastest growing mode then satisfies with growth rate Local linear analysis incorporating magnetism has been conducted by Latter & Papaloizou (2018) under the assumption of ideal MHD. For the VSI modes to fit into the disk, its wavelength should not exceed disk scale height, and this gives an upper limit to field strength where β represents the ratio of thermal pressure to magnetic pressure, and q is a measure of the vertical shear strength. Taking q ∼ H/R = 0.1, we have β 100 (Barker & Latter 2015). Equations (27) and (30) will later be involved to discuss why VSI weakens with increasing disk magnetization ( §6.1).
The above mentioned linear analyses are performed under the local shearing sheet assumption in short wavelength limit. Vertically global analyses demonstrate the emergence of two types of modes, namely the body modes near the midplane which can be further categorized into breathing and corrugation modes, and the high latitude surface modes when applying no-flow boundary conditions (N13; Barker & Latter 2015;McNally & Pessah 2015;Umurhan et al. 2016b). The non-linear evolution of the VSI is eventually taken over by large-amplitude corrugation modes.

FIDUCIAL MODELS
Starting with an analysis on fiducial hydrodynamic model FidH, we demonstrate its consistency to the limited θ-domain simulation ResH ( §4.1). Then we present results from MHD simulations by discussing fiducial model Fid in detail to illustrate the main features of the VSI when incorporating magnetic fields ( §4.2, §4.3, and §4.4). Analyses on model B3 and B5 are detailed in §6.1.
In presenting the results, we use time in units of the innermost orbital time P 0 = P orb (R 0 ) = 2π. If not otherwise noted, we analyze the data at radii between R = 2 and R = 4, and average the results from time 200 to 300P 0 for model Fid and B5, and 600 to 700P 0 for model B3. These time intervals are chosen after the VSI has fully saturated (see Figures 2 and 12).

The hydrodynamic runs
In Figure 1, we show the meridional velocities at time t = 150P 0 for run FidH and ResH. Note the difference between the two runs is that run FidH has an extended domain with a temperature transition at z = ±3.5H d . As discussed in N13, there are higher latitude surface modes and lower latitude body modes in the early stage of their simulations, and the body (corrugation) modes dominate the non-linear state of evolution. Surface modes generally require high resolution to be resolved due to their short radial wavelengths. Our run ResH achieves comparable resolution as in N13, and all these modes are observed where breathing modes show up when initial perturbations are small (∼ 10 −6 c s ). Surface and breathing modes are ultimately taken over by corrugation modes in the nonlinear stage, showing prominent vertical oscillations. While run FidH has lower resolution and does not show surface modes, we see that the corrugation modes at the non-linear state are well captured.
As seen in Figure 1, the two runs share very similar turbulence patterns, except that in run FidH, the VSI is well confined in the bulk disk within the transition region. Turbulence properties in between the two runs also converge, where we find turbulent kinetic energy as defined in Equation (16)  ∼ 40P 0 for both models (see Figure 2). The Reynolds stresses between the two models are also consistent with each other, on the order of α Rφ ∼ 10 −3 .
N13 found turbulent kinetic energy in their simulations rises unphysically after a few hundred orbits, which was attributed to their restricted computational domain and θ boundary conditions. We have also observed this phenomenon in run ResH. While our run FidH with much more extended domain shows fluctuations in kinetic energy after about 300P 0 , the mean value remains similar, and the turbulence structure sustains in the long run.

Time evolution in MHD simulations
The time evolution of our fiducial MHD run Fid is illustrated in four figures. Figure 2 shows the time evolution of kinetic energy fluctuations defined in Equation (16). In Figures 3 and 4, snapshots of vertical velocity v z and magnetic fields are displayed along the evolutionary sequence. In these two Figures, we also show results of the simulation with τ = 1 where the VSI does not develop for comparison. In Figure 12, we provide a space-time plot of vertical velocity at R = 3 to show the development of VSI consecutively.
With poloidal fields, initial stage of the evolution is characterized by wind launching from the surface layer, which occurs around the transition region where gas becomes better coupled to the magnetic field. This process generates some disturbances that traverse the simulation domain in a few tens of orbits that is visible in Figure  2. Once the wind is established in a few tens of P 0 , the overall field configuration remains quasi-steady state as seen in Figure 4. Note that when setting the cooling time to be τ = 1 (last panels of Figure 3 and 4), the wind configurations remain stable and symmetric throughout the simulation time. Therefore, turbulent fluctuations we observe at later times in run Fid should arise from the VSI.
We observe the VSI being developed in the bulk disk body progressively from small to large radii, on top of the magnetized wind. In Figure 3, we see that at t = 25P 0 , wave patterns in a region R < 1.5 are associated with VSI wave modes. Features at R > 1.5 correspond to disturbances from initial relaxation. We mainly focus on regions between R = 2 to 4, for which it takes ∼ 100P 0 for the VSI turbulence to cover this entire region.
The overall development of the VSI in run Fid is similar to that in the hydrodynamic counterpart run FidH. In Figure 2, we see the VSI of run Fid takes longer to develop, which is a general trend in magnetized runs (see §6.1). In Figure 3, it can be seen that regions around R = 2 show breathing modes (body modes with odd-symmetry) at t = 50P 0 , and transitions at t = 60 − 70P 0 to become corrugation modes with evensymmetry. These features are all very similar to the pure hydrodynamic case, though the transition is slightly delayed when comparing the results at the same radial range. Due to the lower resolution, surface modes found in the initial development of the VSI in hydrodynamic simulations of N13 are not observed, though in N13 it was noted that these modes are later overtaken by corrugation modes. As can be seen in Figure 1, finer turbulent structures are present in higher resolution runs. It is unclear whether such finer turbulent structures affect the diffusion magnetic fields, though we expect them to be negligible compared to the much stronger non-ideal MHD effects. By t = 200P 0 and within the entire region of R = 4, the instability is saturated with vigorous VSI turbulence, characterized by strong vertical oscillations. The radial spacing of these oscillation columns is about ∼ 2−3H d , again similar to the pure hydrodynamic case. Oscillation period can be also obtained in Figure  12, which is about 20P 0 at R = 3 corresponding to 3.8 local orbits. The development of the VSI modifies the magnetic field configuration. Around the transition zone, radial magnetic field appears to show some bunching behavior as seen in upper panels of Figure 4. This results from radial shear in vertical oscillations of the VSI. In the disk zone, the τ = 1 simulation shows the standard pattern with smooth toroidal field that flips at the midplane, leading to a current sheet (e.g., Bai & Stone 2017). In the bottom panels of Figure 4, the current sheet corrugates and oscillates around the midplane in accordance with the VSI oscillations. As a consequence, the toroidal field component becomes less smooth showing some banded patterns.
One caveat in our simulations is that the magnetized wind configuration is only quasi-steady. As found in Bai & Stone (2017), disks tend to gradually lose magnetic flux over time. At numerical level, the inner boundary is gradually depleted of magnetic flux and becomes less magnetized, which causes the region to be less stable. This can be seen in the sixth and seventh panels in Figures 4 and 3, and the effect starts to affect the region of interest after about 300P 0 , and is related to the succes- sive decline of the kinetic energy beyond ∼ 300P 0 seen in Figure 2 and the corrugation oscillation period also becomes more irregular as seen in Figure 12. Therefore, in the rest of the paper we restrict our analyses for model Fid in a time interval from 200 to 300P 0 and mainly focus at radius R = 3.

Diagnostic vertical profiles
To examine the gas dynamics with the presence of the VSI, we show in Figure 5 the vertical profiles of major diagnostic physical quantities in a vertical extent z ∈ [−5H d /r, 5H d /r], which covers the entire disk zone and part of the wind zone, at a cylindrical radius R = 3. For comparison, the first and second rows show time averaged results from run Fid, and the corresponding profiles of the model with thermal relaxation time τ = 1.
Dashed curves correspond to profiles measured in initial hydrostatic equilibrium (t = 0).
All profiles are either symmetric or anti-symmetric about the midpalne. In the second column, we see that toroidal fields are the dominant field component which reaches a factor ∼ 5 of midplane field and undergoes a flip of sign in the midplane. As seen in the fourth column, the flip is achieved over a relatively thick layer of z = ±H due to the midplane resistivity added to stabilize the current sheet. Since the transition is at z ∼ ±3.5H d , the bulk disk is still AD dominated that is sufficient for us to study its interplay with the VSI for the outer disk. The relatively low transition height which may be expected for outer disk (e.g., Perez-Becker & Chiang 2011) leads to the fact that the gas is still relatively dense there, and the entire disk is gas pressure dominated. Nevertheless, overall wind profiles remain similar to those obtained in previous wind simulations (e.g., Bai & Stone 2017). The poloidal velocity is accelerated along field lines reaching a substantial fraction of midplane Keplerian speed, whereas rotation remains sub-Keplerian at all heights as seen in the third column due to the pressure gradient. A vertical shear in the rotational velocity is clearly present which is the driving source of the VSI. The wind profiles between run FidH which is turbulent but time averaged and its τ = 1 counterpart which is laminar without VSI closely resemble each other. This is a clear indication that the VSI does not affect the general properties of the disk wind.

Turbulence strengths
To assess the level of VSI turbulence, we calculate the root mean square of velocity fluctuations as δv ≡ (v − v ) rms , where angle brackets denote the temporal and spatial averaging. We average the quantities over a time interval from 200 to 300P 0 and over a spatial extent R = 3 ± 0.1. The results are shown in Figure 6 for runs B3, Fid, and B5.
We focus on run Fid here in the middle panel. The perturbed velocities in τ = 0 run is around 10% local sound speed in the bulk disk, and is dominated by the vertical component as is expected in the VSI. The radial and azimuthal component of the fluctuation rises towards the wind zone. At the wind base around the transition height (z ∼ ±3.5H d ), total velocity fluctuations is at about 20 − 30% midplane sound speed, which leads to fluctuations in wind kinematics (see more discussions in §5.2.3). By contrast, velocity fluctuations are at much lower level by several orders of magnitude in the τ = 1 case.

ANGULAR MOMENTUM TRANSPORT
In this section, we analyze angular momentum transport processes in presence of both the VSI and disk winds through the radial ( §5.1) and vertical ( §5.2) transport.

Radial Transport
In Figure 7, we show components of α Rφ , the time and spatially averaged Reynolds and Maxwell stress, normalized by the midplane gas pressure at R = 3 ± 0.1. Focusing on the middle panel for our fiducial model Fid, the Reynolds stress due to the VSI is appreciable throughout the disk column achieving a maximum of ∼ 5 × 10 −4 at z ∼ ±1.5H d . The τ = 1 model shows a negligible Reynolds stress as expected.
The Maxwell stresses of both the locally isothermal model and thermal relaxation (τ = 1) model drop near the disk midpane (z ∼ ±H d , Figure 5). While this is mainly due to (artificially) enhanced Ohmic resistivity applied to this region, we note that in reality the flip of B r and B φ should give rise to a similar but much narrower feature. The Maxwell stress rises towards disk surface as accompanied by wind launching, and reaches a local α Rφ ∼ 10 −3 near the wind base. With the VSI, additional turbulent fluctuations increase the Maxwell stress in the midplane region so that the overall vertical profile is relatively flat compared to the τ = 1 case.
We note that the Reynolds stresses measured in all three of our magnetized models are noisy and are smaller than those obtained in 3D simulations (α ∼ 10 −3 ; e.g. N13; Manger & Klahr 2018) which allow for better averaging. We find the same holds for hydrodynamic simulations. The higher Reynolds stress in 3D is likely due to the formation of vortices (Richard et al. 2016;Manger & Klahr 2018), though the kinetic energy fluctuations remain very similar between 2D and 3D. Further study with non-ideal MHD in 3D is needed to better constrain the range of Reynolds stress.

Vertical Transport
Without VSI turbulence, previous works (e.g., Bai 2017) show that magnetized disk winds are predominant in driving disk accretion. In this section, we examine how the VSI affects the winds properties.

Flow Structure
As explained earlier, for purely wind-driven accretion, the vertical profile of accretion (radial) velocity should be given by Equation (24). This is clearly seen in the right panel of Figure 8 for simulations with τ = 1, similar to the results in Bai & Stone (2017). The accretion layer is confined to within z = ±2H d where B φ flips, and the mass flux agrees very well with expectations. Gas starts to flow radially outward as the gradient of B φ reverses. Some deviation in the upper layer can be attributed to the fact that the B z can no longer be taken to be constant as assumed.
In the left panel of Figure 8, the onset of the VSI modifies flow structure in both accretion and decretion layers. The accreting/decreting layer becomes narrower/broader, and the mass fluxes of these two layers are both enhanced. The bulk disk is now not in steady state, hence Equation (24) is inapplicable to estimate the flow structure. Sudden drops of the dotted black curve near z = ±2.5H d are associated with the bunching of poloidal fields, and VSI turbulence tends to restrict the accreting flow toward the midplane within z = ±H d . Figure 9 shows the radial profiles of accretion and mass loss rates, computed at R = 3 and averaged over t = 200 − 300P 0 . The right panel shows accretion rate without the VSI. Clearly, accretion is mainly driven by wind, and it dominates the contribution from Rφ component of the Maxwell stress by a factor of ∼ R/H, as expected. There is excessive mass loss, with mass loss rate a factor of several higher than accretion rate. Correspondingly, we find that the Alfvén radius is very close to the wind base, as expected from Equation (26). This is mainly controlled by the location of wind base, and higher wind base would lower the mass loss rate (e.g., Bai et al. 2016). Detailed understanding of wind dynamics requires realistic calculations of heating/cooling process in the disk atmosphere (e.g., Wang et al. 2019) and is beyond the scope of this work. Here, we are mainly concerned the gas dynamics in the bulk disk and the interplay between VSI and winds, and we have experimented to confirm that our overall conclusions are robust when setting the wind base at different heights.

Mass accretion rate
The left panel of Figure 9 shows accretion rates in the presence of the VSI. The bulk accretion and mass loss rates are similar to the τ = 1 case, indicating that the VSI plays a minor role in the overall processes including angular momentum transport and mass loss. Still, magnetized wind remains the dominant mechanism in driving disk accretion. The main differences due to the VSI is that accretion driven by the Rφ component of Reynolds stress is significantly boosted, reaching a level comparable to that from the Maxwell stress. We also note that with the VSI the profiles are noisier even after some timeaveraging. To better quantify the role of radial transport, we use Equation (23) for an order-of-magnitude estimate, taking α ∼ 10 −3 , which is slightly larger than the values obtained from our fiducial run (sum of Maxwell and Reynolds stress), but comparable to those measured in 3D hydrodynamic simulations ( §5.1). Converting to code units, we find a value of ∼ 1.8 × 10 −5 at R = 3, which is again about a factor H/R to wind-driven accretion rates.

Wind variability
The turbulent nature of the bulk disk due to the VSI induces time variabilities in the wind. To show this, we trace a poloidal field line from midplane at R 0 = 3 all the way up to the boundary of the simulation domain. In Figure 10, we present the time evolution of poloidal velocity and mass flux along this field line at cylindrical radii of R = 3.5 and R = 15 over a time interval 200 − 400P 0 . The poloidal velocities measured at the two radii increases with distances, implying the wind is kept accelerating. The wind from the τ = 1 case (no VSI) is fairly steady, as the bulk disk is in quasi-steady-state. With the VSI (τ = 1), wind velocity and wind mass flux are on average about the same as the the τ = 1 case. In the meantime, wind velocity fluctuates at about ±25% at R = 3.5 and ±8.5% at R = 15, whereas the mass flux fluctuates around ±34% at R = 3.5 and ±26% at R = 15. These results suggest that the VSI affects wind kinematics only at a modest level.
6. PARAMETER STUDY In this section, we investigate the effect of magnetic field strength ( §6.1), the thermal relaxation timescale ( §6.2), and the ambipolar diffusion strength ( §6.3) on the onset of the VSI.

Magnetic field strength
We first look at the evolution of kinetic energies for all runs in Figure 2. It is notable that the development of the instability is closely related to the magnetic field strength; stronger fields result in a time delay on the onset of the VSI and weaker kinetic energy fluctuations. Specifically, initial growth proceeds the fastest in the hydrodynamic run FidH. Growths in models B5 and Fid lag a few tens of orbital periods behind. These runs saturate at similar turbulence levels. The growth of model B3 is much more prolonged, and is only noticeable after ∼ 300P 0 and saturates after about 450P 0 . Overall, run B5 resembles the behavior as run Fid except that the body (corrugation) modes show earlier growth and saturate to slightly higher amplitude (Figure 2). Run B3 is more special and we discuss in more detail below.
The vertical velocity fluctuations at different stages of the disk evolution for run B3 are shown in Figure 11. At t = 200P 0 , we can identify the presence of breathing modes with odd symmetry. The amplitude of such fluctuations saturates after about 500P 0 at radii of R = 2−4 and remains till the end of the simulation at 800P 0 (Figure 2). We find that the system is entirely dominated by such low-amplitude breathing mode (top panel of Figure 12). The oscillation period at R = 3 is around ∼ 50P 0 ∼ 10 local orbits, and is not entirely steady as the radial spacing of such mode also varies with time.
Because of the symmetry of the breathing mode, the current sheet remains largely unperturbed in the equatorial plane. We also show results from the simulation with τ = 1 in the last panel of Figure 11. By contrast, there is very little fluctuation in vertical velocity and no sign of such breathing mode, which further supports that what we observe in run B3 is a consequence of the VSI.  . We see that the velocity fluctuations of run B5 generally resemble the behavior of run Fid, which are dominated by vertical motions similar to the pure hydrodynamic case at quantitative level. Velocity fluctuations in run B3 is much weaker. Moreover, as the corrugation mode is suppressed, vertical velocity fluctuation appears to be the weakest, whereas there are stronger velocity fluctuations in v φ . Figure 7 is similar to Figure 6, but for vertical profiles of the Rφ component of Reynolds and Maxwell stresses. For Maxwell stress, the trend is clear that stronger magnetization gives higher Maxwell stress which accompanies wind launching as usual. 3 On the other hand, run B5 yields stronger Reynolds stress than the fiducial run, while run B3 gives much less. Since the Reynolds stress is mainly caused by the VSI, this anti-correlation with field strength again reflects the VSI weakens with stronger magnetization.
The fact that stronger magnetization leads to weaker VSI may have two possible causes. First, wind launching modifies the vertical shear profile and thus likely alters the free energy source that feeds the VSI. This can be seen in Equation (29), where linear growth rate directly scales with this shear rate, though it does not apply to the vertically global body modes. In Figure 13, we show the vertical shear profiles of runs FidH, B5, Fid, B3 in τ = 1 (no VSI) and τ = 0 (with VSI) cases. Comparing the two panels, we see that the VSI modifies the shear profile at modest level. Interestingly, some narrow feature is observed at ∼ 3H d in run B3 that is related to wind launching, in both τ = 1 and τ = 0 cases, whereas below this region, the shear profile remains similar to runs with weaker magnetization. Given the complexity of this feature, however, it is difficult to deduce whether it is responsible for the transition to breathing modes or not.
The perspective above employs the results of linear analysis from pure hydrodynamics on top of a magnetically-modified vertical shear profile, thus ignores the coupling between gas and magnetic field. In the opposite limit, linear analysis by Latter & Papaloizou (2018) mainly assumed ideal MHD limit (i.e., perfect coupling between gas and magnetic field), and found that the MRI generally overwhelms the VSI, and that the VSI modes are stabilized by magnetic tension when β q −2 ∼ 100 (Equation 30). We find that magnetization in our run B3 has reached this limit for z ≥ 2H d , which is consistent with the lower-than-expected growth rate discussed earlier. We will further discuss in §6.3 on the role played by ambipolar diffusion. We also note that the linear analysis in Latter & Papaloizou (2018) was conducted in the short-wavelength limit, applicable mainly to the surface modes, whereas the more relevant is the body modes requiring vertically-global treatment. While this has been done in the hydrodynamic framework (N13; Barker & Latter 2015;McNally & Pessah 2015), it can be highly challenging to further incorporate vertical magnetic fields because of the lack of a welldefined initial equilibrium.

Thermal relaxation timescale
The thermal relaxation timescale is decisive to the onset of the VSI. Short thermal timescales can diminish the stabilizing effect from buoyancy, hence furnish the growth of the instability. Here, we consider the evolution of models by relaxing the locally isothermal assumption. The thermal timescales investigated are fractions of the local orbital periods. We chose τ = 0.1, 0.01, 0.001 corresponding to model t-1, t-2 and t-3, along with the fiducial model Fid where τ = 0. For comparison, each MHD model is accompanied by a hydrodynamic run with the same thermal relaxation prescription.
In Figure 14, we show the evolution of the kinetic energy fluctuations for models with various thermal timescales. In pure hydrodynamic simulations, it is clear that in the τ = 0.1 model, the kinetic energy fluctuations quickly damps and then maintain some constant and very low level, indicating no VSI growth. Model τ = 0.001 shares great similarities with locally isothermal run FidH, whereas Model τ = 0.01 grows slower and saturates at a turbulence level with kinetic energy about twice lower. The critical cooling time of τ = 0.1 is consistent with previous analytical studies and hydrodynamic simulations with either simplified disk models or more realistic radiative transfer (e.g. N13; Stoll & Kley 2014;Lin & Youdin 2015).
The solid curves in Figure 14 show the corresponding MHD runs. Comparing the results with the pure hydrodynamic runs, we see that the threshold thermal relaxation timescale is similar between the two cases. For τ = 10 −3 −10 −2 , the level of saturation is also similar between hydrodynamic and MHD runs, where the system shows prominent vertical oscillations which are characteristic of the corrugation modes.

Ambipolar diffusion strength
The role played by non-ideal MHD effects (especially ambipolar diffusion) in our simulations can be understood from two aspects. First, they suppress the MRI, thus allowing the VSI to stand out. Second, they break the flux freezing condition to make the gas behave close to the unmagnetized case. Latter & Papaloizou (2018) qualitatively discussed the VSI under the Ohmic resistivity, where the VSI should be able to operate completely unimpeded when magnetic diffusion is able to weaken magnetic tension over the timescale of VSI growth, which translates to an Ohmic Elsässer number Λ < q, where q is defined in Equation (30). Ambipolar diffusion shares some similarities with Ohmic resistivity in its dissipative nature, though a more detailed examination is needed. Here, we show the influence of Am from our numerical results below.
We explore the evolution of perturbed kinetic energies as a function of ambipolar diffusion Elsässer number Am which is prescribed as a constant within the bulk disk ( Figure 5). Realistic values of Am is found to be on the order of unity towards outer regions of PPDs (e.g., Bai 2011a). Here, we further increase the value of Am to 0.5 and 1, and the time evolution of kinetic energy fluctuation of these runs are shown in the right panel of Figure 14. We see that the VSI becomes weaker when Am is larger, which corresponds to smaller diffusivity. This result supports the expectation discussed earlier that strong magnetic dissipation can assist the development of the VSI.
From these results, we see that in the limit of strong ambipolar diffusion (Am 1), the disk can be considered as unmagnetized and we expect vigorous VSI turbulence if thermal relaxation in the disk is sufficiently rapid. In the opposite limit (Am 1), the MRI will take over the VSI that again leads to strong turbulence. Hence, weakest turbulence is likely to occur in between these two limits.

CONCLUSIONS AND DISCUSSION
In this paper, we study the onset and non-linear evolution of the vertical shear instability, which has so far been studied mostly in the hydrodynamic framework, in the presence of MHD disk winds. We perform 2D simulations in spherical polar coordinates in the r − θ plane using Athena++ code. We focus on outer regions of PPDs, where thermodynamic conditions favor the development of the VSI. Non-ideal MHD effects including ambipolar diffusion and Ohmic resistivity are incorporated, with prescribed diffusivities that suppress the MRI in the bulk disk and launch MHD winds from the surface. The main results are summarized as follows.
• The VSI generates vigorous turbulence in the presence of magnetized PPD disk winds. The turbulence properties are similar to pure hydrodynamic case, yielding velocity fluctuations on the order of 0.1c s dominated by vertical motions. The Shakura-Sunyaev α parameter is a few times 10 −4 in our 2D simulations.
• Magnetized disk winds persist despite of the VSI turbulence. The bulk wind properties remain similar to the cases without the VSI, with small variability at the level of ∼ 10 − 30%. The midplane current sheet becomes corrugated as a result of the VSI, and poloidal magnetic field shows bunching in the surface layer.
• The growth rate of the VSI depends on disk magnetization. Stronger fields lead to slower growths, and the dominant behavior transitions from corrugation to low-amplitude breathing modes.
• Level of VSI turbulence weakens as gas becomes better coupled to magnetic fields as long as the MRI remains suppressed. We speculate weakest turbulence is achieved at some intermediate values of Am.
• The conclusions above apply to locally isothermal disks. Finite thermal relaxation timescale weakens the VSI in a way similar to the pure hydrodynamic case, and eventually suppresses the VSI for τ 0.1P orb for our fiducial disk parameters.
7.1. Discussion The development of MHD winds and of hydrodynamic instabilities have been mostly studied independently in literature. This work shows for the first time that they can coexist. A profound implication of this result is that unlike the conventional thinking turbulence serves to both transport angular momentum transport and provide particle stirring, the picture now becomes that angular momentum transport is mainly driven by magnetized disk winds, whereas particle stirring is due to turbulence. This means that observational constraints on disk turbulence from either direct measurement (e.g., Flaherty et al. 2017Flaherty et al. , 2018Teague et al. 2018) or indirect inference from the thickness of the dust layer (e.g., HL Tau; Pinte et al. 2016) do not necessarily tell the mechanism of angular momentum transport, as has already been speculated in the case of HL Tau (Hasegawa et al. 2017).
More broadly, magnetic field tend to play a destabilizing role in the ideal MHD limit due to the development of the MRI turbulence. This has been shown by Lyra & Klahr (2011) in the case of subcritical baroclinic instability. For the VSI, we have also seen that it tends to be weaker as ambipolar diffusion weakens, and it is likely eventually taken over by the MRI. Therefore, we may speculate that with a magnetically-driven disk wind launched from disk surface, the largely magneticallyinactive PPD offers a relatively clean environment for hydrodynamic instabilities to develop, and besides the VSI, it is worth investigating the case with other hydrodynamic instabilities.
As the first numerical study incorporating magnetism to investigate the vertical shear instability, we have adopted some simplifications aiming to make the simulations as clean as possible. Meanwhile, our simulations are subject to several limitations both numerically and physically that can be improved in future studies. Numerically, first we employ Ohmic resistivity in all simulations at the midplane to stabilize the current sheet (Appendix B). Properly following the long-term behaviour of the midplane current sheet requires 3D simulations with reasonable resolutions. Second, our 2D simulations inevitably prevent the development of vortices, which have been known to be prevalent in 3D. Third, the inner boundary becomes less stable toward longer time, which is likely related to the secular outward transport (loss) of poloidal magnetic flux, restricting our simulations to relatively short time.
Physically, we have used prescribed vertical profiles of temperature and ambipolar diffusivity, whereas more realistic calculations should properly account for ionization chemistry (e.g., Bai 2017) and radiative processes in the disk (e.g., Stoll & Kley 2016). Moreover, wind properties are sensitive to the physics of the wind launching region, which we have treated using a simple transition calling for further improvement (e.g., Wang et al. 2019). Lastly, we have neglected the Hall effect, which is also expected to play an important role in a broad range of radii, whose behavior depends on the polarity of the net poloidal field with respect to disk rotation (e.g., Wardle 2007). For the aligned case, prior global 2D (Bai 2017;Béthune 2017) and local 2D and 3D (Lesur et al. 2014;Bai 2014Bai , 2015Simon et al. 2015) simulations indicate that the disk is magnetically laminar, with horizontal field substantially amplified by the Hall-shear instability (HSI, Kunz 2008). We might expect the VSI to coexist on top of the wind as found in this work, though it remains to be seen the outcome as the system configuration can be asymmetric about the midplane due to the HSI (Bai 2017). For the anti-aligned case, the disk might undergo bursty turbu-lent motions (Simon et al. 2015), for which further investigation needs to be pursued in global 3D simulations to understand its interplay with the VSI.
Future works are planned to extend the simulations to 3D and to incorporate more physics, which will likely yield more complete and realistic picture of gas dynamics in the outer region of PPDs.
We thank Min-Kai Lin for useful discussions and an anonymous referee for thoughtful comments following a thorough reading of the manuscript. CC acknowledges support from the Natural Science Foundation of China (grants 11573051, 11633006, 11650110427, and 11661161012). XNB acknowledges support from the Youth Thousand Talent program. This work makes use of the High Performance Computing Resource in the Core Facility for Advanced Research Computing at Shanghai Astronomical Observatory.

RESOLUTION STUDY
Firstly, we compare hydrodynamic simulations via four levels of resolution. Model ResH has a domain and resolution similar to those of N13. The limited spatial domain spans over r ∈ [1, 4] and θ ∈ [−5H d /r, 5H d /r], and the resolution achieves 96 cells per H d in r and 108 cells per H d in θ with uniform grid spacing. We do not apply a transition from the disk zone to the wind zone for model ResH. The other three possess large domain size the same as simulations presented in Table 1, reaching 16, 32, and 48 cells per H d at the midplane, respectively. The direct comparison of the evolution of kinetic energies in a box of r ∈ [2, 4] and θ ∈ [−3H d /r, 3H d /r] are shown in the left panel of Figure 15. In the early stages, relaxation from the inner boundary dominates the fluctuations in model Res16, Res32, and Res48. After ∼ 15P 0 , fluctuations by the VSI take over, and the kinetic energies converge for all four runs after ∼ 50P 0 .
We also perform resolution study for MHD simulations but with weak disk magnetization (β 0 = 10 5 ). Three resolutions are tested as in hydrodynamical simulations. The right panel of Figure 15 indicates that the three simulations have their kinetic energies converge after 100P 0 . As experimented, a resolution of 48 cells per H d can be timeconsuming, especially when field strengths are strong. Overall, we conclude that a moderate resolution of 32 cells per H d can be well suited for the purpose of this study.

CORRUGATION OF THE MIDPLANE CURRENT SHEET
In our simulations, the current sheet emerges at the midplane due to the orbital shear generated toroidal magnetic fields. The current sheet dissipates toroidal magnetic field via reconnection. This current sheet is observed to be unstable and leads to fluctuations to the gas motion and magnetic field above/below the midplane. The unstable current sheet can further adversely affect the development of the VSI. As the first numerical experiment to explore the interplay between the VSI and laminar MHD winds, we implement Ohmic resistivity in the midplane to stabilize the current sheet, allowing for the VSI to develop on top of a clean background.
To track the behavior of the current sheet, we conduct test runs without the VSI to explore how this current sheet will evolve and affect its vicinity. Figure 16 shows snapshots of toroidal magnetic fields overlaid with equally spaced contour lines of poloidal magnetic flux with τ = 1 and Am = 0.5. Now the resistivity at the midplane is removed while it is still retained in the buffer zone at r < 1.5. Noticeable corrugations of current sheet in vertical direction is detected after ∼ 200P 0 . As seen in later times, the level of corrugations exacerbate and break up the current sheet for R > 2. We are currently working on 3D simulations to further investigate this issue.