$r$-process Heating Feedback on Disk Outflows from Neutron Star Mergers

Neutron star mergers produce $r$-process elements, with yields that are sensitive to the kinematic and thermodynamic properties of the ejecta. These ejecta properties are potentially affected by dynamically-important feedback from $r$-process heating, which is usually not coupled to the hydrodynamics in post-merger simulations modeling the ejecta launching and expansion. The multi-messenger detection of GW170817 showed the importance of producing reliable ejecta predictions, to maximize the diagnostic potential of future events. In this paper, we develop a prescription for including $r$-process heating as a source term in the hydrodynamic equations. This prescription depends on local fluid properties and on the $Y_{e}$ history as recorded by dedicated tracer particles, which exchange information with the grid using the Cloud-in-Cell method. The method is implemented in long-term viscous hydrodynamic simulations of accretion disk outflows to investigate its feedback on ejecta properties. We find that $r$-process heating can increase the unbound disk ejecta mass by $\sim 10%$ relative to a baseline case that only considers alpha particle recombination. Nuclear heating also enhances the radial velocity of the ejecta with $Y_e < 0.25$ by up to a factor of two, while concurrently suppressing marginally-bound convective ejecta.

💡 Research Summary

**
This paper addresses a long‑standing omission in post‑merger simulations of neutron‑star (NS) mergers: the dynamical feedback of r‑process nuclear heating on the outflows from the remnant accretion disk. While most existing long‑term viscous (or magnetohydrodynamic) simulations include the energy release from α‑particle recombination and, to a lesser extent, from heavy nuclei in nuclear statistical equilibrium (NSE), they neglect the additional heating that occurs when the composition evolves away from NSE during the r‑process. The authors develop a practical prescription to incorporate this heating as a source term in the Eulerian hydrodynamic equations.

Methodology
The heating rates are taken from a library of one‑zone r‑process nucleosynthesis calculations (Lippuner & Roberts 2015) performed on parametrized expansion trajectories with entropy s = 30 k_B per nucleon and expansion timescale τ = 10 ms. Fifty initial electron‑fraction values (Y_e = 0.01–0.50 at 8 GK) are sampled, yielding a two‑dimensional table (\dot\epsilon_{\rm nuc}(T, Y_{e,T6})) where (Y_{e,T6}) denotes the electron fraction when the fluid first reaches 6 GK. Half of the released nuclear energy is assumed to escape as neutrinos; the remaining half is deposited as thermal energy.

Because the hydrodynamics is solved on an Eulerian grid while the heating information is naturally attached to Lagrangian fluid elements, the authors introduce two sets of tracer particles. “Passive tracers” follow the flow for diagnostic purposes, while a much larger set of “memory tracers” (≈ 10⁶ particles) are uniformly seeded throughout the initial torus to record the local Y_e history. At each timestep the grid temperature and Y_e are interpolated onto the tracers. When a tracer first passes through 6 GK, its Y_e is stored as (Y_{e,T6}). During the simulation, if the local temperature drops below a threshold (T_{\rm heat}=4) GK (the authors also test 6 GK), the heating rate for that tracer is obtained by (1) bracketing the stored (Y_{e,T6}) between the two nearest tabulated values, (2) performing a log‑linear interpolation in temperature at each bracketing Y_e, and (3) linearly interpolating between the two resulting rates. The final specific heating rate is then mapped back onto the grid using a Cloud‑in‑Cell (CIC) scheme. Inflowing material (negative radial velocity) receives the same rate with opposite sign, effectively providing a cooling term and preventing double counting when a tracer traverses the same region multiple times.

Simulation Setup
The authors use the FLASH code (v3.2) in 2‑D axisymmetric spherical coordinates, covering radii from 8.8 km to 8.8 × 10⁵ km with 640 logarithmic radial zones and 112 angular zones uniformly spaced in (\cos\theta). The central object is a pseudo‑Newtonian black hole of mass 2.65 M⊙ and spin a = 0.8, representing the remnant of GW170817. The initial torus has mass 0.1 M⊙, constant specific angular momentum, entropy 8 k_B per baryon, and uniform Y_e = 0.1. An α‑viscosity prescription (α = 0.03 or 0.06) drives angular‑momentum transport. The Helmholtz equation of state, a three‑species neutrino leakage scheme, and α‑particle recombination energy are included as in previous works (Fernández & Metzger 2013, etc.).

Results

1. Ejecta Mass: Including r‑process heating raises the total unbound disk‑wind mass by ≈ 10 % compared with the baseline (α‑recombination only). The extra heating inflates the disk atmosphere, allowing more material to become unbound.
2. Velocity: For low‑Y_e ejecta (Y_e < 0.25) the radial velocity can increase by up to a factor of two. The effect is weaker for higher Y_e because the heating rate declines more rapidly with temperature for neutron‑rich compositions.
3. Convective Ejecta: Material that would otherwise be marginally bound and ejected via convective motions is partially suppressed; the additional pressure from nuclear heating keeps it bound, shifting the outflow composition toward more robust, high‑velocity wind components.
4. Threshold Sensitivity: Raising the heating activation temperature to 6 GK reduces both the mass increase and velocity boost by roughly 30 %, indicating that the exact temperature at which r‑process heating becomes effective is a non‑negligible source of uncertainty.

Comparison with Prior Work
Earlier studies (Just et al. 2015; Wu et al. 2016) employed temperature‑only heating prescriptions and reported divergent impacts on ejecta properties. More recent works (Klion et al. 2022; Kawaguchi et al. 2021) used time‑dependent, composition‑independent parametrizations, finding modest velocity enhancements but little mass change. The present study improves upon these by explicitly tracking the Y_e history of each fluid element, thereby capturing the strong composition dependence of the heating rate. The results are broadly consistent with the trends reported by Just et al. (2025), who used a neural‑network surrogate for the heating rate, but the present work provides a transparent, physics‑based implementation.

Limitations and Future Directions

• The heating table is built from a single entropy (30 k_B) and expansion timescale (10 ms); real disk outflows exhibit a broader range of thermodynamic histories. Extending the library to multiple entropies and τ values would improve fidelity.
• Assuming a fixed 50 % neutrino loss fraction is a simplification; a more sophisticated neutrino transport treatment could modify the net heating.
• The simulations are axisymmetric and use an effective α‑viscosity rather than fully resolved magnetorotational turbulence; three‑dimensional MHD simulations with on‑the‑fly nuclear networks would be the ultimate test.

Conclusions
By coupling a composition‑dependent r‑process heating prescription to long‑term viscous disk simulations, the authors demonstrate that nuclear heating beyond α‑recombination can measurably alter both the mass and velocity distribution of post‑merger disk winds. The effect is most pronounced for neutron‑rich ejecta, which are crucial for the synthesis of the heaviest r‑process nuclei and for shaping the early blue kilonova component. The methodology—memory tracers, CIC mapping, and a two‑dimensional heating table—offers a computationally tractable yet physically motivated framework that can be incorporated into future multi‑physics merger simulations, thereby improving the reliability of kilonova light‑curve predictions and nucleosynthetic yield estimates.