© The Authors 2026

1. Introduction

Binary neutron star mergers (BNSMs) are one of the main sources for the production of heavy elements via the rapid neutron capture (r-process) nucleosynthesis (Eichler et al. 1989; Li & Paczynski 1998; Metzger et al. 2010; Cowan et al. 2021; Perego et al. 2021). The first multi-messenger observation of a BSNM was triggered by the gravitational wave (GW) signal GW170817, and associated with the short gamma-ray burst (GRB) GRB 170817A and the kilonova AT2017gfo, i.e., the electromagnetic (EM) counterpart powered by the decay of the freshly synthesized, unstable nuclei (Abbott et al. 2017b,a,c; Coulter et al. 2017; Drout 2017; Savchenko et al. 2017). Kilonova candidates have also been identified in association with other short GRBs (Berger et al. 2013; Gompertz et al. 2018; Rossi et al. 2020) and hybrid or long GRBs (Rastinejad et al. 2022; Troja et al. 2022).

The properties of the matter ejected from a BNSM (mass, temperature, velocity, and electron fraction) are determined by the merger and post-merger dynamics, neutrino interactions, and the nuclear equation of state (EOS) (Rosswog et al. 2014; Radice et al. 2018; Perego et al. 2019; Foucart et al. 2020; Bernuzzi 2020; Fujibayashi et al. 2020; Just et al. 2022; Kiuchi et al. 2023; Lund et al. 2024; Cheong et al. 2025). The nuclear reactions burning in the expanding material depend on these initial conditions and on the thermodynamic evolution of the outflow. The latter is dominated by its hydrodynamic expansion and affected by the released nuclear energy and radiative transport effects (Korobkin et al. 2012; Lippuner & Roberts 2015; Cowan et al. 2021; Fryer et al. 2024; Magistrelli et al. 2024). The final EM (optical, UV, near-infrared, and γ) emission is dictated by the thermal energy of the ejecta, the nuclear power released in the reacting material, the way the products of the nuclear reactions thermalize with the surrounding matter, the ejecta morphology, and the atomic opacities in the outflow (Barnes et al. 2016; Kasen et al. 2017; Wollaeger et al. 2018; Kasen & Barnes 2019; Hotokezaka & Nakar 2020; Tanaka et al. 2020; Korobkin et al. 2020; Fontes et al. 2020; Korobkin et al. 2021; Gillanders et al. 2022; Fontes et al. 2023; Collins et al. 2023; Shingles et al. 2023; Pognan et al. 2024; Sneppen et al. 2024). To understand the nucleosynthesis in BNSM ejecta, how they can participate in the galactic enrichment, and the associated EM counterparts, it is thus necessary to model all of these complex aspects of the merger and post-merger dynamics.

General-relativistic (magneto)-hydrodynamics (GRMHD) simulations implementing candidate nuclear EOSs can follow a binary system from the inspiral phase until hundreds of milliseconds after a merger, collecting information on the thermodynamic properties of the outflows (e.g., Radice et al. 2016; Bauswein et al. 2020; Nedora et al. 2021; Combi & Siegel 2023; Kiuchi et al. 2024; Fields et al. 2025). New moment-based neutrino-transport schemes and Monte Carlo methods have improved our understanding of the role of neutrinos in driving mass ejection (in particular for long-lived remnants) and setting the electron fraction, and thus initial composition, of these outflows (Foucart et al. 2016; Radice et al. 2022; Zappa et al. 2023; Kawaguchi et al. 2025). Usually ejecta properties extracted from numerical relativity (NR) simulations are transferred to other codes for their long-term (approximately days) evolution with the use of tracers or some suitable mapping (Korobkin et al. 2012; Martin et al. 2015; Barnes et al. 2016; Radice et al. 2018; Perego et al. 2022; Curtis et al. 2024; Groenewegen et al. 2025).

As the material expands and cools, the temperature eventually drops below the nuclear statical equilibrium (NSE) threshold. Beyond this point, an additional source of uncertainty in modeling r-process nucleosynthesis arises from the poorly constrained nuclear properties of isotopes far from stability. In particular, simulating ejecta form BNSMs requires knowledge of the properties of nuclei near the neutron drip line. Despite significant theoretical and experimental effort, key quantities such as nuclear masses, β-decay rates, neutron-capture cross sections, and fission yields remain highly uncertain, yet they strongly impact nucleosynthesis pathways and influence EM observables (Mumpower et al. 2016; Barnes et al. 2021; Zhu et al. 2022; Martinez-Pinedo & Langanke 2023; Mumpower et al. 2024; Martinet & Goriely 2025).

Most of the current ejecta models treat nuclear physics and radiation-hydrodynamic separately. For example, Rosswog et al. (2014), Kawaguchi et al. (2021, 2024), and Wu et al. (2022) input precomputed, analytical nuclear powers fits in the energy equation. In Korobkin et al. (2012), Radice et al. (2018), Perego et al. (2022), Gillanders et al. (2022), Just et al. (2022), and Curtis et al. (2024), the nucleosynthesis was calculated via a post-processing method by prescribing a simplified homologous expansion and feeding it to a Nuclear reaction Network (NN). Similar simplifications have been used to perform full Monte Carlo radiative-transport simulations to calculate kilonova light curves and spectra (Collins et al. 2023; Shingles et al. 2023; Pognan et al. 2024). These methods neglect the complex interplay between fluid dynamics, nucleosynthesis, and radiative transport, and they fail to capture the nonlinear feedback of radioactive heating and composition-dependent radiation-matter interactions on the outflow dynamics. To address this limitation and study the interplay between ejecta dynamics and nucleosynthesis, Magistrelli et al. (2024) performed the first analysis of BNSM ejecta evolution that couples ray-by-ray radiation-hydrodynamic expansion with an in situ NN. This more self-consistent modeling revealed qualitative differences in the predicted kilonova emission and nucleosynthesis. The same framework was applied in Bernuzzi et al. (2025) and Jacobi et al. (2026) to investigate outflows from asymmetric BNSMs.

In this work we assess the impact of in situ nuclear-burning models and atomic opacities in radiation-hydrodynamics BNSM ejecta simulations. In order to explore efficiently different models, we focused on a 2D ray-by-ray setup where radiation hydrodynamics was treated in the diffusion limit. In Sect. 2 we describe our frequency-dependent radiation-hydrodynamics setup and its coupling with the in situ NN, which updated the kNECnn code (Magistrelli et al. 2024; see also Morozova et al. 2015; Wu et al. 2022). The code implements simple, commonly used analytical prescriptions for the nuclear burning and opacities as well as different models for composition-informed thermalization and opacities enabled by the NN coupling. We also discuss the extraction of initial ejecta profiles from NR simulations and describe how we incorporated the effects of a jet depositing energy into the polar regions of the ejecta. Further details on the frequency-dependent radiative transport equations and on the initialization of the NN can be found in Appendices A and B. In Sect. 3 we analyze the impact of each physical improvement introduced here relative to previous works. We examine how the inclusion of an in situ NN alters nucleosynthesis and light-curve predictions, and we asses the improvements enabled by the resulting composition-dependent thermalization and opacity models. We also studied the impact of the additional energy deposited by a polar jet. We summarize our findings in Sect. 4.

2. Method

In this work we use kNECnn, a 2D ray-by-ray Lagrangian radiation-hydrodynamic code with in situ NN, composition-dependent thermalization and opacities, and a frequency-dependent radiation transport scheme. We based our development on the original SNEC code from Morozova et al. (2015) and its updated versions from Wu et al. (2022), Magistrelli et al. (2024), extending the treatment of nuclear burning, thermalization, radiation transport, and opacities.

2.1. Ray-by-ray radiation-hydrodynamics

The ejecta are evolved by solving the system of Lagrangian radiation-hydrodynamic equations presented in SNEC (Morozova et al. 2015). Nonspherical effects in the dynamics are approximately incorporated under the axisymmetry and ray-by-ray assumptions. The polar angle dependence of the ejecta properties is taken into account by mapping each angular section into an effective 1D, spherically symmetric problem. The total mass of each section is scaled by the factor λ_θ = 4π/ΔΩ, where ΔΩ is the solid angle subtended by the angular section. The sections are evolved independently, and the hydrodynamic variables and nucleosynthesis results are then recombined via mass weighted averages. The mapping procedure assures that the intensive quantities are always preserved (Magistrelli et al. 2024). Note that non-radial flows of matter and radiation are neglected, as well as higher dimensional fluid instabilities (e.g., Kelvin-Helmholtz), convection and the possibility of shells surpassing each other. These assumptions could in principle lead to overestimating the pressure work exchanged between the radial shells and introduce artificial shocks in the system.

In the 1D effective simulation of each angular section, the energy equation in spherical symmetry reads

$$\begin{array}{c}\hfill \frac{\partial ϵ}{\partial t}=\frac{P}{\rho }\frac{\partial ln\rho }{\partial t}-4\pi r^{2}Q\frac{\partial v}{\partial m}-\frac{\partial L}{\partial m}+{\dot{ϵ}}_{\mathrm{nucl}},\end{array}$$(1)

where ϵ is the specific internal energy, t is time, and P, ρ, r, v and L are the pressure, density, radial position, radial velocity, and luminosity of the fluid element. The Lagrangian coordinate is m, and Q is the von Neumann-Richtmyer artificial viscosity (Von Neumann & Richtmyer 1950). The thermalized nuclear heating rate ${\dot{ϵ}}_{\mathrm{nucl}}$ accounts for the local production and partial thermalization of the energy released by all possible nuclear reactions occurring in the expanding material. As already discussed in Magistrelli et al. (2024), the EOS closing the set of hydrodynamic equations combines of the tabulated Helmholtz EOS (Timmes & Swesty 2000; Lippuner & Roberts 2017) at high temperatures with the Paczynski EOS (Paczynski 1986; Morozova et al. 2015; Wu et al. 2022) at lower temperatures, where the contribution of positrons is negligible.

The complete composition tracking implemented in kNECnn allows us to calculate the optical opacity on the fly accounting for the evolution of the density, temperature and composition profiles. In our most sophisticated opacity model, presented in Sect. 2.5, we make use of frequency-dependent opacities to calculate the luminosity in Eq. (1), the position of the photospheres, and the observed kilonova fluxes. The luminosity appearing on the right-hand side of Eq. 1 is obtained as

$$\begin{array}{c}\hfill L=4\pi r^2{\displaystyle {\int }_0^{\infty }}{F}_{\nu }\mathrm{d}\nu ,\end{array}$$(2)

where the radiative flux for each shell is calculated under the assumption of local thermodynamical equilibrium (LTE) as

$$\begin{array}{c}\hfill F_{\nu }=-{(4\pi r)}^2\frac{{\lambda }_{\nu }}{3{\overline{\kappa }}_{\nu }}\frac{\partial B_{\nu }}{\partial T}\frac{\partial T}{\partial m}.\end{array}$$(3)

Here, B_ν is the Planck spectrum, ${\overline{\kappa }}_{\nu }$ the composition-averaged opacity (see Sect. 2.5 for further details), and λ_ν the flux-limiter

$$\begin{array}{c}\hfill {\lambda }_{\nu }=\frac{6+3R_{\nu }}{6+3R_{\nu }+R_{\nu }^2},\end{array}$$(4)

with

$$\begin{array}{c}\hfill R_{\nu }=\frac{4\pi r^2}{{\overline{\kappa }}_{\nu }}\frac{\partial log{B}_{\nu }}{\partial T}\left|\frac{\partial T}{\partial m}\right|,\end{array}$$(5)

which is proportional to the one defined in Levermore & Pomraning (1981, see our Appendix A for further details).

The numerical solution of Eq. (1) requires an implicit solve involving the EOS and the temperature. This is performed using a Newton-Raphson iteration on the temperature, nested inside a fixed-point iteration for the internal energy. The frequency-dependent generalization follows the same algorithm, but each Newton–Raphson step additionally requires evaluating temperature derivatives of the blackbody function as well as computing integrals over all frequency groups.

To estimate the kilonova light curves, we first calculate the position r_νph of each frequency-dependent photosphere by imposing

$$\begin{array}{c}\hfill {\displaystyle {\int }_{r_{{\nu }_{\mathrm{ph}}}}^{\infty }}{\tau }_{\nu }(r)\mathrm{d}r=\frac{2}{3},\end{array}$$(6)

where ${\tau }_{\nu }={\overline{\kappa }}_{\nu }\rho$ is the optical depth. For thermal emission, the flux emitted by a spherically symmetric source is given by f_ν^S = πB_ν, where B_ν is the blackbody function. The contribution from the photosphere to the flux measured by an observer at a distance D is given by f_ν = f_ν^Sr_νph²/D². Using the Stefan-Boltzmann law and the same arguments as in Wu et al. (2022) for the nuclear contributions from regions outside the photosphere, the observed kilonova fluxes can be expressed as

$$\begin{array}{c}\hfill f_{\nu }=\frac{1}{4\pi D^2}(4\pi r_{{\nu }_{\mathrm{ph}}}^2{F}_{\nu ,{\nu }_{\mathrm{ph}}=\nu }+{\displaystyle {\int }_{r_{{\nu }_{\mathrm{ph}}}}^{r_{\max }}}\frac{\pi {\dot{ϵ}}_{\mathrm{nucl},m}}{\sigma T_m^4}{B}_{\nu }(T_m)dm),\end{array}$$(7)

where F_ν, νph is the flux at frequency ν at the photosphere for photons of frequency ν_ph, r_max is the outer boundary of the outflow, ${\dot{ϵ}}_{\mathrm{nucl},m}$ is the local nuclear heating rate, and σ is the Stefan-Boltzmann constant. The observed AB magnitudes are then calculated as in Wu et al. (2022) and recombined accounting for the viewing angle as in Martin et al. (2015) and Perego et al. (2017).

2.2. Nuclear network coupling

We evolve the matter composition and compute the specific energy deposition ${\dot{ϵ}}_{\mathrm{nucl}}$ self-consistently, by providing each fluid element in our simulation with an in situ NN. In the current version of the code, we use an implementation of the SkyNet NN (Lippuner & Roberts 2017) which we adapted to interact with the kilonova-version of SNEC from Wu et al. (2022). The coupling infrastructure is general and in principle able to host other NNs. The NN used in this paper includes 7836 isotopes up to ³³⁷Cn and uses the JINA REACLIB (Cyburt et al. 2010) database and the same setup as in Lippuner & Roberts (2015), Perego et al. (2022), Magistrelli et al. (2024).

