Vertical Shear Instability: Mechanisms & Implications

• Vertical Shear Instability (VSI) is a hydrodynamic instability in differentially rotating, vertically stratified disks driven by vertical shear and rapid cooling.
• VSI studies employ inertial-wave analyses and global simulations to reveal how vertical thermal structures and dust properties modulate turbulence and angular momentum transport.
• The instability’s effects include layered turbulence, altered dust settling, and observable kinematic signatures with velocity deviations up to ~100 m/s.

Searching arXiv for the cited VSI papers to ground the article in current arXiv records. The vertical shear instability (VSI) is a purely hydrodynamic instability of differentially rotating, vertically stratified disks that taps the free energy stored in a vertical gradient of orbital angular velocity, $\Omega(R,z)$. In protoplanetary disks, it is especially relevant where magnetic turbulence is often suppressed and where thermal relaxation is sufficiently rapid to weaken stabilizing buoyancy (Lin et al., 2015). Recent work has shifted the emphasis from idealized, vertically uniform cooling toward vertically structured thermal relaxation, thermally stratified atmospheres, radiative diffusion, dust–gas thermal coupling, and convergence in fully global simulations. Across these developments, a consistent picture emerges: VSI is best understood as a cooling-sensitive, vertically global inertial-wave instability whose nonlinear saturation can yield vertically extended turbulence, layered turbulence, corrugation-wave zones, and in some studies long-lived vortices, with strong implications for dust settling, angular-momentum transport, and observability (Fukuhara et al., 2022).

1. Historical placement and physical definition

Vertical shear arises because a radially non-isothermal disk is baroclinic: pressure and density gradients are misaligned, so the angular velocity depends on both radius and height. In the simplified vertically isothermal power-law disk used in several studies, this appears as a nonzero $\partial (R\Omega)/\partial z$, with representative expressions such as

$$\frac{\partial (R\Omega)}{\partial z} = \frac{q}{2}\frac{z}{R}\Omega_{\rm K},$$

or equivalently $\partial\Omega/\partial z \neq 0$ in cylindrical coordinates (Fukuhara et al., 2022). The instability is the disk manifestation of the Goldreich–Schubert–Fricke instability and is therefore often described as a buoyancy-modified shear instability: vertical shear is the energy source, while stable vertical stratification is the principal stabilizing influence (Pfeil et al., 2020).

The modern VSI literature distinguishes it from MRI, convective overstability, zombie vortex instability, and gravitational instability by both its trigger and its thermodynamic requirement. MRI requires magnetic coupling; VSI does not. Convective overstability and subcritical baroclinic instability rely on radial entropy structure and finite cooling, whereas VSI relies on vertical shear and fast enough cooling to reduce effective buoyancy (Barraza-Alfaro et al., 2021). This makes VSI particularly relevant in weakly ionized outer disks and dead zones, where magnetic transport may be weak but radiative relaxation can still be short (Lin et al., 2015).

A recurring result is that VSI should not be treated as a purely local diffusive turbulence mechanism. Several recent analyses instead characterize it as an overstability of inertial waves or inertial–gravity waves, with vertically global body modes and, depending on the thermal and boundary structure, surface-localized branches (Svanberg et al., 2022). This suggests that the detailed vertical thermal structure, surface boundary behavior, and radial extent of the modeled disk are not technical details but part of the instability mechanism itself.

2. Cooling, buoyancy, and linear operating criteria

The key linear constraint is rapid thermal relaxation. Lin and Youdin derived the widely used vertically global criterion

$$t_c < \frac{\Omega_{\mathrm K}^{-1} h |q|}{\gamma -1},$$

or, in dimensionless form,

$$\beta \equiv \Omega_{\mathrm K} t_c < \beta_{\rm crit} \equiv \frac{h |q|}{\gamma - 1},$$

for vertically isothermal irradiated disks (Lin et al., 2015). This criterion is central to later simulation work, including studies with vertically varying cooling, because it provides a simple way to label regions as linearly VSI-unstable or linearly VSI-stable (Fukuhara et al., 2022). A related local criterion compares the vertical shear to the vertical Brunt–Väisälä frequency,

$$\beta \leq \beta_{\rm lc} = \frac{\left|\partial_z (R\Omega)\right|}{N_z^2}\Omega_{\rm K},$$

where

