Toward First-Principles Multi-Messenger Predictions: Coupling Nuclear Networks with GR Radiation-MHD in { t Gmunu}

We present a new implementation of nuclear reaction networks in the \texttt{G}eneral-relativistic \texttt{mu}ltigrid \texttt{nu}merical (\texttt{Gmunu}) code, a framework for general relativistic radiation magnetohydrodynamics (GRRMHD). The extended code self-consistently evolves nuclear species coupled to hydrodynamics, magnetic fields, and neutrino radiation transport under the conformal flatness approximation to Einstein’s equations. Four approximate nuclear networks are included, with stiff source terms integrated using implicit-explicit Runge-Kutta schemes. Validation is performed through benchmarks including conserved-to-primitive recovery with a tabulated stellar equation of state, one-zone silicon burning, and hydrodynamic tests of shock tubes, acoustic pulses, and detonation fronts of Type Ia supernovae. These tests confirm accurate coupling between nuclear reactions and fluid dynamics, conserving electron and nuclear mass fractions to machine precision. As an application, we conduct spherically symmetric core-collapse supernova simulations. The models reproduce the expected non-exploding behavior of standard progenitors, while enhanced neutrino heating revives the shock. Including nuclear burning modifies the post-shock composition and dynamics, converting silicon and oxygen layers into iron-group nuclei and strengthening the explosion. This demonstrates the impact of explosive burning on ejecta composition and shock evolution, and establishes the stability of the coupled GR radiation-MHD-nuclear framework. The implementation is fully compatible with multidimensional GRMHD simulations and represents the first GRRMHD code combining M1 neutrino transport with fully coupled nuclear burning.

💡 Research Summary

The paper presents a comprehensive implementation of nuclear reaction networks within the General‑relativistic Multigrid Numerical (Gmunu) code, which already solves general‑relativistic radiation magnetohydrodynamics (GRRMHD) using the conformal‑flatness approximation. Four approximate networks (7‑, 13‑, 19‑, and 21‑species) are incorporated, and the stiff source terms arising from nuclear burning are integrated with implicit‑explicit (IMEX) Runge‑Kutta schemes. The authors detail the mathematical formulation: the mass‑fraction conservation equation is cast in a flux‑source form, with source terms treated implicitly while the hydrodynamic fluxes remain explicit. An analytical Jacobian is supplied by the network modules, and a multidimensional Broyden solver (falling back to Newton‑Raphson) resolves the resulting nonlinear system each timestep.

Thermodynamics are handled by coupling a Helmholtz‑type stellar equation of state (EoS) to a tabulated nuclear statistical equilibrium (NSE) EoS. The internal energy includes both thermal contributions and mass‑excess terms, ensuring that nuclear reactions modify only the thermal part while keeping the total specific internal energy consistent. A temperature‑based blending scheme (T_lo = 5 GK, T_hi = 5.8 GK) smoothly transitions between the stellar and NSE EoS, with consistency checks on electron fraction Y_e and composition X_i to trigger NSE recomputation when necessary.

A hierarchy of verification tests validates the implementation. First, the conserved‑to‑primitive conversion with the tabulated EoS is shown to be accurate. Second, a one‑zone silicon‑burning test reproduces expected temperature and composition evolution while conserving total energy to machine precision. Third, classic Newtonian benchmarks—shock tubes, acoustic pulses, and a Type Ia supernova detonation front—demonstrate that the coupled nuclear‑hydrodynamics correctly captures shock propagation, burning suppression in unresolved shocks, and the formation of a detonation wave. The shock‑detection algorithm uses a pressure‑gradient‑velocity criterion (f_shock = 2/3) to temporarily disable burning in cells identified as part of a numerical shock, preventing spurious energy release.

The authors then apply the framework to spherically symmetric core‑collapse supernova (CCSN) simulations of a 15 M_⊙ progenitor. In the baseline model with standard neutrino heating, the shock stalls and no explosion occurs, reproducing known non‑exploding behavior. When neutrino heating is artificially enhanced to revive the shock, the inclusion of nuclear burning markedly alters the post‑shock region: silicon‑ and oxygen‑rich layers are rapidly processed into iron‑group nuclei (e.g., ^56Ni), releasing additional nuclear binding energy. This extra energy boosts the pressure behind the shock, leading to a stronger, earlier explosion and a higher yield of heavy elements in the ejecta. The results illustrate that explosive burning can materially affect both the dynamics and the nucleosynthetic output of CCSNe, especially when coupled to neutrino‑driven mechanisms.

Although magnetic fields are not evolved in the presented 1‑D runs, the code architecture fully supports multidimensional GRMHD and M1 neutrino transport, making the implementation the first GRRMHD code that couples M1 neutrino radiation with on‑the‑fly nuclear burning. The authors discuss future extensions, including larger reaction networks (thousands of isotopes), full 3‑D simulations of CCSNe, neutron‑star mergers, and Type Ia supernovae, where the interplay of gravity, neutrinos, magnetic fields, and nuclear synthesis will be essential for accurate multi‑messenger predictions. In summary, the paper delivers a technically robust, validated, and extensible platform that integrates stiff nuclear physics into fully relativistic radiation‑magnetohydrodynamics, opening new avenues for predictive, first‑principles modeling of astrophysical explosions.