Within every hydrodynamic time-step Δt, each fluid element independently evolves the associated NN, which automatically defines its own sub-steps based on the nuclear timescales. We call t₀ the initial time of the hydrodynamical time-step. During the NN evolution, the density is prescribed by a log-interpolation between ρ(t₀) and ρ(t₀ + Δt). The latter is calculated from ρ(t₀) and the velocity of the shell boundaries at t = t₀, updated via the momentum transport equation (for which no nuclear input is required). During the timestep, we assume a constant T = T(t₀) for the NN calculations (no self-heating). The temperature will then be updated by the second part of the hydrodynamic step. We collect the energy released by nuclear reactions during each sub-step, and define the (non-thermalized) nuclear power ${\dot{q}}_{\mathrm{nucl}}(t_0)$ as the total energy produced, divided by Δt. This quantity will enter the energy equation, Eq. (1), after being thermalized into ${\dot{ϵ}}_{\mathrm{nucl}}(t_0)$ as described in Sect. 2.4. At the end of the NN calculations, we use the updated isotopic composition to get the mean molecular weight, which enters the EOS, and the opacity profile, needed to calculate the luminosity term in Eq. (1) and the observed fluxes. In other words, we use the abundances at t = t₀ + Δt as a representation of the composition during the whole hydrodynamic step. This procedure produces all the information needed to evolve Eq. (1), and thus to update temperature and velocity of the outflow.

Our ab initio NR simulations do not track nor evolve the detailed matter composition, but always rely on EOS tabulated for matter in nuclear statistical equilibrium (NSE). We must thus prescribe the isotopic composition at the beginning of the radiation-hydrodynamical simulation. At merger, the colliding matter from the two neutron stars typically reaches a high enough temperature (T ≳ 8 GK) to ensure hot NSE conditions. However, expansion usually causes the unbound material to drop out of NSE before reaching the extraction radius. In particular, many fluid elements recorded at the extraction radius have T < 5 GK in the initial ejecta profile (see Fig. 2). If tracers are modeled inside the NR simulation, or if they can be reconstructed in post-processing, one can study the evolution of these fluid elements starting from their reconstructed thermodynamics trajectories. In this case, one can post-process a NN on these trajectories and predict the composition of the fluid elements at the extraction radius properly accounting for the out-of-NSE transition. The operation is not self-consistent, as it does not include out-of-NSE effects in the original dynamic simulation. However, it gives a reliable estimate for the composition of the ejecta at the extraction radius that can be used to initialize kNECnn. If the necessary information to define tracers is missing, our code relies on the NSE assumption to estimate the initial isotopic composition. This is the case for the NR simulations considered in this work.

The simplest initialization method implemented in kNECnn, and already used in Magistrelli et al. (2024), assumes NSE at any temperature T₀ registered at the extraction radius. We refer to this procedure as cold-NSE initialization (CNSE) in the rest of this work. For each tracked isotope, the initial abundance is given by

$$\begin{array}{c}\hfill Y_i=\frac{n_i}{n_b}\propto Y_p^{Z_i}{Y}_n^{(A_i-Z_i)}\frac{{\rho }^{A_i-1}{e}^{{\mathrm{BE}}_i/k_{\mathrm{B}}{T}_0}}{T_0^{3(A_i-1)/2}},\end{array}$$(8)

where n_b is the baryon number density, n_i, A_i, Z_i, and BE_i are the number density, mass number, atomic number and binding energy of the i-th isotope, Y_p and Y_n are the abundances of free protons and neutrons, and k_B is the Boltzmann constant (Cowan et al. 2021; Perego et al. 2021). The initial isotopic mass fraction is given by X_i = A_iY_i. The mass and charge conservation, ∑_iA_iY_i = 1 and ∑_iZ_iY_i = Y_e, with Y_e electron fraction, constrain the NSE composition, which thus depends on ρ, T and Y_e. See Lippuner & Roberts (2017) for details about SkyNet’s implementation of the NSE condition.

For matter that experienced hot NSE during the NR simulation, but underwent NSE freeze-out before reaching the extraction radius, the Boltzmann term in Eq. (8) tends to overestimate the abundances of the iron-group elements. Moreover, the impact of nuclear reactions happening immediately after NSE drop-out is neglected. To correct these effects, we implement a backtracking-NSE initialization similar to the one described in Radice et al. (2016), Perego et al. (2022). Essentially, if T₀, ρ₀, Y_e, 0, s₀ are the properties at the extraction radius, we assume that the fluid element had expanded adiabatically from its NSE configuration at some T_NSE > T₀, ρ_NSE, Y_e, NSE, s_NSE. On the small timescales (≲10 ms) needed for the ejected matter to get to the extraction sphere, we expect the electron fraction and the entropy not to change significantly, so Y_e, NSE ≃ Y_e, 0, s_NSE ≃ s₀. We check a posteriori that the relative changes in the entropy and electron fraction are on the order of 5% or less during the pre-evolution phase in most part of the ejecta (see Fig. B.2). We finally estimate the density ρ_NSE of the fluid element at NSE freeze-out by using the EOS employed in the NR simulation and starting from T_NSE, Y_e, NSE = Y_e, 0, s_NSE = s₀. We calculate the correspondent NSE composition via Eq. (8), and then mimic the abundance evolution up to the extraction radius assuming homologous expansion. In particular, we parametrize the density evolution as in Eq. (1) of Lippuner & Roberts (2015), and estimate the expansion timescale τ_ex as in Radice et al. (2016)¹. To ensure that the evolved thermodynamic trajectory matches the physical conditions of the fluid element at the extraction radius, we also use an analytical prescription for the temperature history inspired by the adiabatic expansion law for an ideal gas, T(ρ) = T_NSE (ρ(t)/ρ_NSE)^{(γ − 1)}. The effective adiabatic factor γ, specific for each fluid element, is determined by γ = 1 + ln(ρ₀/ρ_NSE)/ln(T₀/T_NSE), which ensures T(ρ₀) = T₀. The composition at the end of the backtracking is used together with the originally extracted s to define the initial condition of the fluid element.

Low-entropy tidal ejecta (more efficiently produced for strongly asymmetric binaries) are expelled before the merger can raise the temperature, and are therefore expected to retain a composition close to the one set by the cold-NSE condition, for which T ≲ 1 MeV. In principle, this would suggest using the cold-NSE method. However, in an NR simulation, numerical artifacts (e.g., shock heating propagating inward from the surface of the compact objects) can boost the initial temperature of the neutron stars up to few MeV. In all the NR simulations analyzed in this paper, the neutron stars reach this temperature in ∼1 ms during the inspiral phase. The spurious high temperatures introduce small artifacts in the tracers calculations and in the initialized NSE composition, regardless of the method used. In Appendix B we investigate the systematic uncertainties propagating from the different initialization methods to the nucleosynthesis and kilonova results.

2.3. Nuclear power fits

To investigate the impact of the in situ NN, we need benchmark simulations using only ex situ calculations. For these simulations (labeled in what follows as Apr2) we calculate the nuclear power with an updated version of the analytical fits from Wu et al. (2022). At early times, t ≲ 0.1 days, we consider the fitting function

$$\begin{array}{c}\hfill {\dot{q}}_{\mathrm{nucl}}(t)=q_1{e}^{-t/\beta }+q_0{[\frac{1}{2}-\frac{1}{\pi }arctan\left(\frac{t-t_0}{\sigma }\right)]}^{\alpha },\end{array}$$(9)

where q₁, β, q₀, α, t₀, and σ are fitting parameters. The additional exponential term effectively generalizes the formula proposed by Korobkin et al. (2012) to include initially mildly neutron-rich (0.4 ≲ Y_e, 0 < 0.5) and proton-rich (Y_e, 0 > 0.5) fluid elements. These trajectories, which do not produce r-process elements, typically exhibit an early fast drop in the associated nuclear power, followed by a return to the approximate arctan behavior (although usually shifted down by several orders of magnitude compared to an equivalent more neutron-rich fluid element). At later times, t ≳ 0.1 days, we use the power law fit

$$\begin{array}{c}\hfill {\dot{q}}_{\mathrm{nucl}}(t)=q_0\prime t^{-\alpha \prime },\end{array}$$(10)

where q₀′ and α′ are additional fit parameters. As in Wu et al. (2022), the two nuclear powers are joined together by a log-scaled smoothing procedure applied between 10³ ≤ t [s] ≤ 4 × 10⁴ and centered on t ∼ 0.1 days. The two fits combined describe the nuclear power over times after merger between 0.1 s and 50 days.

The original fits from Wu et al. (2022) relied on the nuclear power grid presented in Perego et al. (2022). We perform our fits on the grid obtained in Chiesa et al. (2024), which improves on the original one presented in Perego et al. (2022) by initializing the NN at a higher NSE temperature, T_NSE = 8 GK, to better capture NSE drop-out. Additionally, specifically for this work, we have extended this grid by increasing the initial electron fraction, entropy, and expansion timescales ranges to 0.01 < Y_e, 0 < 0.6, 1.4925 < s [k_B baryon⁻¹]< 300, 0.5 < τ [ms]< 200, respectively. Note that Wu et al. (2022) ceil the initial electron fraction to Y_e, 0 = 0.48 and use the table boundaries for higher values. The NN is set up as described in Sect. 2.2, consistent with the rest of the calculations in this work.

In Fig. 1 we compare the nuclear power calculated with SkyNet against the predictions from the original fits of Wu et al. (2022) and our updated version. The old fits fail to capture the very early-time energy production of the mildly neutron-rich and proton-rich matter, as well as the late-time nuclear power for proton-rich ejecta, which constitute most of the outflow for angles 30 ≲ θ [degrees]≲60, by several orders of magnitude. The new fits only perform worse than the old ones in localized time intervals (e.g., 1 < t [s]< 10³ for Y_e = 0.46, 0.48). The discrepancy with the original SkyNet calculations remains within a factor of order unity, assuring an overall improvement of the original fits of Wu et al. (2022). While a precise estimate of the nuclear power is not needed at early times, when the dynamics are dominated by the hydrodynamic explosion, it becomes crucial at late times to obtain reliable predictions of kilonova light curves. In 3D ejecta simulations, where nuclear heating can affect the angular dynamics of the ejecta for timescales t ≳ 1 s, more accurate nuclear power fits might be necessary to capture the correct evolution of the outflow.

Fig. 1. Nuclear power as a function of time. All the trajectories have initial entropy s0 ≃ 11 kB baryon−1, expansion timescale τ0 ≃ 8.6 ms, and electron fraction Ye ∈ [0.02, 0.1, 0.2, 0.3, 0.4, 0.46, 0.48, 0.54] indicated by the color map. The vertical dashed lines mark the transition between the early- and late-times fits. Top panel: SkyNet calculations (solid lines) and estimate from the updated (dashed lines) and old (dotted lines) nuclear power fits, respectively from this work and Wu et al. (2022). Bottom panel: Relative differences between the results from SkyNet and from the new (solid lines) and previous (dashed lines) version of the fits. The dashed horizontal lines mark differences of one order of magnitude.

2.4. Thermalization

In the simple thermalization scheme implemented in the kilonova version of SNEC by Wu et al. (2022), the thermalized heating rate is approximated by ${\dot{ϵ}}_{\mathrm{nucl}}=f_{\mathrm{th}}^{\mathtt{kNEC}}{\dot{q}}_{\mathrm{nucl}}$, where ${\dot{q}}_{\mathrm{nucl}}$ is the total energy released by nuclear reactions, and $f_{\mathrm{th}}^{\mathtt{kNEC}}$ is a constant and uniform thermalization coefficient. The standard value $f_{\mathrm{th}}^{\mathtt{kNEC}}=0.5$ is adopted to represent the average behavior of all particles emitted by nuclear decays, ranging from easily escaping neutrinos to heavy, efficiently thermalizing, charged particles (fission fragments). In what follows, we refer to this thermalization scheme as ST05.

In kNECnn, we implement time- and composition-dependent thermalization coefficients separately for electrons and positrons from β^± decays, and for α particles and γ rays emitted from all possible decays, including the m_e ≃ 511 keV photons from electron-positron annihilations. We explicitly remove the energy carried away by (electron) neutrinos and antineutrinos emitted by β^± decays from the thermalized heating rate. We denote this model as DT08 in the rest of the paper.

We calculate the thermalization efficiency $f_{\mathrm{th}}^{\alpha ,\beta }(t)$ for α particles and electrons with the analytical prescriptions of Kasen & Barnes (2019). We assume that the same expression derived for electrons also holds for positrons. The decay energy injected in these particles is calculated as ${\dot{q}}_k(t)={\displaystyle \sum _j}⟨{\dot{q}}_j^k(t)⟩$, where $⟨{\dot{q}}_j^k(t)⟩$ is the average energy injected by the j-th decaying isotope in particles of the radiation type k = α, e⁻, e⁺. We calculate this quantity as

$$\begin{array}{c}\hfill ⟨{\dot{q}}_j^k(t)⟩=\frac{⟨E_j^k⟩Y_j(t)}{m_p{\tau }_j},\end{array}$$(11)

where m_p is the proton mass, τ_j is the average lifetime of the j-th isotope, and ⟨E_j^k⟩ is the mean energy released in particles of type k after a radioactive decay. The average lifetime accounts for all possible decays of the isotope, and the mean released energy averages over all possible decays accounting for their branching ratios. Both these quantities are taken from the ENDF/B-VIII.0 database (Brown et al. 2018). We also include the free-neutron decay information from the IAEA nuclear data service website². The same expressions for $⟨{\dot{q}}_j^k(t)⟩$ and Eq. (11) generalize for k = ν, and prompt γ-rays.

kNECnn internally checks that the total rate of emitted energy $⟨\dot{q}⟩={\displaystyle \sum _{k,j}}⟨{\dot{q}}_j^k(t)⟩$ does not exceed the nuclear power ${\dot{q}}_{\mathrm{nucl}}$ calculated by the NN, and rescales all the ${\dot{q}}_k={\displaystyle \sum _j}⟨{\dot{q}}_j^k(t)⟩$ by the factor ${\dot{q}}_{\mathrm{nucl}}/⟨\dot{q}⟩$ in case this happens. We then calculate the power emitted in γ-rays from positron-electron pair annihilations (not included in the ${\dot{q}}_{\mathrm{nucl}}$ from the NN), as $⟨{\dot{q}}_j^{\mathrm{pair}}(t)⟩=2m_e⟨n_{e^{+}}⟩$, where ⟨n_e⁺⟩ is the average number of emitted positrons reported in the database.

The thermalization factor f_j^γ(t) for prompt γ-rays is calculated, as in Hotokezaka & Nakar (2020), Combi & Siegel (2023), starting from the detailed composition of the ejecta and the same effective opacity tables from the NIST-XCOM catalog (Berger et al. 2010) used in Barnes et al. (2016). The thermalization efficiency associated with the γ particles emitted by nuclei of the j-th isotope through radioactive decays is given by

$$\begin{array}{c}\hfill f_j^{\gamma }(t)=1-exp[-{\kappa }_j^{\gamma }{\mathrm{\Sigma }}_m(t)],\end{array}$$(12)

where κ_j^γ is the absorptive opacity and Σ_m(t) is the mass-weighted column density of the ejecta. In spherical symmetry,

$$\begin{array}{c}\hfill {\mathrm{\Sigma }}_m(t)={\displaystyle {\int }_{m_0}^{m_{\max }}}\frac{\mathrm{d}m}{M_{\mathrm{ej}}}{\displaystyle {\int }_m^{m_{\max }}}\frac{\mathrm{d}m\prime }{4\pi r^2(m\prime )},\end{array}$$(13)