$$N_z^2 = -\frac{1}{\rho_{\rm g} C_P} \frac{\partial P}{\partial z}\frac{\partial s}{\partial z}.$$

This local form is useful for interpreting residual activity near the midplane when the global criterion predicts stabilization (Fukuhara et al., 2022).

Radiative relaxation in protoplanetary disks is not controlled by a single process. Realistic cooling models combine optically thick diffusion, optically thin thermal emission, and gas–dust collisional coupling. A representative formulation is

$$\tau_{\rm relax} = \max\big(\tau_{\rm coll},\,\tau_{\rm diff},\,\tau_{\rm emit}\big),$$

with

$$\tau_{\rm diff} = \frac{1}{f D_E k^2}, \qquad \tau_{\rm thin} = \max(\tau_{\rm emit},\tau_{\rm coll}),$$

and a scale-dependent diffusion time through the perturbation wavenumber $\partial (R\Omega)/\partial z$0 (Pfeil et al., 2020). This has two major consequences. First, the VSI-active layer can be bounded above not by optical depth alone but by dust–gas thermal decoupling in the upper atmosphere. Second, the instability may be strongest only at intermediate radii: Lin and Youdin found it most effective from $\partial (R\Omega)/\partial z$1 AU to $\partial (R\Omega)/\partial z$2 AU with a characteristic growth time of $\partial (R\Omega)/\partial z$3 local orbital periods, and suppressed in both the opaque inner disk and the optically thin outer disk (Lin et al., 2015).

Dust evolution directly modifies these cooling conditions. Dust growth and settling increase cooling times in optically thin outer regions because depletion of small grains weakens thermal emission and gas–dust coupling. One study found that in a default disk model with $\partial (R\Omega)/\partial z$4 solar masses, dust growth from $\partial (R\Omega)/\partial z$5 to $\partial (R\Omega)/\partial z$6 causes a decrease in the VSI growth rate by a factor of more than 10 (Fukuhara et al., 2021). A later semi-analytic treatment found that sufficiently small grains allow a stable equilibrium with VSI-driven diffusion at a level of $\partial (R\Omega)/\partial z$7, whereas increasing the maximum grain size eliminates that equilibrium and leads to runaway settling (Fukuhara et al., 2024). This suggests a two-way feedback: VSI stirs dust, while dust determines whether the disk can cool quickly enough for VSI to persist.

3. Linear mode structure and wave interpretation

A major conceptual advance in recent work is the interpretation of VSI as an inertial-wave instability with distinct mode families. In thermally stratified disks, the instability consists of body modes and surface modes: body modes extend perturbations across the disk, while surface modes are confined to regions of strong shear (Yun et al., 2024). Thermal stratification causes surface modes to bifurcate into two branches, one associated with the strongest shear at mid-height and another near the disk surfaces; the mid-height branch has the higher growth rate (Yun et al., 2024). Surface modes generally require large radial wavenumber $\partial (R\Omega)/\partial z$8, while body-mode growth increases as $\partial (R\Omega)/\partial z$9 decreases (Yun et al., 2024).

The inertial-wave formulation is especially explicit in the wavelike saturation analysis of global 2D disks. There, the relevant $$\frac{\partial (R\Omega)}{\partial z} = \frac{q}{2}\frac{z}{R}\Omega_{\rm K},$$0 corrugation modes obey an inertial-mode dispersion relation

$$\frac{\partial (R\Omega)}{\partial z} = \frac{q}{2}\frac{z}{R}\Omega_{\rm K},$$1

and the VSI is described as selecting and amplifying low-frequency inertial waves with a characteristic preferred frequency

$$\frac{\partial (R\Omega)}{\partial z} = \frac{q}{2}\frac{z}{R}\Omega_{\rm K},$$2

for representative choices of $$\frac{\partial (R\Omega)}{\partial z} = \frac{q}{2}\frac{z}{R}\Omega_{\rm K},$$3 and vertical shear parameter $$\frac{\partial (R\Omega)}{\partial z} = \frac{q}{2}\frac{z}{R}\Omega_{\rm K},$$4 (Svanberg et al., 2022). In that framework, the saturated disk is partitioned into radial wave zones, each associated with an approximately constant frequency and bounded by turning points or corotation-like resonances (Svanberg et al., 2022). This wave-zonal structure reappears in later fully 3D high-resolution models, though with less coherence than in 2D (Lesur et al., 28 Aug 2025).

