High-resolution models of the vertical shear instability

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