where M_ej and m₀ are the total mass of the ejecta and the mass of the remnant (which coincides with the Lagrangian coordinate of the first simulated shell) and m_max = m₀ + M_ej. The absorptive opacity can be written as

$$\begin{array}{c}\hfill {\kappa }_j^{\gamma }={\displaystyle \sum _k}{h}_{k,j}^{\gamma }{\kappa }_{\mathrm{eff}}^{\gamma }(E_{k,j}),\end{array}$$(14)

where h_k, j^γ is the fraction of energy emitted in the k-th γ-ray, of energy E_k, j, by a decay of the j-th isotope (Hotokezaka & Nakar 2020), and κ_eff^γ(E_k, j) is the effective absorptive opacity. We read this information again from the ENDF/B-VIII.0 database (Brown et al. 2018), where the injection energy distributions are already averaged, for each isotope, over all accessible decays, accounting for their branching ratios. The effective opacity is given by the mass-weighted average

$$\begin{array}{c}\hfill {\kappa }_{\mathrm{eff}}^{\gamma }(E_{k,j})={\displaystyle \sum _Z}{X}_Z{\kappa }_{\mathrm{eff},Z}^{\gamma }(E_{k,j}),\end{array}$$(15)

where Z runs over all possible atomic numbers, X_Z = ∑_{s : Zs = Z}X_s, with s running through all the isotopes considered in the NN, and κ_eff, Z^γ(E_k, j) is the effective opacity for a stopping material made only of isotopes with atomic number Z. The latter is included in the NIST-XCOM catalog. Since this database only contains information up to Z = 100, we fix κ_{eff, Z > 100}^γ(E_k, j) = κ₀ and use Eq. (15) to rewrite Eq. (14) as

$$\begin{array}{c}\hfill {\kappa }_j^{\gamma }={\displaystyle \sum _{k,Z\le 100}}{X}_Z{h}_{k,j}^{\gamma }{\kappa }_{\mathrm{eff},Z}^{\gamma }(E_{k,j})+X_{\mathrm{rem}}{\kappa }_0,\end{array}$$(16)

where X_rem = 1 − ∑_{Z ≤ 100}X_Z is the cumulative mass fraction of the elements excluded from the database. Assuming that the energy released in γ-rays is immediately thermalized, one gets

$$\begin{array}{c}\hfill {\dot{ϵ}}_{\mathrm{nucl}}^{\gamma }={\displaystyle \sum _j}{f}_j^{\gamma }(t)⟨{\dot{q}}_j^{\gamma }(t)⟩,\end{array}$$(17)

where $⟨{\dot{q}}_j^{\gamma }(t)⟩$ is given by Eq. (11) with k = γ. For the simulations performed in this work, we fix κ₀ = 0.1 cm² g⁻¹. Equations (12)-(17) also apply for γ-rays for the positron-electrons pair annihilations, for which $h_{k,j}^{\mathrm{pair}}=\delta (E_{k,j}-m_e)$ and ⟨E_j^pair⟩ = 2m_eBR_β⁺, with m_e electron mass and BR_β⁺ branching ratio of the β⁺ decay.

For the remaining nuclear power, associated with the recoil energy of daughter nuclei, fission fragments, (soft) X-rays, proton and neutron emissions, and delayed electron (e.g., Auger) emission, we assume a constant and uniform thermalization factor $f_{\mathrm{th}}^{\mathrm{oth}}$. The heavy charged particles are expected to thermalize efficiently at all times, leaving only to X-rays, injected nucleons and delayed electrons with a nontrivial behavior. In our simulations, we choose, as a standard value, $f_{\mathrm{th}}^{\mathrm{oth}}=0.8>f_{\mathrm{th}}^{\mathtt{kNEC}}$, as we are now explicitly removing the neutrino contributions. The total thermalized heating rate is finally given by

$$\begin{array}{c}\hfill {\dot{ϵ}}_{\mathrm{nucl}}={\dot{ϵ}}_{\mathrm{nucl}}^{\alpha }+{\dot{ϵ}}_{\mathrm{nucl}}^{\beta }+{\dot{ϵ}}_{\mathrm{nucl}}^{\gamma }+{\dot{ϵ}}_{\mathrm{nucl}}^{\mathrm{pair}}+{\dot{ϵ}}_{\mathrm{nucl}}^{\mathrm{oth}},\end{array}$$(18)

where ${\dot{ϵ}}_{\mathrm{nucl}}^k=f_{\mathrm{th}}^k{\dot{q}}_k$ for k ∈ [α, γ, pair, oth], ${\dot{ϵ}}_{\mathrm{nucl}}^{\beta }=f_{\mathrm{th}}^{\beta }({\dot{q}}_{e^{-}}+{\dot{q}}_{e^{+}})$, and ${\dot{q}}_{\mathrm{oth}}={\dot{q}}_{\mathrm{nucl}}-{\dot{q}}_{\alpha }-{\dot{q}}_{e^{-}}-{\dot{q}}_{e^{+}}-{\dot{q}}_{\gamma }-{\dot{q}}_{\nu }$. The e⁺ − e⁻ pair annihilation photons are not included in the heating rate ${\dot{q}}_{\mathrm{nucl}}$ from the NN, and thus their contribute must not be removed from ${\dot{q}}_{\mathrm{oth}}$.

2.5. Opacities

In kNECnn, we implement several approaches to calculate the optical opacity. The simplest model (named W[22]) uses gray opacities. By plugging the standard definition of the Rosseland mean opacity into Eq. (2), the radiative luminosity is

$$\begin{array}{c}\hfill L_R=-{(4\pi r^2)}^2\frac{\lambda ac}{3{\kappa }_R}\frac{\partial T^4}{\partial m},\end{array}$$(19)

where λ is an equivalent flux limiter as the one used in Morozova et al. (2015) and Wu et al. (2022, see our Appendix A for further details). The simplest model is the same gray opacity prescription as in Wu et al. (2022), where κ_R = κ(Y_e, 0) has a constant-in-time profile depending only on the initial electron fraction. The general physical idea behind this assumption is that, within each fluid element, a lower Y_e, 0 will lead to the production of more r-process elements, resulting in a higher opacity. This is a strong approximation, as it does not use information about the actual ejecta composition or its dynamics. This opacity model was employed in Magistrelli et al. (2024), Bernuzzi et al. (2025), and Jacobi et al. (2026).

As a first, direct improvement of the analytical expression from Wu et al. (2022), we implemented in kNECnn the opacity model of Just et al. (2022). In this model (referred to in what follows as J[22]) the opacity is a time-dependent analytical function of the cumulative lanthanide mass fractions and the gas temperature. This models is based on the light curves calculated in Kasen et al. (2017), which rely on radiative-transfer equations and atomic-physics calculations of the opacities of lanthanides and iron-group elements. The prescription of Just et al. (2022) approximately reproduces the qualitative behavior of the kilonova from the original model. It directly connects the production of lanthanides with a boost in the effective opacity, and it approximately includes the effect of electronic recombination at small temperatures (T ≲ 2000 K).

For models based on the full tracking of the matter composition, we rely on the line-binned opacity calculations presented in Fontes et al. (2020, 2023), and available online via the NIST-LANL Opacity Database (Ralchenko et al. 2025). In these publications, the population distribution of the ionization states is computed for a range of temperatures and densities using the Saha equation. Detailed ab initio atomic-structure calculations are then employed to compute the opacity as a function of discretized frequency intervals. No assumptions are needed about the dynamical expansion. kNECnn implements data for Cr, Se, Br, Zr, Pd, Te, Ce, Nd, Sm, and U. In this work we focus on the effect of using atomic data for lanthanides and actinides, and use only the resulting opacities for Nd and U to effectively define a composition-averaged opacity of the ejecta. We discretize the frequency domain ν ∈ [3 × 10⁹,  3.6 × 10¹⁹] Hz in n_gr = 300 groups. We use a piecewise logarithmic frequency grid with coarse spacing at low (ν ≲ 2.42 × 10¹² Hz) and high (ν ≳ 7.25 × 10¹⁵ Hz) energies, and higher resolution in the central energy region to better resolve spectral features. We then calculate the group opacity from the line-binned results as prescribed in (Fontes et al. 2020). The database includes data for 10⁻²⁰ < ρ[g/cm³]< 10⁻⁴ and 116 ≲ T[K]≲5.8 × 10⁴. For densities or temperatures above the limits, we use the analytical model from Just et al. (2022). Close to the upper boundaries of the table, we use an activation function to smoothly switch between the two models. For densities or temperatures below the limits, we use the opacity data at the lower edge of the table. We use Nd and U data as representatives for all the lanthanides or actinides, respectively. All the remaining elements are assigned an opacity κ_oth = 0.2 k_T cm²g⁻¹, with k_T = 1 if T > 2000 K and k_T = (T[K]/2000)⁵ for T < 2000 K, consistent with the no-lanthanides case of the Just et al. (2022) approximation. To get the composition-averaged opacity ${\overline{\kappa }}_{\nu }$ of Eq. (3), we first estimate the opacity κ_ν, i(ρ, T) of each representative element by interpolating its tabulated values in the ρ − T plane. Then, for each frequency group, we average the contributions of the representative elements, weighting each element by the total mass fractions of the species it represents. The composition-averaged opacity is given by

$$\begin{array}{c}\hfill {\overline{\kappa }}_{\nu }(\rho ,T|t)={\displaystyle \sum _i}{\displaystyle \sum _j}{X}_{i,j}(t){\kappa }_{\nu ,i}(\rho ,T|t)+X_{\mathrm{oth}}(t){\kappa }_{\mathrm{oth}}(T),\end{array}$$(20)

where X_oth = 1 − ∑_i, jX_i, j. In the sum, the index i runs over the representative elements, while j runs over all elements that are mapped to the i-th representative.

As a first gray model (LANL-R), we use again the Rosseland luminosity L_R defined in Eq. (19), but we calculate κ_R directly from the line-binned database as

$$\begin{array}{c}\hfill {\kappa }_R={\left[\frac{{\sum }_{\nu }{\overline{\kappa }}_{\nu }^{-1}{\partial }_T{B}_{\nu }\mathrm{\Delta }\nu }{{\sum }_{\nu }{\partial }_T{B}_{\nu }\mathrm{\Delta }\nu }\right]}^{-1}.\end{array}$$(21)

We also implement an analogous model (LANL-P) with the Planck gray opacity

$$\begin{array}{c}\hfill {\kappa }_P=\frac{{\sum }_{\nu }{\overline{\kappa }}_{\nu }{B}_{\nu }(T)\mathrm{\Delta }\nu }{{\sum }_{\nu }{B}_{\nu }(T)\mathrm{\Delta }\nu }.\end{array}$$(22)

The Rosseland approximation, weighting more the transparent frequencies, is more suited for optically thick material. Its definition directly comes from the LTE and diffusion-approximation limits, leading to Eq. (3). The Planck gray opacity, weighting more the opaque frequencies, works better in an optically thin environment. Both approximations are still not realistic, as gray opacities always assume thermal transport of radiation. In BNSM ejecta, the detailed structure of atomic energy levels and the corresponding absorption and emission lines play a crucial role in how radiation interacts with the material. Our most elaborate opacity model (LANL) explicitly includes the information from the opacity database into Eq. (2), discretized as

$$\begin{array}{c}\hfill L=4\pi r^2{\displaystyle \sum _{j_{\mathrm{gr}}}}{F}_{j_{\mathrm{gr}}}\delta {\nu }_{j_{\mathrm{gr}}}+L(\nu <{\nu }_1)+L(\nu >{\nu }_{n_{\mathrm{gr}}}),\end{array}$$(23)

where j_gr = 1, …, n_gr runs over the frequency groups and L(ν < ν₁) and L(ν > ν_ngr) represent the contributes to the integral from frequencies outsides the extreme values of the groups. The terms can affect the results essentially only close to the upper temperature boundary of the opacity tables, where the blackbody function is not yet negligible at high frequencies. To extend the integral out of the group limits, we assume the opacity to be the same as the last available group, and explicitly calculate the other terms. We apply the same method every time an integration over the frequency space is required.

2.6. Jet energy deposition

If a GRB is launched from the polar regions of the remnant, it releases additional energy into the ejecta it traverses. kNECnn can mimic this effect by adding an extra energy term to the right-hand side of Eq. (1) for the innermost shells in the polar regions of the ejecta. We use the “thermal bomb” model already implemented in SNEC (Morozova et al. 2015). The energy of the bomb is released during a chosen time interval and is exponentially attenuated both in time and with distance from the inner boundary. We define two parametrized prescriptions for the energy released by the GRB.

In the “simple jet” model, we assume a jet with an opening angle θ_j to release a uniform (in θ) isotropic energy E_iso(θ) = E₀ in the traversed regions. The total energy released by the jet is given by E_j = (ΔΩ_j/4π)E₀, where ΔΩ_j = 2π(1 − cosθ_j) is the solid angle covered by the jet (or two times this value if one assumes reflection with respect to the equatorial plane).

In the “structured jet” model, the bomb parameters are fixed assuming that the isotropic energy released from the GRB is proportional to the kinetic energy of the structured jet described by Ghirlanda et al. (2019), given by E_iso(θ) = E₀/[1 + (θ/θ_j)^5.5]. The total energy released by the jet can be obtained by dividing E_iso(θ) by the scaling factor λ_θ introduced in Sect. 2.1, and then integrating over the solid angle. The energy release is taken to be constant over the duration Δt_jet of the GRB emission, which starts at t_jet after the start of the simulation.

2.7. Initial profiles

Initializing kNECnn requires specifying an initial Lagrangian profile. This profile must include the density ρ, temperature T, radial velocity v, entropy s and electron fraction Y_e as a function of the enclosed mass m (including the mass of the remnant), and the position of the innermost fluid element r₀. An estimate for the expansion timescale is also needed if the heating rate fits from Wu et al. (2022) are employed. If the online NN is activated, the nuclear composition for each shell is initialized from the other thermodynamic quantities as described in Sect. 2.2. From these initial data, kNECnn automatically defines a new profile with a prescribed number of mass shells n_sh from a selected mass distribution function (i.e., the relative amount of mass δm_i included in each shell i = 1, …, n_sh, rescaled by the total mass of the ejecta M_ej) by performing a linear interpolation. The initial positions r_i of the fluid elements are calculated by progressively layering the mass shells as prescribed by the constraint m(r) = 4π∫_r0^r_iρ(r) r²dr introduced in Wu et al. (2022).

Our simulations are initialized with initial ejecta profiles extracted from ab initio NR simulations performed with the THC code and including a microphysical EOS, M1 neutrino transport, and magnetic-field induced turbulence (Radice et al. 2018; Perego et al. 2019; Bernuzzi 2020; Radice et al. 2022; Zappa et al. 2023; Bernuzzi et al. 2025). The main parameters of the considered binaries are listed in Table 1. All models have already been presented in Bernuzzi et al. (2025) and Gutierrez et al. (2025). For details on the EOS used, see Hempel & Schaffner-Bielich (2010), Typel et al. (2010), Bombaci & Logoteta (2018), and Logoteta et al. (2021). Both selected models BLh_1.43 and DD2_1.67 describe asymmetric BNSM. The chirp mass of the BLh_1.43 configuration is compatible with GW170817, and the mass ratio lies at the upper bound of the constraints inferred with low spin priors Abbott et al. (2019). The DD2_1.67 represent a comparable but more extreme mass-ratio scenario. The neutron stars composing the binary systems are all nonrotating. The central remnant is a massive neutron star for the duration of the simulation in both models.

