Beyond GRMHD: A Robust Numerical Scheme for Extended, Non-Ideal General Relativistic Multifluid Simulations

TL;DR

This paper introduces a general relativistic multifluid framework that extends GRMHD by allowing multiple relativistic fluid species to interact with a shared electromagnetic field through explicit source-term couplings, thereby incorporating non-ideal effects such as electron inertia and Hall-like dynamics. The authors develop a robust tetrad-based finite-volume solver, using modified HLLC/Roe Riemann solvers and Strang splitting to evolve fluid and field variables in curved spacetime, with a strong emphasis on primitive variable reconstruction and divergence control. They validate the approach across 1D, 2D, and 3D tests in both black hole and neutron star spacetimes, demonstrating convergence to SRMHD and GRMHD limits as appropriate, and showing the ability to sustain high Lorentz factors and magnetizations where GRMHD struggles. The results reveal non-ideal phenomena such as charge separation, large parallel electric fields, and enhanced poloidal fields and reconnection in ergospheres, highlighting the potential for self-consistent jet formation and more realistic electromagnetic modeling around compact objects. The work paves the way for fully dynamic spacetime coupling, implicit/IMEX schemes for stiff source terms, and inclusion of pair production, enabling more accurate simulations of black hole jets, neutron star magnetospheres, and related high-energy phenomena.

Abstract

The equations of general relativistic magnetohydrodynamics (GRMHD) have become the standard mathematical framework for modeling high-energy plasmas in curved spacetimes. However, the fragility of the primitive variable reconstruction operation in GRMHD, as well as the difficulties in maintaining strong hyperbolicity of the equations, sharply limit the applicability of the GRMHD model in scenarios involving large Lorentz factors and high magnetizations, such as around neutron stars. Non-ideal effects, such as electron inertia and Hall terms, are also neglected, and the absence of an explicitly evolved electric field precludes the self-consistent modeling of the strong poloidal fields found around spinning black holes, which are known to be crucial for jet formation. Here, we present a general relativistic multifluid model which strictly generalizes the GRMHD equations, consisting of an arbitrary number of relativistic fluid species interacting with a shared electromagnetic field via an explicit coupling of their source terms, thus allowing for the incorporation of non-ideal effects. We sketch how our model may be derived from general relativistic kinetics (via moments of the relativistic Boltzmann-Vlasov equation), as well as how GRMHD may be recovered in the single-fluid limit as the mobility of charge carriers goes to infinity. We present a numerical scheme for solving the general relativistic multifluid equations, and validate it against the analogous scheme for the GRMHD equations. Since the primitive variable reconstruction operation for our multifluid model is purely hydrodynamic, and therefore independent of the magnetic field, the resulting solver is highly robust, and able to simulate significantly larger Lorentz factors and higher magnetizations (across both black hole and neutron star spacetimes) than GRMHD without loss of either accuracy or stability.

Figures (22)

• Figure 1: The relativistic mass density ${\rho W}$ at time ${t = 0.4}$ for the relativistic Brio-Wu shock tube problem, showing the two asymptotic reference solutions in the ideal MHD (${r_L \to 0}$) limit and the gas dynamic (${r_L \to \infty}$) limit. The wave structure of the MHD solution has been labeled to show the left-moving fast rarefaction wave (FR), the left-moving slow compound wave (SC), the contact discontinuity (CD), the right-moving slow shock (SS), and the right-moving fast rarefaction wave (FR).
• Figure 2: The total relativistic mass densities ${\rho W}$ at time ${t = 0.2}$ for the relativistic Brio-Wu shock tube problem, with ion Larmor radii ${r_L = 10, 1, 0.1, 0.01, 0.001, 0.0001}$. The ideal MHD (${r_L \to 0}$) and gas dynamic (${r_L \to \infty}$) asymptotic reference solutions are also shown in each case. At ${r_L = 10}$, the electrostatic interactions between electron and ion fluids have caused the two-fluid solution to depart slightly from the gas dynamic solution, marginally increasing the speed of the left-moving rarefaction wave and introducing an additional left-moving wave which propagates faster than the MHD waves. At ${r_L = 1}$, the right-moving shock from the gas dynamic solution has slowed considerably and is now much closer in speed to the right-moving slow MHD shock. The slow compound wave is beginning to form, and there is additional left-moving wave structure which is propagating faster than the MHD waves. At ${r_L = 0.1}$, the left-moving wave structure is propagating slower than for ${r_L = 1}$ but still faster than the MHD waves, and the speed of the left-moving fast rarefaction wave has increased even further. The slow compound wave is now approaching the MHD solution. At ${r_L = 0.01}$, both the left-moving fast rarefaction wave and the right-moving slow shock have begun to slow down, with both now approaching the MHD solution. The left-propagating wave structure has slowed further, though is still faster than the MHD waves, as compared to ${r_L = 0.1}$. By ${r_L = 0.001}$ and ${r_L = 0.0001}$, the two-fluid solution has now almost fully converged to the MHD solution.
• Figure 3: The total relativistic mass density ${\rho W}$ at time ${t = 0.2}$ for the relativistic Noh shock tube problem, with ion Larmor radii ${r_L = 0.1}$ and ${r_L = 0.01}$. The ideal MHD (${r_L \to 0}$) asymptotic reference solution is also shown. Due to the high Lorentz factor, the 2D Noble et al. primitive variable reconstruction scheme for relativistic MHD periodically fails, and the solver falls back to the effective 1D Newman-Hamlin method.
• Figure 4: The total relativistic mass density ${\rho W}$ at time ${t = 0.2}$ for the highly magnetized variant of the relativistic Noh shock tube problem, with ion Larmor radii ${r_L = 0.1}$, ${r_L = 0.01}$, and ${r_L = 0.001}$. Due to the combination of high magnetization and moderately high initial Lorentz factor, even the effective 1D Newman-Hamlin primitive variable reconstruction scheme for relativistic MHD fails, so no numerical MHD solution exists for this problem. However, by setting the ion Larmor radius to ${r_L = 0.001}$, we are nevertheless able to ascertain approximately what the numerical relativistic MHD solution would be, if it could be obtained.
• Figure 5: Magnetic flux surfaces at time ${t = 50}$ for the magnetospheric Wald problem in axisymmetry for a spinning (Kerr) black hole with ${a = 0.95}$, obtained using the GRMHD solver (i.e. ${r_L \to 0}$) on the left, and the general relativistic multifluid solver (with initial ion Larmor radius ${r_L = 0.1}$) on the right. We see a larger induced Poynting flux in the multifluid solution, particularly in the poloidal direction, with a larger proportion of the flux surfaces intersecting the horizon of the black hole as compared to the MHD solution.