High-resolution models of the vertical shear instability

TL;DR

The paper investigates the nonlinear saturation and turbulence of the vertical shear instability (VSI) in radially extended, fully 3D global protoplanetary-disc models using high-resolution GPU simulations. It demonstrates convergence of large-scale transport at about 100 points per scale height with an outward angular-momentum transport parameter of $\alpha \approx 1.3\times 10^{-3}$, and finds no persistent zonal flows or vortices in 3D, though inner-boundary artefacts can bias results near the inner edge. The turbulent cascade is interpreted through the framework of critically balanced rotating turbulence, with strong anisotropy at large scales and a steep high-$k$ tail likely affected by numerical diffusion; small-scale non-axisymmetric motions can be represented as an effective anisotropic viscosity on the axisymmetric flow. These findings highlight the importance of radially extended domains and motivate 2D closures that capture 3D VSI dynamics, enabling more efficient modelling of disc evolution and angular-momentum transport.

Abstract

(abridged) Context: The vertical shear instability (VSI) is a promising mechanism to generate turbulence and transport angular momentum in protoplanetary discs. While most recent work has focused on adding more complex physics, the saturation properties of the instability in radially extended discs and its convergence as a function of resolution are still largely unknown. We tackle the question of VSI saturation and associated turbulence using radially extended fully 3D global disc models at very high resolution so as to capture both the largest VSI scales and the small-scale turbulent cascade. We use the GPU-accelerated code Idefix to achieve resolutions of up to 200 points per scale height in the 3 spatial directions, with a full 2pi azimuthal extent and disc aspect ratio H/R=0.1. Results: We demonstrate that large-scale transport properties are converged with 100 points per scale height, leading to a Shakura-Sunyaev alpha=1.3e-3. Inner boundary condition artifacts propagate deep inside the computational domain, leading to reduced alpha in these regions. The large-scale corrugation wave zones identified in 2D models survive in 3D, albeit with less coherence. Our models show no sign of long-lived zonal flows, pressure bumps or vortices, in contrast to lower-resolution simulations. Finally, we show that the turbulent cascade resulting from VSI saturation can be interpreted in the framework of critically balanced rotating turbulence. Conclusion: The VSI leads to vigorous turbulence in protoplanetary discs, associated with outward angular momentum transport but without any significant long-lived features that could enhance planet formation. The innermost regions of VSI simulations are always polluted by boundary-condition artifacts affecting the first VSI wave train, so radially extended domains should be used in a more systematic manner.

Figures (20)

• Figure 1: Three-dimensional rendering of the meridional velocity sonic Mach number $v_\theta/c_s$ in simulation EX (200 points per $H$) at $t=800$ orbits after restart. The front part shows the flow in the midplane, while the back part shows the flow at $3H$ above the midplane, with numerical boundaries at $4H$. The domain has been truncated at $R_\mathrm{out}=5.6$ to remove the outer wave-killing zone.
• Figure 2: Radial angular momentum coefficient ($\alpha$) and poloidal kinetic energy ($e_k$) measured at $R=4$ in run HX (solid line) and run LX (dashed line). For clarity, we show the instantaneous $\alpha$ and the time-averaged $\alpha$ over a window of 20 local orbits, $\langle \alpha\rangle_{t,20}$. Run HX was initialized from the saturated state of LX, which explains the absence of an initial transient in HX.
• Figure 3: Angular momentum transport and poloidal kinetic energy ($e_k$) averaged over the last 400 orbits of the simulations EX (solid line), HX (dashed line), and LX (dotted line). The predicted amplification and saturation radius for a single linear wave train with frequency $\omega=0.075\Omega_0$ and $n=1$ from Ogilvie.Latter.ea25 is shown as a dashed green line (eq. \ref{['eq:lin_wave train']}).
• Figure 4: Instantaneous $(R,z)$ cut of the HX model at $t=900$ orbits. Large-scale $n=1$ corrugation modes are evident in $v_z$, reaching sonic velocities (top panel), along with shock fronts propagating radially in the density field (bottom panel).
• Figure 5: Instantaneous $(R,\phi)$ cuts of vertical velocity ($v_z$) and density ($\rho$) in the HX model at $t=900$ orbits. Top: Midplane cut of $v_z$. Middle: $z=0.3R$ (corresponding to 3 pressure scale heights) cut of $v_z$. Bottom: $z=0.3R$ cut of the gas density $\rho$.