Throughout each simulation, the unbound material is identified via the Bernoulli criterion (e.g., Nedora et al. 2021) and collected at an extraction radius r_ext ≃ 443 km. Both models eject matter for the full duration of the simulations until ∼100 ms post merger. In order to prepare the subsequent long-term evolution, the ejecta properties are averaged over the azimuthal coordinate, while the polar dependence is retained by discretizing the polar angle into 51 uniformly spaced angular sections. The entropy is recalculated from the original EOS to avoid averaging artifacts and ensure thermodynamic consistency. The extracted time-dependent ejecta profile is then mapped into a Lagrangian profile by positioning the latest shell at r₀ = r_ext, and progressively layering previously ejected shells on top.

In constructing the Lagrangian ejecta profiles, we neglect non-radial components of the ejecta velocity. Moreover, owing to the ray-by-ray nature of the simulations, possible fallback at later times is largely suppressed, since outer shells can only fallback if the inner shells are also infalling. We do not expect these assumptions to introduce strong qualitative corrections to the results discussed in this paper, and defer a quantitative assessment of their impact on the ejecta evolution and final observables to future work.

The initial profile extracted from the BLh_1.43 simulation is shown in Fig. 2. This model produces a massive, neutron-rich (Y_e, 0 ≲ 0.22) tidal tail along the equatorial regions. Emissions from the long-lived remnant drive a proton-rich (Y_e, 0 > 0.5) neutrino-driven wind in the polar regions above the disk persisting until the end of the simulations. The DD2 case is qualitatively similar, but the higher mass ratio strongly enhances the tidal-tail dynamical ejection, while less matter is ejected in the polar neutrino-driven wind. In what follows, we mainly focus on the BLh profile to systematically asses the impact of refined nuclear burning and atomic opacities models on the ejecta dynamics and nucleosynthesis. The DD2 case is also considered to verify the robustness of our results.

Fig. 2. Thermodynamic profile of the ejecta at the beginning of the kNECnn simulation for the BLh_1.43 binary. From left to right, the quadrants show the initial density, temperature, entropy, radial velocity, and electron fraction.

3. Results

We analyze a series of simulations performed on the BLh and DD2 ejecta profiles to explore the physical impact on the nucleosynthesis and light-curve calculations from the introduction of the in situ NN, as well as the effects of switching between the ST05 and DT08 thermalization schemes or among the opacity models described in Sect. 2.5. The simulations considered in this section and their respective configurations are listed in Table 2.

3.1. Effects of dynamics on nucleosynthesis

The online nuclear calculations implemented in kNECnn allow for a coupled evolution between the ejecta hydrodynamics and its composition. In Magistrelli et al. (2024), we showed that prescribing the density evolution with the analytical model of Lippuner & Roberts (2015) can lead to inaccurate nucleosynthesis predictions. Capturing the detailed early (≲500 ms) dynamics of the outflow is essential for reliably determining the nuclear yields and the abundance evolution across the entire nuclide chart. On the other hand, the self-consistently computed nuclear heating can, in principle, feed back on the ejecta evolution. Access to the instantaneous isotopic composition is also required to include detailed physical models for the thermalization of the released nuclear power and the optical opacity. In this and the next sections, we quantify the influence of the dynamics on nuclear calculations, and assess the impact of the improved nuclear power alone on the ejecta dynamics and on the resulting light curves. We deliberately exclude here the composition-dependent thermalization and opacity models (made possible only by the NN coupling), and defer the discussion of their effects to Sects. 3.3 and 3.4. We therefore only focus in this section on the simulations labeled ThK-S and Apr2 in Table 2, both for the BLh and DD2 ejecta profiles.

For the runs without a coupled NN, we compute the nucleosynthesis using two different post-processing procedures. In the first approach (pp_grid), following the same method explored in Magistrelli et al. (2024), we take the Y_e, s, τ grid defined in Perego et al. (2022) and add one dimension for the temperature. We include 15 linearly spaced points between the minimum and maximum temperatures recorded in the initial kNECnn profile. To account for proton-rich winds, we also extend the Y_e dimension with 11 additional points uniformly distributed between 0.5 ≤ Y_e ≤ 0.7. We then map all of the initial fluid elements onto this grid and run an independent instance of SkyNet for each of the grid points, prescribing an exponential-plus-homologous expansion as in Lippuner & Roberts (2015), and defined from the initial ejecta profile as in Perego et al. (2022). The NSE initialization employs the same backtracking prescription used in the full kNECnn run. The final results are obtained by averaging the isotopic abundances over the grid, weighting each point by the total mass of the represented fluid elements. In the second post-processing method (pp_traj), we instead evolve a NN along the thermodynamic trajectory of each individual fluid element from the original kNECnn simulation, and then again mass-weight the abundances.

Comparing the pp_grid procedure against the online results highlights the effects of the detailed early (≲500 ms) dynamics on the nucleosynthesis. During this phase, the expansion of the outflow is shaped by the details of the initial explosion and is dynamically influenced by the released nuclear energy as well as by interactions among fluid elements. Comparing instead the pp_traj method against the online results allows us to isolate the impact of self-consistent nuclear-energy feedback on the dynamics, and consequently again on the nuclear calculations, relative to a simpler analytical prescription for the nuclear power.

We first investigate how the detailed hydrodynamic evolution impacts the nuclear calculations. In Fig. 3 we compare the final nucleosynthesis yields obtained with in situ NN against those predicted using nuclear-power fits and the two post-processing methods. We also test the pp_grid method by doubling the number of points in the Y_e, s, τ dimensions to assess the impact of grid resolution on the final abundances. The two asymmetric simulations give similar nucleosynthesis, whose general features are confirmed between among the different procedures. Strong r-process reactions dominate the equatorial, low-Y_e regions of the tidal tails. The light (A ≲ 110) element peaks, mainly produced in the ejecta emitted above the remnant disk (θ ≳ 30 degrees), are more pronounced in the BLh profile, which features a less prominent tidal tail.

Fig. 3. Mass-weighted nucleosynthesis yields at t ≃ 5 × 104 s for the BLh (left column) and DD2 (right column) ejecta profiles. Top panels: Final yields for the simulations performed with in situ NN (blue) or using analytical fits for the nuclear power and post-processing the nucleosynthesis either on the original thermodynamic trajectories (orange), or on analytical density evolutions prescribed from a grid of initial thermodynamic conditions, with a standard (green) or double (red) resolution on the Ye, τ, s grid. The black dots represent the Solar r-process abundances (Prantzos et al. 2020). All the abundances are scaled to produce a unitary total abundance of the elements with 170 ≤ A ≤ 200. The histogram shows, for the in situ NN run, the global cumulative mass fractions of the first (blue), second (orange) and third (purple) r-process peaks, and the rare-earths (brown). Bottom panel: relative differences against the simulation with in situ NN. The upward [downward] triangles indicate discrepancies greater than two orders of magnitude [smaller than 5 × 10−4]. The horizontal dashed lines highlight the 1% and factors of 1 and 10 discrepancies.

The pp_grid procedure leads to inaccurate final nuclear yields (even when the effect of the released nuclear power is included in the temperature evolution). Relative differences of ≳1 appear across the entire range of mass numbers, with deviations approaching an order of magnitude for several isotopes within the second and third r-process peaks. A noticeable shift of the third peak is also evident when comparing these results to those from the simulation employing in situ NN. The grid discretization has only a secondary effect, as doubling the resolution on the Y_e, τ, s dimensions does not significantly affect the results.

The pp_traj analysis reproduces most of the final abundances within a ∼20% error. This improvement comes from the more accurate representation of the density (and temperature) evolution along each trajectory. If (opacity and) thermalization effects are neglected, and if reliable analytical fits for the nuclear powers are employed, a sufficiently large sample of tracers, well resolved in time up to ∼1 s, can therefore yield reliable nucleosynthesis predictions (if non-radial motion is neglected). The remaining discrepancies originate from corrections to the thermodynamic evolution of the fluid elements that come exclusively from the improved nuclear-power calculation provided by the in situ NN (see the next section for details). The effect on the final mass fraction of lanthanides is on the order of 10%, as reported in Table 2.

In Fig. 4 we compare the time evolution of selected species obtained from the runs with in situ NN to those from the post-processing tracers procedure. We focus on neutrons and protons as diagnostic isotopes, ⁵⁶Ni (half-life τ ≃ 6.077 days), ⁹²Sr (with τ ≃ 2.611 hr), ¹²⁸Sb (with τ ≃ 9.05 hr) and ¹³¹I (with τ ≃ 8.02 days) for their potential contributions to γ-ray emission (Korobkin et al. 2020); Bernuzzi:2024mfx; Jacobi:2025eak, ⁶⁰Fe, ¹²⁹I, ²⁴⁴Pu and ²⁴⁷Cm for their relevance in galactic chemical evolution studies (Côté et al. 2021; Wallner et al. 2021; Davis 2022; Chiesa et al. 2024), and the cumulative mass fractions of lanthanides and actinides to monitor strong r-process nucleosynthesis. At t  ∼  1 s, the system departs from (n, γ)↔(γ, n) equilibrium. The neutron-to-seed ratio consequently drops, triggering neutron freeze-out. The remaining free neutrons undergo β-decay on the t ∼ 10 minutes timescale. ⁵⁶Ni and lanthanides are already present at t = 0. ⁵⁶Ni is further synthesized on the t ∼ 10 ms timescale, while lanthanides are produced for t ≳ 100 ms together with actinides. The isotopes ¹²⁸Sb, ¹²⁹I, and ¹³¹I are produced at significant abundances on a timescale of t ∼ 1 hour, with ¹²⁸Sb already decaying at t ∼ 10 hours.

Fig. 4. Time evolution of the mass-weighted abundances for selected isotopes and cumulative abundances of lanthanides and actinides. The first and second columns refer to simulations ThK-S_BLh and ThK-S_DD2, respectively. The top and bottom rows show the results obtained with in situ NN and their relative differences with respect to the post-processed abundances computed with the tracer method. The horizontal dashed lines in the bottom panels indicate factors-of-unity deviations and the case ⟨Yi⟩=⟨Yipp⟩. Differences are displayed only at times for which ⟨Yi⟩> 10−12.

The neutron freeze-out occurs earlier in the online calculations than in the post-process predictions for the BLh profile, and later in the DD2 case. In the two cases, the post-processing approach overestimates or underestimates, respectively, the amount of remaining free neutrons by ∼10%. The final neutron abundances agree within a few tens of percent, with the online NN consuming more neutrons during the last stages of the evolution. For both profiles, lanthanides and actinides are overproduced by the post-processing calculations by a few tens of percent. The same level of discrepancy represents an approximate upper limit for the other selected elements, but significantly larger deviations (up to two or more orders of magnitude) can be observed at early times (t ≲ 1 s). The abundance evolution of ¹²⁹I seems unaffected by the online calculations.

3.2. Effects of nucleosynthesis on dynamics and light curves

We check the impact of the in situ calculation of the nuclear power on the ejecta hydrodynamics by analyzing the density, temperature, and thermalized heating rate evolution using the histograms in Fig. 5. At early times, the thermodynamic evolution of the outflow is dominated by the initial explosion and is only weakly affected by nuclear burning. At later times, nuclear heating becomes the main source of energy in the ejecta, increasing the temperature of material with initial electron fraction Y_e, 0 ≲ 0.1 by a factor of two or more due to the decay of freshly synthesized neutron-rich isotopes. For t ≳ 1 day, nuclear heating also raises the temperature in regions with Y_e, 0 ≳ 0.45, when ⁵⁶Ni begins to undergo β-decay.

Fig. 5. Mass-weighted histograms of the density, temperature, and instantaneous thermalized heating rate obtained from the BLh runs using the in situ NN (ThK-S_BLh) or the analytical fits for the nuclear power (Apr2_BLh). The top and bottom groups of plots show the ejecta at early (pre neutron freeze-out) and late (t ≃ 10 days) times. The left [right] column shows the profiles as functions of the polar angle [initial electron fraction]. The gray lines in the heating rate panels separate negative and positive values.

The analytical fits fail to reproduce the detailed structure and spread of the nuclear power distribution (especially at early times) and tend to cluster fluid elements around the regions already highly populated in the in situ NN case. Nevertheless, they globally recover the correct orders of magnitudes of the heating rates. Consequently, the impact on the hydrodynamic evolution is limited, yielding only minor deviations in the density and temperature profiles over time. As a result, the hydrodynamics corrections introduced by the online NN have a negligible effect on the nucleosynthesis predictions.

The differences in the predicted evolution of the photosphere temperature and the nuclear power released in optically thin regions, however, have a qualitative impact on the light curves. In Fig. 6 we compare kilonova predictions computed with and without in situ NN for the BLh and DD2 ejecta profiles. The light curves are brighter, peak later, and last longer in the DD2 case, which produces more than double the amount of ejected mass, as expected from the Arnett’s law (Arnett 1982). Both simulations show a blue peak at t ∼ 1 hour, when the photosphere enters regions previously heated up by the decay of freshly synthesized r-process isotopes (Bernuzzi et al. 2025). The late-time red kilonova is powered both by the ongoing decays of these neutron-rich isotopes, and by the ⁵⁶Ni →⁵⁶Co →⁵⁶Fe decay chain (Jacobi et al. 2026).

Fig. 6. Predicted AB apparent magnitudes in the Gemini u, g, i, and Ks bands for an observer at a polar angle of 30 degrees and a distance of 40 Mpc. Comparison between the simulations using in situ NN (solid lines) and analytical nuclear power fits (dashed lines). All models adopt simple (no composition-dependent) thermalization and opacity models. The left and right panels show results for the BLh and DD2 ejecta profiles, respectively.

Calculating the nuclear power on the fly generally shifts the emission peak to later times and lower frequencies (see also Table 2). The lower luminosity at t ∼ 1 hour is caused by the colder temperature predicted by the online NN for the photosphere, which most contributes to the light curves on this timescale. We show this effect in Fig. 7, which displays the heating rate and temperature profiles for the BLh simulation at t ∼ 1 hour. At this time, nuclear contributions from the optically thin layers of the ejecta provide only a minor correction to the luminosity, and differences between the in situ and ex situ simulations arise only in low-activity regions. Similar considerations hold for the DD2 profile.

Fig. 7. Thermalized heating rate and temperature profiles of the BLh ejecta at t ∼ 1 hour. The top panels show the results obtained with in situ NN, and the bottom panels show the relative differences with respect to the run using analytical nuclear power fits. The dashed black and the dotted gray lines represent the position of the photosphere in the runs with in situ NN and fits, respectively.

