High-resolution models of the vertical shear instability

¹ Univ. Grenoble Alpes, CNRS, IPAG, 38000 Grenoble, France
² Department of Applied Mathematics and Theoretical Physics, Centre for Mathematical Sciences, University of Cambridge, Wilberforce Road, Cambridge CB3 0WA, UK

^★ Corresponding author: geoffroy.lesur@univ-grenoble-alpes.fr

Received: 13 June 2025
Accepted: 8 August 2025

Abstract

Context. The vertical shear instability (VSI) is a promising mechanism for generating turbulence and transporting angular momentum in magnetically decoupled regions of 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.

Aims. We address the question of VSI saturation and associated turbulence using radially extended, fully 3D global disc models with very high resolution in the locally isothermal approximation, to capture both the largest VSI scales and the small-scale turbulent cascade.

Methods. We used the GPU-accelerated code Idefix to achieve resolutions of up to 200 points per scale height in the three spatial directions. We chose numerical techniques that minimise numerical diffusion as much as possible: third-order reconstruction schemes, orbital advection, and a third-order time integrator. We modelled the VSI in disc domains extending up to R_out/R_in = 7 in the highest resolution case and R_out/R_in = 25 in our intermediate model (100 points per scale height), with a full 2 π azimuthal extent and a disc aspect ratio H/R = 0.1.

Results. We demonstrate that large-scale transport properties converge with 100 points per scale height, leading to a Shakura-Sunyaev α=1.3 × 10⁻³ in the bulk of the computational domain. Inner boundary condition artefacts propagate deep inside the computational domain (typically over Δ R/R ∼ 2−3), leading to a reduced α in these regions. The large-scale corrugation wave zones identified in 2D models persist in 3D, albeit with less coherence. Our models show no signs 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 within the framework of critically balanced rotating turbulence. We propose that the small non-axisymmetric scales could be modelled with an effective anisotropic viscosity in 2D simulations, significantly reducing the computational cost of these models while still capturing the important physics.

Conclusions. 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 consistently polluted by boundary-condition artefacts influencing the first VSI wave train. Radially extended domains should therefore be used in a more systematic manner, and more realistic inner boundaries should be explored to mimic the radial structure of real protoplanetary discs.