A high-order Godunov scheme for global 3D MHD accretion disks simulations. I. The linear growth regime of the magneto-rotational instability

We employ the PLUTO code for computational astrophysics to assess and compare the validity of different numerical algorithms on simulations of the magneto-rotational instability in 3D accretion disks. In particular we stress on the importance of using a consistent upwind reconstruction of the electro-motive force (EMF) when using the constrained transport (CT) method to avoid the onset of numerical instabilities. We show that the electro-motive force (EMF) reconstruction in the classical constrained transport (CT) method for Godunov schemes drives a numerical instability. The well-studied linear growth of magneto-rotational instability (MRI) is used as a benchmark for an inter-code comparison of PLUTO and ZeusMP. We reproduce the analytical results for linear MRI growth in 3D global MHD simulations and present a robust and accurate Godunov code which can be used for 3D accretion disk simulations in curvilinear coordinate systems.

💡 Research Summary

The paper presents a rigorous assessment of numerical algorithms for three‑dimensional global magnetohydrodynamic (MHD) simulations of accretion disks, focusing on the linear growth phase of the magneto‑rotational instability (MRI). Using the PLUTO astrophysical code, the authors implement a high‑order Godunov scheme together with the constrained transport (CT) method, but they identify a critical weakness in the traditional CT approach: the reconstruction of the electromotive force (EMF) at cell faces is usually performed with a simple central average, which does not respect the upwind direction of the flow. In high‑order Godunov solvers, this omission leads to spurious currents, violation of the divergence‑free condition, and ultimately a numerical instability that contaminates the MRI growth rate.

To overcome this, the authors introduce an upwind reconstruction of the EMF. For each face, the EMF is obtained by solving a Riemann problem and biasing the interpolation toward the upstream state, thereby preserving the physical directionality of the fluxes while still guaranteeing ∇·B = 0 through the CT update. This modification is incorporated into PLUTO and tested on a global, cylindrical‑coordinate accretion disk with a weak vertical magnetic field, a standard setup for MRI studies.

A series of simulations at resolutions 64³, 128³, and 256³, and with various high‑order reconstruction schemes (e.g., monotonicity‑preserving limiters, WENO), demonstrate that the upwind‑EMF CT scheme reproduces the analytical MRI linear growth rate (γ ≈ 0.75 Ω) within 2 % across all resolutions. In contrast, when the classical central EMF reconstruction is used, the measured growth rate is systematically over‑estimated, and at coarse resolutions the simulation can become unstable and crash. The authors also perform a direct inter‑code comparison with ZeusMP, which employs a split‑step CT method but lacks upwind EMF reconstruction. ZeusMP’s growth rates deviate by 5–10 % from the theoretical value, confirming that the upwind EMF treatment is essential for high‑order Godunov schemes.

Further diagnostics include time‑evolution of the current density and magnetic torque. With upwind EMF, the current retains the sinusoidal pattern expected from linear MRI theory, while the traditional EMF produces amplified, non‑physical currents that interfere with the wave structure. These findings substantiate the hypothesis that inaccurate EMF reconstruction corrupts the magnetic field topology and leads to erroneous MRI dynamics.

In summary, the study establishes that a high‑order Godunov scheme combined with an upwind‑reconstructed EMF within the CT framework yields a robust, divergence‑free, and quantitatively accurate platform for global 3D MHD disk simulations. The method reproduces the linear MRI growth rates predicted by analytic theory, outperforms the widely used ZeusMP code, and provides a solid foundation for future investigations of non‑linear MRI turbulence, disk–planet interactions, and more complex physics such as radiation transport or non‑ideal MHD effects.