The plateau in the blue, green, and yellow filters at t ∼ 1 day is sustained by the decays of ⁵⁶Ni and, more generally, by the reactions occurring in the initially mildly neutron-rich and proton-rich fluid elements. Their effect is overestimated by the simulations relying on the analytical nuclear power fits. At later times, the decays of r-process isotopes also contribute significantly to the EM emission. The online NN slightly enhances the light curves across all frequencies with respect to the offline calculations, but predicts fainter red kilonovae at very late times (t ≳ 10 days). The dimmer luminosity at t ∼ 10 days is consistent with the heating rates presented in the bottom part of Fig. 5, where the offline calculations are shown to slightly underestimate heating rates all across the ejecta. Our comparative analysis is not affected by non-LTE effects, which become relevant at these late times, when the ejecta enters the nebular phase (e.g., Pognan et al. 2022; Ricigliano et al. 2025).

While the light curves appear reasonably accurate at late times in this stage of the analysis, tracking the ejecta composition is needed for the detailed calculations of the improved thermalization factors in the DT08 scheme and for the JK, LANL-R, LANL-p, and LANL opacity models. Their relevant qualitative impacts on the light curve predictions are discussed separately in Sects. 3.3 and 3.4, while the overall improvements enabled by the online NN are discussed in Sect. 3.6.

3.3. Effects of thermalization treatment

The thermalization prescription directly enters the nuclear term in Eq. (1). However, nuclear heating remains a subdominant contribution to the dynamics for most of the ejecta evolution. We show this behavior in Fig. 8, which displays the instantaneous average contributions to the internal energy of the fluid for the BLh profile under different thermalization and opacity prescriptions. On the timescales relevant for the nucleosynthesis, the effects of the different thermalization schemes on the temperature and density evolution are almost negligible. As a result, for the BLh profile, the final nuclear yields agree within roughly ∼10%, and the final mass fraction of lanthanides is consistent to within a few percent (see Table 2). The isotopic evolution is similar for both thermalization schemes. Neutron freeze-out occurs slightly earlier (by ∼200 ms) in the simple thermalization case, with no qualitative change in the rest of the isotopic evolution.

Fig. 8. Contributions to the internal energy of the BLh ejecta for the simple and detailed thermalization schemes described in Sect. 2.4 (left column), and for selected opacity models from Sect. 2.5 (right column). Top: Average pressure work (blue), nuclear heating (orange), and radiation energy (red) per unit time and per fluid element, mass-weighted over the ejecta. Bottom: Ratio of nuclear heating and radiation energy over pressure work. The horizontal dashed line marks a 30% ratio.

At t ≳ 1 day, the simple ST05 thermalization scheme significantly overestimates the fraction of nuclear power that is actually thermalized and thus contributes to the temperature and EM emission of the ejecta. This result is independent of the specific choice of the constant $f_{\mathrm{th}}^{\mathtt{kNEC}}$, as the thermalization efficiency of all decay products drops exponentially on the days timescales. At these times, the outflow begins to become optically thin, and the thermalized energy is promptly emitted as thermal radiation from the layers above the photosphere. As shown in Fig. 8, the radiative power is predominantly supplied by nuclear heating at these timescales, while the photospheric contribution becomes negligible. The effect on the light curves is evident from Fig. 9, which displays the kilonova predicted by the two thermalization schemes. The simple ST05 prescription boosts the luminosity in all filters at late times compared to DT08, giving a very long-lived kilonova. The effect is particularly pronounced in the low-frequency bands, as the ejecta has cooled enough at these times (T ∼ 2000 K) to have its blackbody spectrum peaking in the infrared.

Fig. 9. Predicted AB apparent magnitudes in the Gemini u, g, i, and Ks bands for an observer at a polar angle of 30 degrees and a distance of 40 Mpc. Comparison between the BLh simulations using in situ NN and either the time-dependent, composition-based thermalization scheme (solid lines) or a constant thermalization factor (dashed lines).

At early times, the thermalization efficiency is close to unity for all decay products. Any constant, fractional thermalization factor therefore underestimates the nuclear contribution to the internal energy of the ejecta. The resulting lower temperature of the photosphere explains the dimmer and redder light curves produced by the simple thermalization scheme at t ∼ hours and the lower bolometric luminosity reported in Table 2. The two models cross between roughly ∼ 6 hours and ∼ 1 day, when the effective averaged thermalization factor from the online calculations drops to ⟨f_th⟩∼0.5.

Figure 10 shows the mass-weighted contributions to the nuclear power and thermalized heating rate from the species considered in Sect. 2.4, and their thermalization efficiency over time as calculated by the DT08 model. At the t ≲ 1 s timescale, the nuclear power is dominated by the residual term. In the neutron-rich parts of the outflow, contributions to this term come from neutron emissions from the (γ, n)↔(n, γ) equilibrium and β decays of very neutron-rich isotopes not included in the ENDF/B-VIII.0 database. In the proton-rich regions, proton emissions contribute to the residual nuclear power on the other side of the valley of stability. The same considerations hold for the heating rate, as all injected particles thermalize efficiently at very early times. At 1 s ≲ t ≲ 1 hour, the nuclear power is approximately equally injected in electrons, positrons, prompt γ-rays and partially neutrinos. The latter are removed from the heating rates, while the hierarchy of the other contributions is influenced by the heating rates from t ≳ 1 hour, when the specific thermalization factors begin to show qualitatively different behaviors. At t ≳ 1 hour, α decays also become relevant. Their contribution is partially overestimated, as some α emitters are not correctly included in the REACLIB database used by SkyNet, but they are accounted in the ENDF/B-VIII.0 database. The effect is a boost of the relative fraction of energy associated with α particles (the total nuclear power is not affected, as it is directly obtained as an output from the NN). The effects on the total heating rates are expected to be minor, as while an overestimated fraction of the nuclear power will be thermalized as α particles, their thermalization efficiency drops exponentially in time.

Fig. 10. Mass-weighted contributions to the nuclear power (top panel) and heating rates (central panel) from α particles, electrons and positrons, γ-rays (direct emission and e−/e+ annihilation), neutrinos, and all remaining sources of injected energy in addition to their thermalization efficiency (bottom panel) over time for the GK simulation.

At t ≳ 1 day, the residual term begins again to be comparable with the other sources of thermalized heating rates. Only part of this behavior is due to physical contributions from the recoil energy of daughter nuclei and fission fragments, which are expected to always thermalized efficiently. Part of the late-time residual nuclear power is due to small discrepancies between the nuclear inputs in SkyNet and the ENDF/B-VIII.0 database. This component becomes relevant in the heating rate because of its constant thermalization efficiency ($f_{\mathrm{th}}^{\mathrm{oth}}=0.8$ in our simulation), whereas the other energy contributions are either strongly suppressed by thermalization factors that are ≲5 × 10⁻² for charged particles and ≲5 × 10⁻³ for γ-rays, or lost entirely (e.g., neutrinos). If a small fraction of these suppressed contributions is incorrectly reassigned to the residual heating channel, it will thermalize too efficiently and consequently lead to an overestimate of the late-time luminosity. The same limitations exist in the simpler thermalization method of Wu et al. (2022), which assumes a constant efficiency $f_{\mathrm{th}}^{\mathtt{kNEC}}=0.5$ for all channels. Our improved DT08 model addresses this issue to a large extent, but a residual uncertainty still persists at late times. For t ≲ 3 days, the light curves are overestimated by no more than a factor of 40%. This upper limit increases at later times. The systematic uncertainty can in principle be reduced by introducing a time-dependent thermalization factor for the residual term, which must take into account the fission contributions and must decrease monotonically to suppress the deviations between SkyNet and ENDF/B-VIII.0. Nonetheless, by these epochs the LTE assumption underlying our radiative-transport model becomes unreliable, and the thermalization details become a minor source of uncertainty. This uncertainty does not affect the comparison carried out in the other sections of this paper, as the systematics are the same in all simulations.

3.4. Effects of opacity models

The (optical) opacity enters the luminosity term in Eq. (1), and can influence the ejecta dynamics. Consequently, a different opacity model could also indirectly affect the composition evolution. However, since radiation contributes only a subdominant fraction of the energy balance (see Fig. 8), a detailed multigroup treatment of the opacity is not necessary to accurately capture the hydrodynamics. The radiation term becomes relevant only at t ≳ 1 days, when the expansion is already fully homologous and changes in the composition of the outflow are driven only by nuclear decays. An order-of-magnitude estimate of the gray opacity is therefore sufficient to correctly reproduce the ejecta dynamics and perform reliable nucleosynthesis calculations. The relative difference between the abundances predicted by different opacity models always remain below a few 0.1%, and the final mass fraction of lanthanides is essentially unchanged (see Table 2).

Figure 11 shows the evolution of the effective gray opacity for a representative set of fluid elements spanning the full BLh ejecta for the opacity models introduced in Sect. 2.5. The W[22] model (used for the SK simulation) assigns to each fluid element a constant opacity determined only by its initial electron fraction. At early times, corresponding to temperatures T ≳ 2 × 10³ K, this prescription overestimates the effective gray opacity in fluid elements that are lanthanide-poor and underestimates it in lanthanide-rich regions when compared to the multi–group results. At late times, W[22] always overestimates the optical opacity.

Fig. 11. Evolution of the gray opacity κgray for a sampling of fluid elements across all angles and depths, from simulations using different opacity models. Top row: Runs using the analytical expression from Wu et al. (2022, left) or the prescription from Just et al. (2022, right). Central row: Simulations with Rosseland (left) or Planck (right) gray opacities computed from the LANL data for neodymium as representative of the lanthanides. Bottom row: Rosseland mean opacity from the multigroup models using only neodymium to represent lanthanides (left) or also uranium as representative of actinides (right). The color bar is cut at κgray = 10−2 cm2g−1 for easier visualization. For each simulation, the top and bottom subpanels show the trajectories of the fluid elements in the density-temperature and density-cumulative lanthanides mass fraction (Xlan) planes, respectively. Markers indicate different orders of magnitudes of Xlan, with circles for Xlan < 5 × 10−5, downward triangles for 10−5 < Xlan < 10−3, upward triangles for 10−3 < Xlan < 5 × 10−2, and stars for Xlan > 5 × 10−2. For the simulations using only neodymium data, the background shows the single-element mean opacity κel calculated from the LANL tables (Rosseland for GK and Nd300, Planck for GKPl). Similarly, the uranium Rosseland mean opacity is shown in the NdU300 panel.

The same reasoning applies when comparing the multigroup calculations with the J[22] model, employed for the JK run. Although this prescription adapts better to the multigroup results based on neodymium (and uranium) data (simulations Nd300 and NdU300), it still overestimates the opacity in the (equatorial) lanthanide-active (X_lan ≳ 10⁻⁵) regions at temperatures around T ∼ 10³ K.

By construction, the Rosseland (LANL-R) and Planck (LANL-P) mean opacities (simulations GK and GKPl) tend to respectively underestimate and overestimate the effective opacity of the medium when compared against the frequency-dependent (LANL) model. As expected, the effective gray opacities of the lanthanide-producing fluid elements from the GK and GKPl simulations follow the structure of the underlying neodymium data. The trajectories with high initial electron fraction do not synthesize lanthanides, and their opacities are therefore set by the time-dependent opacity floor, consistent with the prescription of Just et al. (2022). The effective opacity of the frequency-dependent Nd300 and NdU300 simulations have a qualitatively different behavior with respect to the gray GK and GKPl calculations, showing a strongly enhanced peak around T ∼ 5 × 10³ K, corresponding to 0.1 ≲ t [days]≲1 depending on the fluid element.

The inclusion of uranium data (simulation NdU300) reduces the effective gray opacity of the r-process fluid elements, as the Rosseland opacity of the actinide outside the resonant region lies below our floor value of κ_oth = 0.2 κ_T cm²g⁻¹ (used for the actinides in Nd300, where uranium is not included). In the NdU300 run, the opacity peak at T ∼ 5 × 10³ K is more localized, with a more gradual rise and an earlier decline, with respect to the Nd300 simulation.

The opacity model shapes the evolution of the photosphere. Figure 12 shows the effective gray photosphere for selected angular sections as a function of time. The W[22] model always places the photosphere farther out in the ejecta than the LANL models (simulations Nd300 and NdU300), resulting in a systematically colder photosphere. At early times (T ≳ 2 × 10³ K), this occurs because W[22] overestimates the opacity in the lanthanide-poor regions, where the photosphere is initially located. At later times, as the photosphere recedes into the r-process–active, lanthanide-rich (equatorial) layers, this model continues to overestimate the opacity across all initial electron fractions.

Fig. 12. Position (top), temperature (middle), and spherically symmetry-equivalent luminosity (bottom) of the effective gray photosphere for the BLh profile. Line styles and colors distinguish between different opacity models and angular sections, respectively.

The J[22] prescription reproduces the early opacity behavior more accurately, though it still fails to capture the opacity drop around T ∼ 2 × 10³ K. Consequently, the differences in the photosphere evolution between this model and NdU300 are more pronounced in the equatorial regions (which contribute the most to the EM emission) than in the polar ones.

The nonanalytical models exhibit non-monotonic opacity evolution in the lanthanide-rich elements. This behavior arises either from the opacity peak around T ∼ 5 × 10³ K in the LANL-R and LANL cases, or from the activation of the atomic data at T ≃ 6 × 10⁴ K in the LANL-P model (see Fig. 11). As a result, the photosphere is observed to bounce toward the outer regions of the ejecta at t ∼ 1 hour in the GKPl simulation, and later at t ∼ 5 hours in the GK, Nd300, and NdU300 runs. The photosphere location predicted by the multigroup calculations always lies between those obtained with LANL-R and LANL-P, which are known to respectively under- and overestimate the effective gray opacity.

In Fig. 13 we examine the impact of the opacity model on the kilonova light curves. The colder photosphere predicted by W[22] with respect to the LANL model results in a suppressed emission at t ≲ 1 hour (in agreement with the bolometric luminosities reported in Table 2), with a stronger suppression at higher frequencies that makes the early-time signal appear redder. Nevertheless, the photospheric luminosity and the optical-band magnitudes remain broadly comparable across models, as the larger photospheric radius partially compensates for the lower emissivity of the colder photosphere. At later times, the artificially high constant gray opacity causes the ejecta to become completely optically thin only at comparatively late stages. Therefore, the photospheric contribution to the luminosity persists for longer, thus delaying the decay of the infrared emission at t ≳ 1 day.

Fig. 13. Predicted AB apparent magnitudes in the Gemini u, g, i, and Ks bands for an observer at a polar angle of 30 degrees and a distance of 40 Mpc. Comparison between a series of opacity models applied to the BLh profile. Left: Frequency- and composition-dependent opacities based on atomic calculations for neodymium and uranium (solid lines), versus the Just et al. (2022) prescription (dashed) and the simple Wu et al. (2022) opacity (dotted). Center: Opacities based on the neodymium atomic calculations, multigroup (solid lines) versus gray models (Rosseland or Planck, dashed and dotted). Right: Frequency-dependent LANL opacities including uranium (solid lines) or only neodymium (dashed) compared against the analytical Just et al. (2022) model (dotted).