Boundary behavior matters because the unstable modes are vertically global waves. Studies with partially reflecting upper and lower boundaries show that growth rates of unstable modes diminish when the outgoing component of the flow is greater than the incoming one for high-order body modes, and that when the outgoing wave component dominates, the growth of VSI is notably suppressed (Wu et al., 2024). The same work finds unstable non-axisymmetric modes whose growth rate decreases with azimuthal wavenumber, and unequal boundary conditions at the two disk surfaces naturally produce non-symmetric modes relative to the midplane (Wu et al., 2024). This supports the view that VSI is fundamentally a wave-trapping problem, not simply a local shear instability.

4. Nonlinear saturation and vertically structured states

The nonlinear outcome depends strongly on vertical thermal structure. In disks with vertically varying cooling times, 2D global simulations identify two saturated states (Fukuhara et al., 2022). The first are T states, characterized by vertical turbulent motion penetrating into the linearly VSI-stable midplane layer. The second are pT states, characterized by turbulent motion confined in the unstable layers above and below the midplane. The transition is controlled by the thicknesses of the stable and unstable layers: pT states are realized when the midplane VSI-stable layer is thicker than two gas scale heights, and VSI-driven turbulence is largely suppressed at all heights when the unstable region above and below the midplane is thinner than two gas scale heights (Fukuhara et al., 2022).

The same study formalizes this using the thicknesses

$$\frac{\partial (R\Omega)}{\partial z} = \frac{q}{2}\frac{z}{R}\Omega_{\rm K},$$5

where $$\frac{\partial (R\Omega)}{\partial z} = \frac{q}{2}\frac{z}{R}\Omega_{\rm K},$$6 is the top of the stable midplane layer and $$\frac{\partial (R\Omega)}{\partial z} = \frac{q}{2}\frac{z}{R}\Omega_{\rm K},$$7 is the top of the full unstable region (Fukuhara et al., 2022). The practical thresholds are simple: $$\frac{\partial (R\Omega)}{\partial z} = \frac{q}{2}\frac{z}{R}\Omega_{\rm K},$$8 yields pT states, while $$\frac{\partial (R\Omega)}{\partial z} = \frac{q}{2}\frac{z}{R}\Omega_{\rm K},$$9 suppresses turbulence globally (Fukuhara et al., 2022). This extends earlier “sandwich mode” results, in which realistic thermal relaxation around the water ice line produced a turbulent VSI-active layer roughly within $\partial\Omega/\partial z \neq 0$0, bounded above by a quiescent atmosphere where dust–gas collisional decoupling makes cooling too slow (Pfeil et al., 2020).

Thermal stratification in full 3D simulations can strengthen the instability considerably. In vertically isothermal and thermally stratified disks, the latter develop enhanced shear and correspondingly stronger turbulence. One recent study reports that in thermally stratified disks the turbulence stress reaches $\partial\Omega/\partial z \neq 0$1, more than an order of magnitude stronger than the isothermal case, with saturation more pronounced near the disk surfaces than at the midplane (Yun et al., 2024). A plausible implication is that the classical isothermal VSI picture underestimates the transport and observability of the instability where irradiated atmospheres steepen vertical shear.

Saturation mechanisms remain an active topic. Two-moment radiation-hydrodynamical simulations identify secondary instabilities parasitic to VSI. In those calculations, axisymmetric VSI modes and baroclinic torque create bands of nearly uniform specific angular momentum, while high-shear layers between bands trigger Kelvin–Helmholtz instability. The resulting transfer of kinetic energy to small-scale eddies limits the maximum energy of the VSI modes and likely causes saturation (Fuksman et al., 2023). The same work identifies a third mechanism, amplification of eddies by baroclinic torques, producing meridional vortices with Mach numbers up to $\partial\Omega/\partial z \neq 0$2 (Fuksman et al., 2023). This suggests that VSI saturation can differ qualitatively depending on whether one emphasizes inertial-wave organization, parasitic KHI, or baroclinic amplification; these are not necessarily contradictory, since they may operate in different thermal and numerical regimes.

5. Turbulence, stresses, convergence, and interaction with other physics

