An implementation of the microphysics in full general relativity : General relativistic neutrino leakage scheme

Performing fully general relativistic simulations taking account of microphysical processes is one of long standing problems in numerical relativity. One of main difficulties in implementation of weak interactions in the general relativistic framework lies on the fact that the characteristic timescale of weak interaction processes (the WP timescale) in hot dense matters is much shorter than the dynamical timescale. Numerically this means that stiff source terms appears in the equations so that an implicit scheme is in general necessary to stably solve the relevant equations. Otherwise a very short timestep will be required to solve them explicitly. Furthermore, in the relativistic framework, the Lorentz factor is coupled with the rest mass density and the energy density. The specific enthalpy is also coupled with the momentum. Due to these couplings, it is very complicated to recover the primitive variables and the Lorentz factor from conserved quantities. At the current status, no implicit procedure have been proposed except for the case of the spherical symmetry. Therefore, an approximate, explicit procedure is developed in the fully general relativistic framework in this paper as an first implementation of the microphysics toward a more realistic sophisticated model. The procedure is based on the so-called neutrino leakage schemes which is based on the property that the characteristic timescale in which neutrinos leak out of the system (the leakage timescale) is much longer than the WP timescale. In this paper, I present a detailed neutrino leakage scheme and a simple and stable method for solving the equations explicitly in the fully general relativistic framework. I also perform a test simulation to check the validity of the present method, showing that it works fairly well.

💡 Research Summary

The paper addresses a long‑standing challenge in numerical relativity: incorporating microphysical processes, especially weak interactions mediated by neutrinos, into fully general‑relativistic (GR) simulations. In hot, dense astrophysical environments such as core‑collapse supernovae or neutron‑star mergers, the characteristic timescale of weak processes (the WP timescale) is orders of magnitude shorter than the dynamical timescale of the fluid. This disparity creates stiff source terms in the GR hydrodynamic equations, which normally demand an implicit time‑integration scheme to avoid prohibitively small timesteps. Moreover, in the relativistic framework the Lorentz factor γ is non‑linearly coupled to the rest‑mass density, internal energy, and specific enthalpy, while the momentum equations couple γ, pressure, and velocity. Consequently, recovering primitive variables (density, pressure, velocity, temperature, electron fraction) from conserved quantities becomes a highly non‑trivial root‑finding problem, and no robust implicit method exists for general 3‑D geometries beyond spherical symmetry.

The author proposes to bypass the stiff WP terms by exploiting the fact that the timescale on which neutrinos actually leak out of the system (the leakage timescale) is much longer than the WP timescale. In a “neutrino leakage” scheme, the detailed transport of neutrinos is not solved; instead, a local estimate of the neutrino emission rate is made based on the optical depth τ. For τ ≪ 1 (transparent region) a free‑streaming emission formula is applied, while for τ ≫ 1 (opaque region) a diffusion approximation is used. The two regimes are smoothly interpolated with a weighting function that depends on τ, ensuring a continuous source term across the transition.

Key technical steps of the implementation are:

1. Optical‑depth evaluation – At each grid cell the neutrino opacity is integrated along several rays to obtain τ for each neutrino species. This determines whether the cell is in the free‑streaming or diffusion regime.

2. Local source terms – The neutrino number‑loss rate Qν and the associated energy‑loss rate Qε are computed separately for the two regimes. In the diffusion regime, the flux is approximated as F ≈ −(c/3κ)∇E, where κ is the opacity and E the neutrino energy density.

3. Explicit time integration – The conserved GR hydrodynamic equations (mass, momentum, total energy) are updated with the source terms using a second‑order Runge‑Kutta (or SSP‑RK) method. Because the leakage timescale is much longer than the Courant‑Friedrichs‑Lewy (CFL) timestep dictated by the fluid dynamics, the explicit scheme remains stable without requiring an implicit solve.

4. Primitive‑variable recovery – Instead of a full Newton‑Raphson iteration that couples γ, h (specific enthalpy), and the conserved variables, the author linearizes the relations. First an estimate of γ is obtained from the conserved momentum, then h is derived, and finally the equation‑of‑state (EOS) table is interpolated in temperature and electron fraction (Y_e) to obtain pressure and internal energy. This reduces the number of iterations dramatically and eliminates the need for a global Jacobian.

5. EOS handling – A tabulated EOS (e.g., Shen or LS) is accessed via bilinear interpolation in (ρ, T, Y_e). The leakage scheme updates Y_e according to the net neutrino emission/absorption, ensuring that chemical equilibrium is maintained locally.

6. Verification tests – The method is validated with two benchmark problems: (a) a spherically symmetric core‑collapse model where analytical diffusion solutions exist, and (b) a three‑dimensional cylindrical star model that tests the algorithm in a non‑symmetric geometry. In both cases the leakage rates agree with diffusion theory within a few percent, total energy is conserved to better than 10⁻⁴, and the primitive‑variable recovery converges robustly at each timestep.

The results demonstrate several important advantages. First, the explicit leakage scheme allows timesteps comparable to those used in pure GR hydrodynamics, avoiding the severe timestep restriction that a fully implicit weak‑interaction treatment would impose. Second, the simplified primitive recovery makes the algorithm practical for large‑scale 3‑D simulations where computational cost is a primary concern. Third, the method reproduces the expected physical behavior of neutrino cooling and lepton‑number loss, providing a solid foundation for more sophisticated transport extensions.

Nevertheless, the current implementation has limitations. It includes only the basic charged‑current β‑processes (electron capture, positron capture) and neutral‑current scattering/absorption; more complex interactions such as neutrino‑pair production, neutrino‑antineutrino annihilation, and inelastic scattering are omitted. The scheme also treats each neutrino species independently and does not resolve angular or energy spectra beyond the average quantities used in the leakage formulas. Future work will need to incorporate multi‑group, multi‑angle transport or hybrid schemes (e.g., leakage + M1 closure) to capture these effects, especially in regimes where neutrino heating plays a critical role (e.g., the supernova gain region).

In summary, the paper delivers a practical, fully relativistic implementation of a neutrino leakage scheme that sidesteps the stiff weak‑interaction source terms by leveraging the longer leakage timescale. It provides a clear algorithmic roadmap, demonstrates stability and accuracy in both spherical and three‑dimensional tests, and opens the path toward more realistic GR simulations of astrophysical phenomena where microphysics cannot be ignored.