As discussed above, J[22] performs better than W[22], but still overestimates the opacity of the lanthanide-rich ejecta. As a result, the light curves from J[22] more closely resemble those produced with the LANL models than those from W[22], but they still exhibit a lower peak bolometric luminosity and a delayed shut-off of the photospheric emission and thus a more prolonged signal.

The central panel of Fig. 13 compares simulations using gray (GK and GKPl) and multigroup (Nd300) opacities based on neodymium atomic calculations. As expected, the Rosseland mean opacity performs better in the light-curve calculations at early times, when the ejecta is still optically thick. At later times, the Planck opacity provides a better approximation for the kilonova emission, especially in the higher-frequency bands. The two gray approximations bracket the actual effective opacity, and therefore generally bracket the multigroup light curves across all filters over most of the evolution. In particular, the Planck mean misses the early blue/UV peak at t ∼ 0.1 days, whereas the Rosseland opacity predicts too bright and delayed magnitude peaks, especially in the higher-frequency filters. The peak bolometric luminosity is consequently enhanced and delayed by the LANL-R model (see Table 2). At late times (t ≳ 10 days), the Planck opacity performs worse again because of the delayed recession of the photosphere in a more opaque material (see Fig. 12). This effect dominates over the direct suppression of the emergent flux (e.g., Eqs. (3) and (19)), leading to an inversion of the hierarchy of the models in the red filter at t ≳ 5 days.

The right panel of Fig. 13 shows that J[22] reproduces multigroup results based on neodymium data more accurately than our gray, atomic-physics-based models. This result highlights the necessity of employing frequency-dependent transport when using detailed atomic opacities, as compressing the information from line-resolved calculations into effective gray averages can be misleading to the point of being outperformed by analytical fitting models.

The original work of Kasen et al. (2017) employed atomic calculations only up to the lanthanides, which at least partially explains the disagreement with our most sophisticated opacity model that also includes uranium (NdU300). Fontes et al. (2020) demonstrated differences between semi-relativistic and fully relativistic neodymium opacities. Additional improvements in opacity modeling may be achieved by replacing theoretically calculated energy levels with NIST-calibrated values. Incorporating these calibrated energies may affects the corresponding light curves and spectra at late times.

When only neodymium is incorporated, the effective opacities predicted by the frequency-dependent simulations lie between the Rosseland and Planck mean values, but are qualitatively closer to the latter (see Fig. 11). Correspondingly, the J[22] approximation produces light curves more similar to LANL-P, particularly at late times, when the ejecta becomes optically thin. However, the analytical prescription fails to capture the drop in effective opacity for the r-process fluid elements at T ∼ 2 − 3 × 10⁴ K. As a result, it misses the luminosity peak at t ∼ 1 − 3 hours, caused by regions of the outflow that are hotter than their surroundings due to early (t ∼ 1 s) decays of freshly synthesized r-process elements (Bernuzzi et al. 2025). This missing peak, already visible in the original publication by Just et al. (2022), is more relevant for trajectories with higher lanthanide fractions. Including uranium (NdU300) gives a brighter, earlier and, bluer kilonova peak (see also Table 2), a smoother decay in the blue/UV filters (u and g), a prolonged emission in the red (i) filter, and a faster drop in the infrared (K_s), driven again by the more rapid recession of the photosphere.

3.5. Impact of a polar jet

As a representative example of the effect that a jet as described in Sect. 2.6 can have on our results, we run again the GK model including a polar, structured jet. We fix θ_j = 15 degrees for the opening angle, within the range also explored in Hamidani et al. (2020) for the GRB associated with GW170817. We start the energy injection at t_jet = 200 ms and choose for the jet the same duration as the observed GRB from Abbott et al. (2017a), Δt_j = 100 ms. We impose for the isotropic energy parameter E₀ = 10⁵¹ erg, compatible with the isotropic luminosity range given by Hamidani et al. (2020).

In Fig. 14 we show the evolution of the radial velocity and temperature, each rescaled to the interval [0, 1] over the plotted time range, for a few representative fluid elements on the timescales relevant for the GRB. The activation of the thermal bomb generates a shock visible in the radial velocity of the innermost mass shells. In the most polar regions, the same fluid elements also experience a rapid increase in temperature. The impact of the jet becomes progressively weaker at larger polar angles, where the deposited energy is smaller. For outer mass shells, the response is delayed, reflecting the finite propagation time of the shock through the ejecta.

Fig. 14. Radial velocity (top) and temperature (bottom) for a set of fluid elements from the GK simulation including the polar jet. Both quantities are rescaled to the interval [0, 1] for visualization purposes. Colors and styles correspond to different angular sections or mass-shell indices, respectively. The vertical black lines indicate the duration of the jet.

As already commented in Magistrelli et al. (2024), the impact of the jet on the globally averaged nucleosynthesis is negligible. Figure 15 shows the nuclear yields from the polar regions (θ ≲ 15 degrees) of the BLh ejecta in the structured-jet and no-jet configurations. Because this angular cut excludes the neutron-rich equatorial outflow, the resulting nucleosynthesis is dominated by the production of elements up to the first peak. The nuclear yields are essentially unchanged by the introduction of the jet, with typical discrepancies on the order of a few tens of percent or less.

Fig. 15. Same as Fig. 3, but now comparing the GK simulation against the same run but including a polar jet (jet), and considering only the material ejected within θ ≲ 15 degrees. The abundances are normalized such that the cumulative abundance of the second r-process peak (111 ≤ A ≤ 145) is unity. The histogram shows the cumulative mass fractions for the simulation with the polar jet.

Figure 16 displays time evolution of selected isotopes, lanthanides and actinides. The overall behavior is characteristic of an only mildly neutron-rich environment, with no significant production of lanthanides, negligible actinides, and only minor synthesis of elements heavier than A ≳ 90. The evolution is similar irrespectively of the jet. Minor differences include a boost in the lanthanides production during the jet activity, and an earlier appearance of ¹²⁸Sb and ¹³¹I on longer timescales. The production of the small amount of actinides is faster in the jet case, while their abundance drops at a comparable rate in the both simulations.

Fig. 16. Time evolution of the mass-weighted abundances for selected isotopes and cumulative lanthanides’ and actinides’ abundances in the GK simulation, considering only matter ejected at latitudes’ θ ≲ 15 degrees. The vertical black lines mark the time interval during which the thermal bomb is active. Top panel: Absolute abundances from the simulation including a polar jet (jet). Bottom panel: Relative differences with respect to the same simulation, but without polar jet. The horizontal dashed lines indicate factors-of-unity deviations and the case ⟨Yi⟩=⟨Yipp⟩. Differences are shown only at times when ⟨Yi⟩> 10−12.

The energy deposited by the thermal bomb in our simulation is not enough to significantly affect the kilonova light curves. Based on the small effects observed in our test, we expect that a more efficient energy deposition or a more powerful jet would enhance the early (t∼ hours), high-frequency optical emission, while producing a dimmer kilonova at later times (t ≳ 1 day), consistent with the findings of Nativi et al. (2020), Magistrelli et al. (2024).

3.6. Best model comparison

In this section we summarize the combined effects of all the improvements that we implemented on top of the code from Wu et al. (2022). We compare our most advanced setup (NdU300), featuring the in situ NN, detailed thermalization, multigroup radiative transfer, and atomic-physics based opacities, against the results obtained using the 2D ray-by-ray extension of original version of the code (Apr2_BLh). For the latter, we also implement the generalized fits for the nuclear power discussed in Sect. 2.3.

As already mentioned in the previous sections, the improved thermalization and opacity prescriptions do not significantly affect the composition evolution of the ejecta. Slightly larger discrepancies appear in the nucleosynthesis yields and abundances evolution with respect to those described in Sect. 3.1, but the physical considerations remain the same. The introduction of the in situ NN is the dominant source of variation in the nucleosynthesis results, and a detailed description of the hydrodynamic evolution of the ejecta is essential for reliable nucleosynthesis predictions. In particular, assuming homologous expansion is not safe on r-process timescales, leading to deviations greater than unity for most mass numbers and a shifted third r-process peaks (see the green and red curves in the bottom panels of Fig. 3). On the other hand, calculating the nuclear power on the fly and including its back-reaction on the hydrodynamics has a comparatively minor impact, typically at the ≲20% level (see the orange lines in the same panels). For the purpose of studying element formation, the hydrodynamics of the ejecta can therefore be reproduced well enough by performing simulations with accurately calibrated nuclear power fits. A post-processing NN can then recover the results of the online calculations within a few 10%, provided that the set of tracers is sufficiently large to properly resolve both the dynamics and the initial composition distribution of the ejecta.

The LANL opacity model has a negligible impact on the hydrodynamics of the ejecta, as already discussed in Sect. 3.4. However, once the composition-dependent thermalization made available by the in situ NN is activated, its combination with the online nuclear calculations leads to qualitative differences in the thermodynamic evolution of the fluid elements. Figure 17 shows the ratios between the two simulations for the histograms of density, temperature and thermalized heating rate at t ∼ 3 days, plotted as functions of polar angle and initial electron fraction. The deviations in the density profile develop during the transition to the homologous expanding phase, at times 0.1 ≲ t [s]≲1. These differences are slightly amplified, but remain qualitatively comparable to those caused by enabling the in situ NN alone (see Fig. 5), since the average thermalization efficiency is of the same order of the simple $f_{\mathrm{th}}^{\mathtt{kNEC}}$ at early times. Once homologous expansion is established, the deviations in the density profiles freeze and remain fixed until the end of the simulation.

Fig. 17. Late-time (t ≃ 3 days) ratios between the mass-weighted histograms of the density, temperature, and instantaneous thermalized heating rate obtained from the BLh runs using the in situ NN and our most advanced thermalization and opacity prescriptions (NdU300), or from the 2D ray-by-ray extension of the configuration of Wu et al. (2022) (Apr2_BLh). The left (right) column shows the ratios as functions of the polar angle (initial electron fraction). Bins containing less than 10−5 M⊙ in both simulations are plotted with increasing transparency for decreasing mass content (considering the run with most mass in that bin). Thus, the color encodes the value of the ratio in ΔM, while the saturation reflects the highest between the two original values ΔM (cf. Fig. 5)

The impact of the in situ NN and the DT08 thermalization scheme is more apparent in the temperature evolution. The two temperature profiles begin to differ after a few hours, and their ratios are initially comparable to those arising solely from the online nuclear calculations. However, at late times (t≳ day), the thermalization efficiency predicted by the advanced model becomes orders of magnitude smaller than the fixed $f_{\mathrm{th}}^{\mathtt{kNEC}}$. Although the nuclear-power fits generally reproduce the correct order of magnitude of the released energy, the simpler thermalization scheme significantly overestimates the heating rate at late times. As a result, it yields temperatures that are a factor of a few higher than those obtained with the composition-dependent thermalization prescription.

The differences in temperature evolution and thermalized heating rates have an immediate impact on the predicted light curves. Figure 18 displays the kilonovae associated with the Apr2_BLh and NdU300 simulations, both showing the characteristic blue/UV peak at t ∼ 1 hour and the later red emission. At early times (t ≲ 0.1 day), the average thermalization efficiency is comparable between the two simulations, and higher in the DT08 scheme. At the same time, the multigroup opacities drive a faster recession of the photosphere. Both effects contribute to a hotter photosphere than in the simpler model. The early kilonova, which is dominated by the photospheric emission, is therefore brighter and peaks at higher frequencies. Moreover, the more efficient thermalization keeps the ejecta hotter for longer, shifting the peak at later times, around t ∼ 2 hours. The thermalization efficiency drops exponentially for t ≳ 0.1 day. At late times (t ≳ 1 day), the advanced model predicts a smaller contribution to the light curves from the energy thermalized in the optically thin layers of the outflow. The material is overall colder, and so is the photosphere (despite its faster recession due to the multigroup opacities). Both the photospheric emission and the contribution to the luminosity from the optically thin layers are therefore suppressed by the composition-dependent thermalization, explaining the steeper late-time decline. The LANL opacity model further reinforces this suppression because of the faster recession of the photosphere, which leads to an earlier drop in the infrared bands.

Fig. 18. Predicted AB apparent magnitudes in the Gemini u, g, i, and Ks bands for an observer at a polar angle of 30 degrees and a distance of 40 Mpc. Comparison between the BLh simulations using in situ NN and our most advanced models for thermalization and opacity (solid lines), or the 2D ray-by-ray extension of the setup of Wu et al. (2022), with their constant prescription for thermalization and opacity, and our updated nuclear power fits (dashed lines).

We stress again that the use of our most advanced models for thermalization and opacity relies on the nuclear information provided by the in situ NN. Our results can be employed to construct effective fits for the heating rate or analytical opacity prescriptions, following the approach of Just et al. (2022), for example. However, because the thermodynamics of the ejecta strongly depends on the properties of the original binary system, a posteriori consistency checks are recommended when possible (e.g., on a selected region of the outflow).

4. Conclusions

In this paper, we carried out a systematic analysis of the impact of fully coupling nuclear calculations to the hydrodynamic evolution of ejecta from BNSMs on their dynamics, nucleosynthesis, and kilonova light curves. We evolved ejecta profiles extracted from ab initio NR simulations up to t ∼ 30 days, and explored different nuclear burning, thermalization, and opacity models. Our most advanced setup coupled the information provided by the in situ NN to composition-dependent, nuclear-informed thermalization calculations and atomic-physics based, frequency-dependent opacities.

By fully coupling the hydrodynamics and nuclear evolution, we have shown that prescribing a simplified exponential-plus-homologous expansion for the ejecta at early times (t ≲ 500 ms) leads to inaccurate predictions of the composition evolution and final yields. However, the back-reaction of the nuclear power calculated on the fly only gives minor corrections to the density evolution, as the radial explosion is dominated by the kinetic energy of the outflow. Simulations in higher dimensions would be needed to asses whether more accurate nuclear-power calculations could be relevant for the motion along the polar and azimuthal directions.

The in situ NN qualitatively impacts the early-time (t ∼ 1 hour) kilonova light curves by altering the location and temperature of the photosphere. As a result, the emission is dimmer, redder, and peaks later compared to models using analytical fits for the nuclear power. These simpler prescriptions perform relatively better in approximating the light curves at later times. However, employing an in situ NN also enables the use of composition-dependent thermalization and opacity models, whose significant qualitative effects on the EM emission are summarized in the remainder of this section. The global corrections introduced by the online nuclear calculations to the density evolution of the outflow have negligible dynamical effects.

Composition-informed thermalization efficiencies are essential for accurately modeling the kilonova brightness and color evolution. Using a constant thermalization factor typically predicts an insufficiently dim and red kilonova at t∼ hours, and strongly overestimates the late-time EM emission from the optically thin regions of the outflow. The thermalization treatment does not significantly affect the nucleosynthesis predictions. A constant thermalization factor only gives ∼10% deviations on the final abundance distribution.

