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

TL;DR

This work implements in-situ nuclear reaction networks within the GRRMHD code Gmunu, enabling self-consistent evolution of composition, gravity, magnetic fields, and neutrino transport under a conformally flat spacetime. It introduces IMEX Runge–Kutta schemes to handle stiff nuclear source terms, a bridging scheme between NSE and stellar equations of state, and an NSE solver, validated by a sequence of tests and a 1D CCSN application. The results show that explosive burning can modify post-shock composition and contribute meaningfully to explosion energetics, while preserving numerical stability and conservation to machine precision. This framework establishes a first, fully coupled GRRMHD platform with M1 neutrino transport and in-situ nuclear burning, enabling more realistic multi-messenger predictions and paving the way for multidimensional CCSN and related explosive phenomena.

Abstract

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.

Figures (10)

• Figure 1: Relative errors of temperature $T$ (top) and specific energy $\epsilon$ (bottom) from the round-trip tests for five electron fractions $Y_{\rm e}$. The upper-right regions enclosed by the white dashed lines denote where the nuclear EoS is applied, while the stellar or transitional (mixed) EoS is used elsewhere. The mixed regime corresponds to $5 < T < 5.8~{\rm GK}$ where both the stellar and nuclear EoS tables overlap in density. At densities $\rho \gtrsim 10^{7}~{\rm g~cm^{-3}}$ and temperatures $T \lesssim 10^{10}~{\rm K}$, the temperature recovery is imperfect because the code prioritizes solutions from the nuclear EoS (see Section \ref{['sec:eos_bridging']}). Nevertheless, the relative error in the specific energy $\epsilon$, which serves as a primitive variable in hydrodynamical evolution, remains below $10^{-14}$ across the explored parameter space.
• Figure 2: Averaged relative error in conserved-to-primitive tests for selected Lorentz factors $W$, electron fractions $Y_{\rm e}$, and magnetization $\beta_{\rm mag} = 1000$ (weakly magnetized case). The upper-right regions enclosed by the white dashed lines denote where the nuclear EoS is applied, while the stellar or transitional (mixed) EoS is used elsewhere.
• Figure 3: One-zone pure silicon burning test at fixed $\rho = 10^{7}~\rm{g \cdot cm^{-3}}$ and $T = 6\times10^{9}~{\rm K}$, starting from pure $\rm ^{28}Si$. The network evolves the composition toward an iron-group dominated state consistent with nuclear statistical equilibrium. Top panel: Mass fractions of selected species (solid) and NSE predictions (dashed). Middle panel: Energy generation rate. Bottom panels: Deviations from conservation: $\left\vert 1 - \sum X_l \right\rvert$ and $\left\vert 1 - Y_{\rm e}/Y_{\rm e}(0) \right\rvert$.
• Figure 4: Profiles of rest-mass density $\rho$ (upper left), velocity $v$ (upper right), pressure $P$ (lower left), and temperature $T$ (lower right) at $t = 8\times10^{-4}~{\rm s}$ for the Newtonian real gas shock tube test 1 of 2015ApJS..216...31Z. Red dots: numerical results from Gmunu; black solid lines: exact reference solution. The excellent agreement validates the stellar EoS implementation.
• Figure 5: Snapshots of the rest-mass density $\rho$ at six times for the Newtonian real gas acoustic pulse test 2019ApJ...886..105Z2025ApJ...981...63H. The upper panel shows pure hydrodynamics, and the lower panel includes nuclear reactions via the 13-species network aprox13. When nuclear burning is included, the expansion dynamics are visibly altered.