Reported transport levels vary, but several consistent ranges appear. In a global 3D locally isothermal simulation used for synthetic observations, the steady state gives $\partial\Omega/\partial z \neq 0$3 (Barraza-Alfaro et al., 2021). In models with realistic thermal relaxation around 5 au, volume-averaged $\partial\Omega/\partial z \neq 0$4 and peak $\partial\Omega/\partial z \neq 0$5 at $\partial\Omega/\partial z \neq 0$6 are reported (Pfeil et al., 2020). In the vertically varying cooling study, the strongest cases reach $\partial\Omega/\partial z \neq 0$7, while fully developed T states have typical midplane vertical velocity dispersion $\partial\Omega/\partial z \neq 0$8 (Fukuhara et al., 2022). High-resolution parameter studies further show that the saturated stress depends strongly on aspect ratio, with $\partial\Omega/\partial z \neq 0$9 in full $$t_c < \frac{\Omega_{\mathrm K}^{-1} h |q|}{\gamma -1},$$0 3D simulations (Manger et al., 2020).

A central recent issue has been convergence. Fully 3D radially extended global models at up to 200 points per scale height find that large-scale transport properties are converged with 100 points per scale height, leading to a Shakura-Sunyaev alpha $$t_c < \frac{\Omega_{\mathrm K}^{-1} h |q|}{\gamma -1},$$1 in the self-similar region of the disk (Lesur et al., 28 Aug 2025). The same work emphasizes that inner boundary condition artifacts propagate deep inside the computational domain and reduce $$t_c < \frac{\Omega_{\mathrm K}^{-1} h |q|}{\gamma -1},$$2 in those regions, implying that short radial domains can underestimate or radially bias VSI transport (Lesur et al., 28 Aug 2025). This result is directly relevant to earlier disagreements in the literature over stress levels and the presence of long-lived coherent structures.

Whether VSI robustly generates large vortices remains unsettled. A high-resolution parameter study with 18 cells per $$t_c < \frac{\Omega_{\mathrm K}^{-1} h |q|}{\gamma -1},$$3 and full $$t_c < \frac{\Omega_{\mathrm K}^{-1} h |q|}{\gamma -1},$$4 azimuth found consistent formation of large scale vortices across all investigated parameters, with uniformly aspect ratios of $$t_c < \frac{\Omega_{\mathrm K}^{-1} h |q|}{\gamma -1},$$5 and radial widths of approximately $$t_c < \frac{\Omega_{\mathrm K}^{-1} h |q|}{\gamma -1},$$6 (Manger et al., 2020). Earlier global 3D work likewise found that the VSI is capable of forming large vortices which can exist at least several hundred orbits in simulations covering a disk with $$t_c < \frac{\Omega_{\mathrm K}^{-1} h |q|}{\gamma -1},$$7, suggesting that VSI can trigger vortex formation via a secondary instability such as Rossby Wave Instability or Kelvin Helmholtz Instability (Manger et al., 2018). By contrast, the 2025 high-resolution global models report no sign of long-lived zonal flows, pressure bumps or vortices, in contrast to lower-resolution simulations (Lesur et al., 28 Aug 2025). This discrepancy is one of the main current controversies. A plausible reading of the published results is that vortex survival is sensitive to resolution, thermodynamics, and perhaps the role of finite cooling times, since long-lived vortices are more prominent in some non-isothermal or moderate-resolution studies than in the highest-resolution locally isothermal calculations.

VSI also interacts with embedded planets and with magnetic fields. In disks with planets, the instability can contribute to the accretion stress, but depending on disk conditions, an embedded planet can coexist with or suppress VSI turbulent stress (Ziampras et al., 2022). Spiral shocks and planet-generated vortices interfere with VSI near the planet, with the instability recovering at large enough distances (Ziampras et al., 2022). In magnetized disks, toroidal or poloidal fields stabilize VSI, with surface modes vanishing before body modes as magnetization increases; Ohmic resistivity can revive subdued modes, and in poloidal-field comparisons the VSI dominates over MRI for Ohmic Elsässer numbers $$t_c < \frac{\Omega_{\mathrm K}^{-1} h |q|}{\gamma -1},$$8 at plasma beta parameter $$t_c < \frac{\Omega_{\mathrm K}^{-1} h |q|}{\gamma -1},$$9 (Cui et al., 2021). This implies that VSI remains relevant even in weakly magnetized disks, provided non-ideal MHD sufficiently diffuses field-line tension.