The opacity model also qualitatively affects the kilonova emission. Analytical prescriptions fail to capture the non-monotonic evolution of the effective opacities and typically overestimate them in the regions that determine the position of the photosphere. As a result, they predict a more extended (and thus colder) photosphere. The kilonova is therefore slightly dimmer and shifted toward redder frequencies at early times (t ≲ 1 hour), while the delayed recession of the photosphere prolongs the late-time red emission at t ≳ 5 days. The work from Just et al. (2022) shows that composition-dependent opacity fits can qualitatively reproduce the bolometric luminosity and the late-time red kilonova, but their analytical prescription overestimates the opacity at early times and is thus prone to miss early (t∼ hours) blue peaks. When actinides (represented here by uranium) are included in addition to lanthanides (represented by neodymium), the kilonova exhibits a brighter and bluer early peak (t∼ hours), a smoother decay in the high-frequency bands, and a faster drop in the late-time red emission. The opacity model has a negligible impact on the nucleosynthesis calculations, with discrepancies of ≲10⁻³ between different prescriptions. Since the radiation term is subdominant in the energy equation, the hydrodynamic evolution of the ejecta is already well captured by gray, analytical approximations.

The effect of a polar jet injecting additional energy into the outflow is negligible for all global observables for a realistic set of parameters. More energetic jets or a more efficient release of their energy into the ejecta are expected to produce a brighter and bluer kilonova at early times (t∼ hours), followed by a faster decay on timescales of days.

The overall improvements enabled by the introduction of an in situ NN can be summarized as follows. Resolving the detailed hydrodynamical evolution of the outflow during the first several hundred milliseconds is essential for obtaining reliable nucleosynthesis predictions. Calibrated analytical fits for the nuclear power, together with constant prescriptions for the thermalization efficiency and optical opacity, are sufficient to capture the bulk dynamics of the outflow. However, a more accurate thermalization scheme is required to correctly reproduce its temperature evolution. Moreover, we find that both the composition-dependent thermalization scheme and the multigroup opacity treatment are necessary to accurately model the kilonova emission.

Due to the strong velocity gradients observed in the expanding ejecta, Doppler corrections are expected to modify the observed kilonova. Such corrections also enter in the luminosity term in the energy equation. However, due to the subdominant contribution of the luminosity term to the dynamics (see Fig. 8), we expect only a minor impact on the dynamical evolution of the outflow. In kNECnn, we implemented Doppler corrections in the luminosity term of the energy equation, in the photosphere tracking, and in the light-curve calculations. These corrections are not included in the results presented in this paper. Preliminary tests show a significant increase in execution time mainly in the polar regions, where localized shocks appear and slow down calculation. As expected, our initial investigation indicates an almost negligible impact on the thermodynamic evolution of the outflow and on the nucleosynthesis. The kilonova emission, however, is generally suppressed when Doppler effects are included. The early (t∼ hours) blue peak is damped, particularly at high frequencies, while the infrared emission is only mildly reduced, apart for short-lived luminosity bursts at t ≲ 15 hours. We defer a more detailed analysis of the Doppler effects to a future publication.

In this paper, we used a Lagrangian approach and we did not investigate the role of composition mixing. Some degree of mixing is physically expected, and its relevance is set by the Damköler number D_α, expressing the ratio between the mixing and the (nuclear) reaction timescales (Dimotakis 2005). For D_α ∼ 1, composition mixing and nuclear reactions cannot be decoupled and might affect each other. Lagrangian codes usually assume D_α ≫ 1 (mixing is not relevant), while Eulerian calculations typically overestimate the impact of mixing due to numerical dissipation. In Zhai et al. (2026), the effect of radial mixing is found to negligibly affect the global nucleosynthesis and light curves. In the future, we plan to deepen this analysis using the more advanced opacity models introduced in this paper and to extend it to simulations in higher dimensions.

Acknowledgments

We thank Leonardo Chiesa for making available for this work the updated version of the SkyNet grids used to define the nuclear power fits, and Ryan Wollaeger and Stephan Fritzsche for valuable discussions on atomic opacities. FM acknowledges support from the Deutsche Forschungsgemeinschaft (DFG) under Grant No. 406116891 within the Research Training Group RTG 2522/1. FM also thanks the Trento Institute for Fundamental Physics and Applications (TIFPA) for its hospitality during the development of the project. SB acknowledges funding from the EU Horizon under ERC Consolidator Grant, no. InspiReM-101043372 and from the Deutsche Forschungsgemeinschaft, DFG, project MEMI number BE 6301/2-1. AP is supported by the European Union under NextGenerationEU, PRIN 2022 Project No. 2022KX2Z3B. MJ acknowledges from the Deutsche Forschungsgemeinschaft, DFG, project MEMI number BE 6301/2-1. Simulations were performed on the ARA and DRACO clusters at Friedrich Schiller University Jena, on the supercomputer SuperMUC-NG at the Leibniz- Rechenzentrum (LRZ, www.lrz.de) Munich, and on the national HPE Apollo Hawk at the High Performance Computing Center Stuttgart (HLRS). The ARA cluster is funded in part by DFG grants INST 275/334-1 FUGG and INST 275/363-1 FUGG, and ERC Starting Grant, grant agreement no. BinGraSp-714626. The authors acknowledge the Gauss Centre for Supercomputing e.V. (www.gauss-centre.eu/) for funding this project by providing computing time on the GCS Supercomputer SuperMUC-NG at LRZ (allocations pn36ge and pn36jo). The authors acknowledge HLRS for funding this project by pro- viding access to the supercomputer HPE Apollo Hawk under the grant number INTRHYGUE/44215.

References

Appendix A: Radiative transport equations

The transport equation for the specific intensity I_ν(r,Ω,t) of radiation has the form

$$\begin{array}{c}\hfill \frac{1}{c}\frac{\partial I_{\nu }}{\partial t}+\mathit{n}·\mathrm{\nabla }{I}_{\nu }+{\sigma }_{\nu }{I}_{\nu }=\frac{c}{4\pi }({\sigma }_{\nu }^a{ϵ}_{\nu }^P-{\sigma }_{\nu }^s{ϵ}_{\nu }^r),\end{array}$$(A.1)

where c is the speed of light, n is the direction of propagation, and ${\sigma }_{\nu }^a(\mathit{r},t)={\kappa }_{\nu }\rho$, with ρ density and κ_ν opacity, and ${\sigma }_{\nu }^s(\mathit{r},t)$ are the absorption and scattering coefficients, with σ_ν = σ_ν^a + σ_ν^s. The emissivity is ${ϵ}_{\nu }^r(\mathit{r},t)\equiv c^{-1}{\displaystyle {\int }_{4\pi }}d\mathrm{\Omega }{I}_{\nu }$, with Ω solid angle, and ϵ_ν^P = 4πB_ν/c is the blackbody energy density. Integrating Eq. (A.1) over dΩ gives

$$\begin{array}{c}\hfill \frac{\partial {ϵ}_{\nu }^r}{\partial t}+\mathrm{\nabla }·{\mathit{F}}_{\nu }=c{\sigma }_{\nu }^a({ϵ}_{\nu }^P-{ϵ}_{\nu }^r),\end{array}$$(A.2)

where ${\mathit{F}}_{\nu }(\mathit{r},t)\equiv {\displaystyle {\int }_{4\pi }}d\mathrm{\Omega }\mathit{n}{I}_{\nu }$. If scattering is neglected, σ_ν^s = 0, and Eq. (A.1) reduces to the Kirchhoff law for radiative transport

$$\begin{array}{c}\hfill \frac{1}{c}\frac{\partial I_{\nu }}{\partial t}+\mathit{n}·\mathrm{\nabla }{I}_{\nu }={\kappa }_{\nu }\rho (B_{\nu }-I_{\nu }),\end{array}$$(A.3)

Eq. (A.1) is exact and more general, but integrating Eq. (A.3) over dΩ still leads to the same Eq. (A.2).

Assuming isotropic and quasi-stationary radiation (i.e., radiation locally described by the time evolution of T and ρ), one can manipulate Eq. (A.2) to obtain the isotropic diffusion approximation (IDA) limit:

$$\begin{array}{cc}\hfill {\mathit{F}}_{\nu }& =-\frac{c}{3{\sigma }_{\nu }^a}\mathrm{\nabla }{ϵ}_{\nu }^r,\hfill \\ \hfill \mathrm{\nabla }·{\mathit{F}}_{\nu }& =c{\sigma }_{\nu }^a({ϵ}_{\nu }^P-{ϵ}_{\nu }^r).\hfill \end{array}$$(A.4)

Local thermodynamical equilibrium (LTE) is a sufficient (since it implies IDA) but not necessary condition for this approximation, for which only isotropic radiation on distances smaller than the mean free path (σ_ν^a)⁻¹ of the radiation is required. If LTE, then ϵ_ν^P = ϵ_ν^r and thus

$$\begin{array}{c}\hfill {\mathit{F}}_{\nu }=-\frac{c}{3{\kappa }_{\nu }\rho }\mathrm{\nabla }{ϵ}_{\nu }^P=-\frac{4\pi }{3{\kappa }_{\nu }\rho }\frac{\partial B_{\nu }}{\partial T}\mathrm{\nabla }T.\end{array}$$(A.5)

The treatment in Levermore & Pomraning (1981) is more general than IDA, as it develops the calculation starting from Eq. (A.2) without assuming isotropic diffusion. In particular, they introduce the normalized specific intensity ψ_ν such that I_ν = cϵ_ν^rψ_ν. In the general case, Eq. (A.1) becomes

$$\begin{array}{c}\hfill \frac{1}{c}\frac{\partial ({ϵ}_{\nu }^r{\psi }_{\nu })}{\partial t}+\mathit{n}·\mathrm{\nabla }({ϵ}_{\nu }^r{\psi }_{\nu })+{\sigma }_{\nu }{ϵ}_{\nu }^r{\psi }_{\nu }=\frac{1}{4\pi }({\sigma }_{\nu }^a{ϵ}_{\nu }^P-{\sigma }_{\nu }^s{ϵ}_{\nu }^r).\end{array}$$(A.6)

Its angular integral still gives Eq. (A.2), as ψ_ν is normalized as ∫_4πdΩψ_ν = 1. Using only the assumption of a slowing varying normalized intensity, ${\partial }_t{\psi }_{\nu }+c\mathit{n}·\mathrm{\nabla }{\psi }_{\nu }=0$, one can find

$$\begin{array}{c}\hfill {\psi }_{\nu }=\frac{1}{4\pi }\frac{1}{R_{\nu }coth{R}_{\nu }-\mathit{n}·{\mathit{r}}_{\nu }},\end{array}$$(A.7)

where

$$\begin{array}{c}\hfill {\mathit{r}}_{\nu }=-\frac{\mathrm{\nabla }{ϵ}_{\nu }^r}{{\sigma }_{\nu }\omega {ϵ}_{\nu }^r},\end{array}$$(A.8)

$R_{\nu }=|{\mathit{r}}_{\nu }|$, and

$$\begin{array}{c}\hfill {\omega }_{\nu }=\frac{{\sigma }_{\nu }^a{ϵ}_{\nu }^P+{\sigma }_{\nu }^s{ϵ}_{\nu }^r}{{\sigma }_{\nu }{ϵ}_{\nu }^r}\end{array}$$(A.9)

is the effective albedo. The flux is given by

$$\begin{array}{c}\hfill {\mathit{F}}_{\nu }=-\frac{{\lambda }_{\nu }c}{3{\sigma }_{\nu }{\omega }_{\nu }}\mathrm{\nabla }{ϵ}_{\nu }^r,\end{array}$$(A.10)

with λ_ν the flux limiting factor rationally approximated by Eq. (4).

In IDA, R_ν → 0, ψ_ν → (4π)⁻¹ and λ_ν → 1. Under the LTE assumption, ϵ_ν^r = ϵ_ν^P, ω = 1, and

$$\begin{array}{c}\hfill {\mathit{r}}_{\nu }=-\frac{\mathrm{\nabla }{B}_{\nu }}{{\kappa }_{\nu }\rho B_{\nu }}.\end{array}$$(A.11)

If the flux limiting factor is kept explicit, one gets

$$\begin{array}{c}\hfill {\mathit{F}}_{\nu }=-\frac{{\lambda }_{\nu }c}{{\sigma }_{\nu }}\mathrm{\nabla }{ϵ}_{\nu }^P=-\frac{4\pi {\lambda }_{\nu }}{3{\kappa }_{\nu }\rho }\frac{\partial B_{\nu }}{\partial T}\mathrm{\nabla }T,\end{array}$$(A.12)

where σ_ν = κ_νρ, with κ_ν total opacity. In spherical symmetry (a single ray-by-ray section is an effective 1D problem), ∇T = dT/dr = 4πr²ρ dT/dm, and thus the flux and flux limiter can be written as in Eqs. (3)-(5).

In the LTE gray radiative transport approximation of Bersten et al. (2011), Morozova et al. (2015), the integrated flux can be written as in Eq. (19) by formally defining the average

$$\begin{array}{c}\hfill R_g=\frac{\int B_{\nu }{R}_{\nu }d\nu }{\int B_{\nu }d\nu }=\frac{\left|\mathrm{\nabla }{T}^4\right|}{{\kappa }_R\rho T^4}=\frac{4\pi r^2}{{\kappa }_R{T}^4}\left|\frac{\partial T^4}{\partial m}\right|,\end{array}$$(A.13)

where we used the formal definition of the Rosseland mean opacity, ${\kappa }_R^{-1}={\displaystyle {\int }_0^{\infty }}{\kappa }_{\nu }^{-1}{\partial }_T{B}_{\nu }d\nu /{\displaystyle {\int }_0^{\infty }}{\partial }_T{B}_{\nu }d\nu$, and the normalization $4\pi {\displaystyle {\int }_0^{\infty }}{B}_{\nu }d\nu =ac{T}^4$. The rightmost expression in Eq. (A.14) only holds in spherical symmetry. Inserting Eq. (A.14) in place of R_ν into Eq. (4) and then plugging the resulting gray λ into Eq. (3) in place of λ_ν gives Eq. (19). For our frequency-dependent radiation-hydrodynamic models, we get Eq. (5) from Eq. (A.12), and

$$\begin{array}{c}\hfill L=-{(4\pi r^2)}^2\left({\displaystyle {\int }_0^{\infty }}\frac{4\pi {\lambda }_{\nu }}{3{\kappa }_{\nu }}\frac{\partial B_{\nu }}{\partial T}d\nu \right)\frac{\partial T}{\partial m}\end{array}$$(A.14)

from Eq. (3) and (2).

Appendix B: Network initialization

To investigate the systematic uncertainties on the nucleosynthesis arising from the NSE-based initialization procedures, we perform the same simulation three times, initializing the nuclear composition either with the backtracking method at T_NSE = 6, 8 GK (BK6 and GK simulations), or directly at the temperature recorded at extraction (CNSE).

