Technical Summary: Finite-volume scheme for first-order viscoresistive relativistic magnetohydrodynamics

Problem Statement

Relativistic magnetohydrodynamics (MHD) is essential for modeling conducting fluids in high-energy and astrophysical environments, such as quark-gluon plasmas and accretion flows. A persistent challenge in this field is the incorporation of dissipative transport coefficients (viscosity, thermal conduction, and resistivity) in a manner that preserves causality and stability. Traditional first-order theories yield parabolic equations that violate causality and induce instabilities at high velocities. While Müller–Israel–Stewart (MIS) theories resolve this by promoting dissipative currents to independent dynamical variables, this approach significantly increases complexity and the number of degrees of freedom.

Recent developments in Bemfica–Disconzi–Noronha–Kovtun (BDNK) theory offer a minimal, first-order alternative that restores causality and stability by including higher-derivative corrections to conserved currents without introducing new variables. However, extending BDNK theory to resistive MHD presents a specific difficulty: while the "one-form" formulation of MHD (treating magnetic flux as a conserved charge) removes the stiffness associated with Ampère's law, it renders the resistive dynamics parabolic. Incorporating BDNK corrections into this coupled, dissipative, one-form MHD framework within a stable and efficient numerical scheme has remained an open gap.

Methodology

The authors present a numerical implementation of first-order viscoresistive relativistic MHD based on the BDNK framework, utilizing the BHAC code. The methodology involves:

1. Theoretical Formulation:

   • The authors adopt the ultra-relativistic limit (neglecting mass density) and the one-form approach, where magnetic field lines are treated as conserved string-like objects.
   • Constitutive relations for the stress-energy tensor ($T^{\mu\nu}$) and dual field strength tensor ($J^{\mu\nu}$) include first-order dissipative corrections for thermal conductivity, viscosity, and resistivity.
   • To ensure causality, BDNK corrections are applied via a zeroth-order on-shell transformation. This introduces relaxation times ($\tau_\epsilon, \tau_u, \tau_b$) that convert parabolic diffusion equations into hyperbolic telegrapher-type equations.
   • Crucially, the authors identify that three BDNK terms are insufficient for the coupled ultra-relativistic energy-momentum sector. A fourth term, controlled by the coefficient $\tau_X$, is required to ensure the system remains hyperbolic and causal, preventing elliptic behavior in the sound sector.

2. Numerical Scheme:

   • The equations are projected onto space-like hypersurfaces and decomposed into conserved variables ($U_I$) and fluxes ($F_I$), alongside primitive variables ($P_I$) and source terms ($S_I$).
   • The scheme employs a finite-volume discretization with second-order total variation diminishing (TVD) time-stepping and spatial reconstruction.
   • Dissipative fluxes are evaluated using second-order central differencing of reconstructed interface states, with upwinded numerical fluxes computed via an approximate Riemann solver.
   • Primitive Variable Recovery: A key feature of the implementation is an efficient algebraic method for recovering primitive variables. Because the BDNK corrections make the conserved variables depend linearly on the time derivatives of the primitive variables, the recovery reduces to a local matrix inversion ($S_I = \sum (M^{-1})_{IJ} [U'_J - U_J(0)]$). This avoids the need for iterative non-linear solvers typically required in dissipative schemes.
   • The code utilizes constrained transport (FCT) algorithms to maintain the divergence-free condition of the magnetic field ($\partial_i J^{it} = 0$) to machine precision.

Key Contributions

• Extension of BDNK to Resistive MHD: The paper constructs a numerical scheme that combines viscous and resistive dissipation within the BDNK framework for relativistic MHD.
• Identification of a Necessary Fourth Term: The authors demonstrate that while minimal BDNK terms can create telegrapher-type equations for individual variables, a fourth term ($\tau_X$) is strictly necessary to guarantee the hyperbolicity of the coupled energy-momentum system in the ultra-relativistic limit.
• Efficient Primitive Recovery: The work demonstrates that BDNK corrections allow for a non-iterative, algebraic primitive-variable recovery procedure, significantly enhancing computational efficiency compared to traditional dissipative schemes.
• Robust Implementation: The scheme is implemented in the BHAC code and validated across a range of dissipative regimes, from near-ideal to diffusion-dominated.

Results

The implementation was validated through a series of one- and two-dimensional benchmarks:

1. Analytical Benchmark: A boosted telegrapher solution for magnetic field evolution was used to verify correctness. The numerical results matched the analytical solution with second-order convergence.
2. 1D Shock Tubes: Simulations of shock tubes with varying resistivity ($r_b$) and diffusion parameters ($D_u, D_\epsilon$) showed the scheme remains stable even with discontinuous initial data. The results correctly reproduced ideal MHD behavior for low dissipation and diffusion-dominated dynamics for high dissipation.
3. Kelvin-Helmholtz (KH) Instability: 2D simulations demonstrated that increasing momentum diffusion ($D_u$) suppresses the instability, while increasing resistivity ($r_b$) enhances it by relaxing the coupling between fluid velocity and magnetic field lines. The scheme maintained causal front velocities throughout the evolution.
4. Orszag-Tang Vortex: This test highlighted the critical role of the fourth BDNK term ($\tau_X$). When $\tau_X = 0$, the simulation exhibited anti-diffusive behavior and crashed due to complex front velocities. With $\tau_X > 0$, the system remained stable with real, subluminal front velocities.
5. Harris Sheet: A perturbed double Harris sheet test, featuring thin magnetic reversal layers and large gradients, was simulated using adaptive mesh refinement. The scheme successfully resolved current sheet formation and deformation, maintaining a maximum front velocity of $v_{max} \approx 0.999$, confirming causal evolution.

Significance

The paper claims that this work provides a concrete bridge between recent theoretical developments in relativistic hydrodynamics (specifically BDNK theory) and their application to realistic simulations of magnetized relativistic flows. By demonstrating that dissipative relativistic MHD can be implemented in a numerically robust, causal, and efficient manner without introducing additional dynamical variables, the authors establish a foundation for simulating complex phenomena in environments such as weakly coupled quark-gluon plasmas, the early universe, and magnetar magnetospheres. The ability to independently control resistive and viscous scales (fixing the magnetic Prandtl number) is highlighted as a crucial capability for studying phenomena like magnetic dynamos and magnetar flares. The authors also note that the algebraic structure of their primitive recovery offers a promising pathway for extending these causal viscoresistive corrections to general-relativistic MHD simulations.