6. Dust evolution, observability, and broader significance

Because VSI turbulence is strongly anisotropic, vertical diffusion is often more consequential than the nominal radial stress. In the vertically varying cooling study, the inferred vertical diffusion coefficient is written as

$$\beta \equiv \Omega_{\mathrm K} t_c < \beta_{\rm crit} \equiv \frac{h |q|}{\gamma - 1},$$0

and previous simulations imply a wide range of correlation times, yielding $$\beta \equiv \Omega_{\mathrm K} t_c < \beta_{\rm crit} \equiv \frac{h |q|}{\gamma - 1},$$1 for fully developed T-state VSI at the midplane (Fukuhara et al., 2022). In the self-consistent dust-settling model, the equilibrium active state corresponds to $$\beta \equiv \Omega_{\mathrm K} t_c < \beta_{\rm crit} \equiv \frac{h |q|}{\gamma - 1},$$2, while larger grains eliminate the equilibrium and permit runaway settling (Fukuhara et al., 2024). The practical consequence is that VSI can either maintain a vertically extended small-grain layer or collapse into a state that favors very thin dust sublayers and possibly planetesimal formation.

Observationally, the best-developed signatures are kinematic. Synthetic CO observations of a global 3D VSI-unstable disk show a steady-state stress of $$\beta \equiv \Omega_{\mathrm K} t_c < \beta_{\rm crit} \equiv \frac{h |q|}{\gamma - 1},$$3 and large-scale velocity deviations of $$\beta \equiv \Omega_{\mathrm K} t_c < \beta_{\rm crit} \equiv \frac{h |q|}{\gamma - 1},$$4 as axisymmetric rings (Barraza-Alfaro et al., 2021). Optimal conditions are at inclinations $$\beta \equiv \Omega_{\mathrm K} t_c < \beta_{\rm crit} \equiv \frac{h |q|}{\gamma - 1},$$5, and the work argues that standard diagnostics based on non-thermal broadening are not applicable to anisotropic VSI turbulence (Barraza-Alfaro et al., 2021). Thermally stratified disks strengthen these signatures: synthetic velocity residual maps show axisymmetric rings in isothermal disks and ring segments in thermally stratified disks, visible at inclinations as high as $$\beta \equiv \Omega_{\mathrm K} t_c < \beta_{\rm crit} \equiv \frac{h |q|}{\gamma - 1},$$6, with amplitudes from $$\beta \equiv \Omega_{\mathrm K} t_c < \beta_{\rm crit} \equiv \frac{h |q|}{\gamma - 1},$$7 to $$\beta \equiv \Omega_{\mathrm K} t_c < \beta_{\rm crit} \equiv \frac{h |q|}{\gamma - 1},$$8 at a $$\beta \equiv \Omega_{\mathrm K} t_c < \beta_{\rm crit} \equiv \frac{h |q|}{\gamma - 1},$$9 inclination (Yun et al., 2024). Since optically thick tracers probe the disk surface where thermally stratified VSI is strongest, the observed residual grows with tracer optical depth (Yun et al., 2024).

At the same time, other work cautions that VSI signatures are unlikely in disks that contain massive non-axisymmetric features. In simulations with embedded planets, spiral shocks and vortices both weaken or quench the instability and also dominate the observable kinematics, making a clean VSI diagnosis difficult (Ziampras et al., 2022). This suggests that the most favorable targets are smooth, approximately axisymmetric disks without strong planet-driven spirals.

The broader significance of VSI lies in its role as a coupling point among thermal physics, dust evolution, and hydrodynamic transport. Rapidly cooling, weakly magnetized regions can sustain vertically global turbulence and wave activity; dust growth and settling can then reduce cooling efficiency and confine or suppress that turbulence; the resulting changes in diffusion alter dust scale heights and planetesimal-forming conditions; and the kinematic outcome is, in principle, observable with sufficiently high spectral and angular resolution. Across current arXiv work, the most defensible synthesis is therefore not that VSI has a single universal saturated state, but that it occupies a family of regimes—globally turbulent, partially turbulent, layered, wave-dominated, or suppressed—set primarily by vertical thermal structure, cooling physics, and the efficiency with which waves and secondary instabilities redistribute energy (Fukuhara et al., 2022).