In Fig. B.1 we show the composition for the BLh ejecta profile at the beginning of the simulation, during the first phase of r-process nucleosynthesis, and once the production of the heavy nuclei is complete. The cold-NSE initialization misses the early reactions after NSE drop-out (physically happening for many fluid elements during NR simulation times, t < 0). Moreover, the low temperatures at extraction cause the Boltzmann term in the NSE equation to be overestimated (see Eq. (8)). Together, these two effects suppress or inflate the initial abundances of the isotopes around A ∼ 65 and A ∼ 80, respectively, and fail to capture the very early production of elements above the first peak. The BK6 and BK8 initializations are considerably more similar to each other, with differences from the lower temperature qualitatively similar to those discussed for the cold-NSE approach. After initialization and during the r-process nucleosynthesis (1 ms ≲ t ≲ 1 s), the rates of the dominant nuclear reactions are set by the local composition. These rates evolve in a way that counteracts the initial errors, progressively reducing the systematic uncertainty over time. As a result, the system gradually loses memory of the initialization method, with relative discrepancies in the late-time yields falling below unity for most mass numbers. This effect is reflected in the final lanthanides mass fractions reported in Table 2, which agree within a few percent across the GK, BK6, and CNSE simulations. As a notable exception, a group of isotopes slightly lighter than the iron group is overproduced by the more approximated initialization methods. Single isotopes can experience more significant discrepancies. The cold-NSE method tends to leave less free neutrons after freeze-out, due to a stronger initial presence of seed nuclei. In our case, this introduces only minor corrections to the light curves, and only for t ≲ 0.5 days.

Fig. B.1. Top panels: Mass-weighted composition of the BLh ejecta profile after initialization (left column), at t = 0.1 s (pre neutron freeze-out, central column), and t = 5 × 104 s (right column). Blue, orange, and green lines show results respectively from the backtrack (TNSE = 8 GK and TNSE = 6 GK) and the cold-NSE initialization procedures. The black dots represent the Solar r-abundances (Prantzos et al. 2020). The abundances are scaled to match a unitary cumulative abundance for the elements of the first r-process peak, with 72 ≤ A ≤ 94. The vertical dashed lines highlight the position of the elements A = 56. The blue, orange, and purple shaded regions remark the position of the first, second and third r-process peaks, respectively. Bottom panel: relative differences against the BK8 run. The upward [downward] triangles indicate discrepancies greater than two orders of magnitude [smaller than 5 × 10−4]. The horizontal dashed lines mark the factors 0.01, 1, and 10 discrepancies.

The backtracking at the higher temperature T_NSE = 8 GK has in general to be preferred over the alternatives, as it approximately takes into account the otherwise overlooked early out-of-NSE processes. However, backtracking initialization methods assume that the pre-evolution is short enough (t ≪ 1 s) to have an approximately negligible effect on the electron fraction and entropy. We check these hypotheses in Fig. B.2. For the high initialization temperature, around 90% of the ejected mass shows relative discrepancies respectively on the electron fraction and entropy below 10% or 3% after the pre-evolution. For around 75% of the mass, this evolution lasts less than 30 ms. The backtracking assumptions are better satisfied by the low-temperature case. For most of the ejecta, the backtracking is safe and closes before the ignition of the core reactions of the r-process nucleosynthesis. Few percent of the fluid elements register pre-evolution times of 30 ms or more, and a change in the Y_e higher than 10%. Such a discrepancy can be crucial for accessing the strong r-process nucleosynthesis for Y_e ∼ 0.22 and to distinguish between neutron and proton rich material at Y_e ∼ 0.5. Long pre-evolution times also mean missing the NN-hydrodynamics coupling in the first phases of the r-process nucleosynthesis. Figure B.1 provides an estimate of the systematic uncertainties propagated from our initialization methods. In the future, we plan to reduce these uncertainties by initializing the fluid element compositions using tracer particles extracted from the original NR data and preprocessed up to the initial time of the hydrodynamic simulation.

Fig. B.2. Histograms of the discrepancies between the extraction and post-backtracking electron fractions (first column) and entropies (second column), and pre-evolution time tbt (third column) for the BLh ejecta profile. The colors represent different prescriptions for TNSE. The histograms are normalized to the value of the highest bin. The vertical, colored lines mark the cumulative 0.75, 0.9 and 0.99 (mass) fraction of the ejecta. The black dashed line correspond to a 5% discrepancy in the first two columns, or a 30 ms pre-evolution time.

All Tables

All Figures

Fig. 1. Nuclear power as a function of time. All the trajectories have initial entropy s0 ≃ 11 kB baryon−1, expansion timescale τ0 ≃ 8.6 ms, and electron fraction Ye ∈ [0.02, 0.1, 0.2, 0.3, 0.4, 0.46, 0.48, 0.54] indicated by the color map. The vertical dashed lines mark the transition between the early- and late-times fits. Top panel: SkyNet calculations (solid lines) and estimate from the updated (dashed lines) and old (dotted lines) nuclear power fits, respectively from this work and Wu et al. (2022). Bottom panel: Relative differences between the results from SkyNet and from the new (solid lines) and previous (dashed lines) version of the fits. The dashed horizontal lines mark differences of one order of magnitude.

In the text

	Fig. 2. Thermodynamic profile of the ejecta at the beginning of the kNECnn simulation for the BLh_1.43 binary. From left to right, the quadrants show the initial density, temperature, entropy, radial velocity, and electron fraction.
In the text

Fig. 3. Mass-weighted nucleosynthesis yields at t ≃ 5 × 104 s for the BLh (left column) and DD2 (right column) ejecta profiles. Top panels: Final yields for the simulations performed with in situ NN (blue) or using analytical fits for the nuclear power and post-processing the nucleosynthesis either on the original thermodynamic trajectories (orange), or on analytical density evolutions prescribed from a grid of initial thermodynamic conditions, with a standard (green) or double (red) resolution on the Ye, τ, s grid. The black dots represent the Solar r-process abundances (Prantzos et al. 2020). All the abundances are scaled to produce a unitary total abundance of the elements with 170 ≤ A ≤ 200. The histogram shows, for the in situ NN run, the global cumulative mass fractions of the first (blue), second (orange) and third (purple) r-process peaks, and the rare-earths (brown). Bottom panel: relative differences against the simulation with in situ NN. The upward [downward] triangles indicate discrepancies greater than two orders of magnitude [smaller than 5 × 10−4]. The horizontal dashed lines highlight the 1% and factors of 1 and 10 discrepancies.

In the text

Fig. 4. Time evolution of the mass-weighted abundances for selected isotopes and cumulative abundances of lanthanides and actinides. The first and second columns refer to simulations ThK-S_BLh and ThK-S_DD2, respectively. The top and bottom rows show the results obtained with in situ NN and their relative differences with respect to the post-processed abundances computed with the tracer method. The horizontal dashed lines in the bottom panels indicate factors-of-unity deviations and the case ⟨Yi⟩=⟨Yipp⟩. Differences are displayed only at times for which ⟨Yi⟩> 10−12.

In the text

Fig. 5. Mass-weighted histograms of the density, temperature, and instantaneous thermalized heating rate obtained from the BLh runs using the in situ NN (ThK-S_BLh) or the analytical fits for the nuclear power (Apr2_BLh). The top and bottom groups of plots show the ejecta at early (pre neutron freeze-out) and late (t ≃ 10 days) times. The left [right] column shows the profiles as functions of the polar angle [initial electron fraction]. The gray lines in the heating rate panels separate negative and positive values.

In the text

Fig. 6. Predicted AB apparent magnitudes in the Gemini u, g, i, and Ks bands for an observer at a polar angle of 30 degrees and a distance of 40 Mpc. Comparison between the simulations using in situ NN (solid lines) and analytical nuclear power fits (dashed lines). All models adopt simple (no composition-dependent) thermalization and opacity models. The left and right panels show results for the BLh and DD2 ejecta profiles, respectively.

In the text

	Fig. 7. Thermalized heating rate and temperature profiles of the BLh ejecta at t ∼ 1 hour. The top panels show the results obtained with in situ NN, and the bottom panels show the relative differences with respect to the run using analytical nuclear power fits. The dashed black and the dotted gray lines represent the position of the photosphere in the runs with in situ NN and fits, respectively.
In the text

Fig. 8. Contributions to the internal energy of the BLh ejecta for the simple and detailed thermalization schemes described in Sect. 2.4 (left column), and for selected opacity models from Sect. 2.5 (right column). Top: Average pressure work (blue), nuclear heating (orange), and radiation energy (red) per unit time and per fluid element, mass-weighted over the ejecta. Bottom: Ratio of nuclear heating and radiation energy over pressure work. The horizontal dashed line marks a 30% ratio.

In the text

	Fig. 9. Predicted AB apparent magnitudes in the Gemini u, g, i, and Ks bands for an observer at a polar angle of 30 degrees and a distance of 40 Mpc. Comparison between the BLh simulations using in situ NN and either the time-dependent, composition-based thermalization scheme (solid lines) or a constant thermalization factor (dashed lines).
In the text

	Fig. 10. Mass-weighted contributions to the nuclear power (top panel) and heating rates (central panel) from α particles, electrons and positrons, γ-rays (direct emission and e−/e+ annihilation), neutrinos, and all remaining sources of injected energy in addition to their thermalization efficiency (bottom panel) over time for the GK simulation.
In the text

Fig. 11. Evolution of the gray opacity κgray for a sampling of fluid elements across all angles and depths, from simulations using different opacity models. Top row: Runs using the analytical expression from Wu et al. (2022, left) or the prescription from Just et al. (2022, right). Central row: Simulations with Rosseland (left) or Planck (right) gray opacities computed from the LANL data for neodymium as representative of the lanthanides. Bottom row: Rosseland mean opacity from the multigroup models using only neodymium to represent lanthanides (left) or also uranium as representative of actinides (right). The color bar is cut at κgray = 10−2 cm2g−1 for easier visualization. For each simulation, the top and bottom subpanels show the trajectories of the fluid elements in the density-temperature and density-cumulative lanthanides mass fraction (Xlan) planes, respectively. Markers indicate different orders of magnitudes of Xlan, with circles for Xlan < 5 × 10−5, downward triangles for 10−5 < Xlan < 10−3, upward triangles for 10−3 < Xlan < 5 × 10−2, and stars for Xlan > 5 × 10−2. For the simulations using only neodymium data, the background shows the single-element mean opacity κel calculated from the LANL tables (Rosseland for GK and Nd300, Planck for GKPl). Similarly, the uranium Rosseland mean opacity is shown in the NdU300 panel.

In the text

	Fig. 12. Position (top), temperature (middle), and spherically symmetry-equivalent luminosity (bottom) of the effective gray photosphere for the BLh profile. Line styles and colors distinguish between different opacity models and angular sections, respectively.
In the text

Fig. 13. Predicted AB apparent magnitudes in the Gemini u, g, i, and Ks bands for an observer at a polar angle of 30 degrees and a distance of 40 Mpc. Comparison between a series of opacity models applied to the BLh profile. Left: Frequency- and composition-dependent opacities based on atomic calculations for neodymium and uranium (solid lines), versus the Just et al. (2022) prescription (dashed) and the simple Wu et al. (2022) opacity (dotted). Center: Opacities based on the neodymium atomic calculations, multigroup (solid lines) versus gray models (Rosseland or Planck, dashed and dotted). Right: Frequency-dependent LANL opacities including uranium (solid lines) or only neodymium (dashed) compared against the analytical Just et al. (2022) model (dotted).

In the text

	Fig. 14. Radial velocity (top) and temperature (bottom) for a set of fluid elements from the GK simulation including the polar jet. Both quantities are rescaled to the interval [0, 1] for visualization purposes. Colors and styles correspond to different angular sections or mass-shell indices, respectively. The vertical black lines indicate the duration of the jet.
In the text

	Fig. 15. Same as Fig. 3, but now comparing the GK simulation against the same run but including a polar jet (jet), and considering only the material ejected within θ ≲ 15 degrees. The abundances are normalized such that the cumulative abundance of the second r-process peak (111 ≤ A ≤ 145) is unity. The histogram shows the cumulative mass fractions for the simulation with the polar jet.
In the text

Fig. 16. Time evolution of the mass-weighted abundances for selected isotopes and cumulative lanthanides’ and actinides’ abundances in the GK simulation, considering only matter ejected at latitudes’ θ ≲ 15 degrees. The vertical black lines mark the time interval during which the thermal bomb is active. Top panel: Absolute abundances from the simulation including a polar jet (jet). Bottom panel: Relative differences with respect to the same simulation, but without polar jet. The horizontal dashed lines indicate factors-of-unity deviations and the case ⟨Yi⟩=⟨Yipp⟩. Differences are shown only at times when ⟨Yi⟩> 10−12.

In the text

Fig. 17. Late-time (t ≃ 3 days) ratios between the mass-weighted histograms of the density, temperature, and instantaneous thermalized heating rate obtained from the BLh runs using the in situ NN and our most advanced thermalization and opacity prescriptions (NdU300), or from the 2D ray-by-ray extension of the configuration of Wu et al. (2022) (Apr2_BLh). The left (right) column shows the ratios as functions of the polar angle (initial electron fraction). Bins containing less than 10−5 M⊙ in both simulations are plotted with increasing transparency for decreasing mass content (considering the run with most mass in that bin). Thus, the color encodes the value of the ratio in ΔM, while the saturation reflects the highest between the two original values ΔM (cf. Fig. 5)

In the text

Fig. 18. Predicted AB apparent magnitudes in the Gemini u, g, i, and Ks bands for an observer at a polar angle of 30 degrees and a distance of 40 Mpc. Comparison between the BLh simulations using in situ NN and our most advanced models for thermalization and opacity (solid lines), or the 2D ray-by-ray extension of the setup of Wu et al. (2022), with their constant prescription for thermalization and opacity, and our updated nuclear power fits (dashed lines).

In the text

Fig. B.1. Top panels: Mass-weighted composition of the BLh ejecta profile after initialization (left column), at t = 0.1 s (pre neutron freeze-out, central column), and t = 5 × 104 s (right column). Blue, orange, and green lines show results respectively from the backtrack (TNSE = 8 GK and TNSE = 6 GK) and the cold-NSE initialization procedures. The black dots represent the Solar r-abundances (Prantzos et al. 2020). The abundances are scaled to match a unitary cumulative abundance for the elements of the first r-process peak, with 72 ≤ A ≤ 94. The vertical dashed lines highlight the position of the elements A = 56. The blue, orange, and purple shaded regions remark the position of the first, second and third r-process peaks, respectively. Bottom panel: relative differences against the BK8 run. The upward [downward] triangles indicate discrepancies greater than two orders of magnitude [smaller than 5 × 10−4]. The horizontal dashed lines mark the factors 0.01, 1, and 10 discrepancies.

In the text

Fig. B.2. Histograms of the discrepancies between the extraction and post-backtracking electron fractions (first column) and entropies (second column), and pre-evolution time tbt (third column) for the BLh ejecta profile. The colors represent different prescriptions for TNSE. The histograms are normalized to the value of the highest bin. The vertical, colored lines mark the cumulative 0.75, 0.9 and 0.99 (mass) fraction of the ejecta. The black dashed line correspond to a 5% discrepancy in the first two columns, or a 30 ms pre-evolution time.

In the text