Multi-messengers from the radioactive decay of $r$-process nuclei

TL;DR

This paper calculates the emission spectra from the radioactive decay of r-process nuclei, revealing their potential as multi-messenger signals to probe heavy element formation in astrophysical events.

Contribution

It provides the first detailed, first-principles calculations of multi-messenger emission spectra from r-process decay, incorporating nuclear structure and decay data.

Findings

Emission spectra vary significantly over time.

Particles have non-thermal energy distributions.

Signals are complementary to kilonova observations.

Abstract

The radioactive $\beta$-decay of nuclei synthesized in the rapid neutron capture process ($r$-process) releases a variety of particles, including electrons, $\gamma$-rays, neutrinos, and neutrons. These particles provide a rich set of multi-messenger signals that carry information about the astrophysical environments where neutron-rich nucleosynthesis occurs. In this work, we calculate from first principles the emission spectra resulting from the $\beta$-decay of $r$-process nuclei. Our approach incorporates detailed nuclear structure and decay data to model the energy distributions of each particle species. We couple the spectra with a nuclear reaction network simulation to obtain the temporal evolution of these distributions. We find that the emission distributions vary significantly in time and are non-thermal, with substantial average energies. We investigate these nuclear signals as a direct probe of heavy element formation and show that they are complementary observables to kilonova.

Figures

Equations

$${}_Z^{A}X{\to }_{Z+1}^A{X}^{\ast }+e^{-}+{\bar{\nu }}_e$$

$${}_Z^{A}X{\to }_{Z+1}^{A-C}{X}^{\ast }+e^{-}+{\bar{\nu }}_e+Cn$$

$$N^{(i)}(t)=\sum _j{S}_j^{(i)}(E,t)w_j,$$

$${\Phi }^{(i)}(t)=\sum _j{S}_j^{(i)}(E,t)w_j{E}_j,$$

$${\bar{E}}^{(i)}(t)={\Phi }^{(i)}(t)/N^{(i)}(t)=\frac{{\sum }_j{S}_j^{(i)}(E,t)w_j{E}_j}{{\sum }_j{S}_j^{(i)}(E,t)w_j}$$

$$S^{(i)}(E,t)=\sum _{{}_Z^{A}X}{F}_{{}_Z^{A}X}(t)S_{{}_Z^{A}X}^{(i)}(E)$$

$$p(Q_{\beta })\sim p_e^2{p}_{\nu }^{2}F(Z+1,A,E_e)(1-f(E_e))\delta (Q_{\beta }-E_e-E_v)$$

$$\rho (t)={\rho }_0{\left(1+\frac{t}{{\tau }_1}+{\left(\frac{t}{{\tau }_2}\right)}^{\xi }\right)}^{-1},$$

$$T(t)=T_0{\left(\frac{\rho (t)}{{\rho }_0}\right)}^{\gamma -1},$$

$$F^{(i)}(E,t^{\ast })={\Phi }^{(i)}(E,t)(\frac{M_{ejecta}}{M_n})\frac{1}{4\pi D^2},$$

$$\mathcal{N}\approx N_p{\bar{\sigma }}_{\text{IBD}}{F}_{total}^{\nu }(t)$$

$$\frac{dN}{dE}=\frac{3N_{total}{E}_0}{4E^{2}v}(1-\frac{(1-E_0/E{)}^2}{v^2})$$

$$\frac{dm}{dv}=m_{ex}\frac{3v^2}{v_{ex}^2},$$

$$\lambda (v_z)={\lambda }_0(1-v_z),$$

$$dm=\frac{3m_{ex}\dot{ϵ}}{4v_{ex}^3}(v_{ex}^2-v_z^2)d{v}_z$$

$$\frac{dN}{d\lambda }=\frac{dm}{d\lambda }\dot{ϵ}=\frac{3m_{ex}\dot{ϵ}}{4v_{ex}^3}(v_{ex}^2-v_z^2)\frac{d{v}_z}{d\lambda }$$

$$\frac{dN}{d\lambda }=\frac{3m_{ex}\dot{ϵ}}{4v_{ex}{\lambda }_0}\left(1-\frac{(1-\lambda /{\lambda }_0{)}^2}{v_{ex}^2}\right)$$

$$\frac{dN}{dE}=\frac{3N_{total}{E}_0}{4E^2{v}_{ex}}(1-\frac{(1-E_0/E{)}^2}{v_{ex}^2}),$$

Full text

\uatGamma-ray bursts629, \uatNuclear astrophysics1129, \uatNucleosynthesis1131, \uatR-process1324, \uatCompact objects288

show]agross@lanl.gov

I Introduction

It has long been proposed that the rapid neutron capture process (r-process) is responsible for the synthesis of a significant fraction of the elements heavier than iron (see e.g., Burbidge et al. (1957) , Arcones and Thielemann (2023)). The r-process occurs in astrophysical environments with large numbers of free neutrons, and proposed sites include supernovae (see e.g. Woosley et al. 1994; Qian and Woosley 1996; Wanajo et al. 2001; Fryer et al. 2006; Winteler et al. 2012; Nishimura et al. 2015; Siegel et al. 2019; Yong et al. 2021) as well as compact object mergers, e.g. the collision of two neutron stars (NSMs) (see Lattimer and Schramm 1974; Eichler et al. 1989; Rosswog et al. 1999; Freiburghaus et al. 1999; Fryer et al. 2006; Goriely et al. 2011; Korobkin et al. 2012; Hotokezaka et al. 2013; Bauswein et al. 2013; Wanajo et al. 2014). Due to the extreme conditions present in these environments, as well as limited experimental data for many neutron-rich isotopes far from stability where the r-process occurs, it is difficult to obtain precise predictions for r-process production yields; and therefore, the relative contributions of these sites is not well understood.

In recent years, there has been an increased focus on the use of astrophysical observations to provide constraints on the astrophysical location of the $r$-process. The landmark multi-messenger observations of GW170817 (see Abbott et al. (2017)) have given credence to NSMs as being a site of robust r-process nucleosynthesis. This conclusion stems from the observed kilonova, a radioactive transient thought to be powered by the decay of $r$-process nuclei (see Drout et al. 2017; Kasen et al. 2017; Kasliwal et al. 2017; Murguia-Berthier et al. 2017; Tanaka et al. 2017; Waxman et al. 2018; Brethauer et al. 2024). The production of lanthanide elements which have high opacity (e.g., Kasen et al. 2013; Tanaka and Hotokezaka 2013) has been shown to result in a spectrum that peaks in the infrared, and therefore, the observed optical and near-infrared peak was interpreted as evidence that NSMs undergo a strong $r$-process. In contrast, a weaker $r$-process is projected to yield a peak in the blue and optical bands (e.g., Metzger et al. 2010; Roberts et al. 2011; Metzger and Fernández 2014).

Although kilonova are of great interest and importance to the astrophysical community, the necessary multi-physics involved makes the modeling challenging. The numerous poorly constrained parameters in kilonova models can lead to degenerate light curves, obscuring the physical interpretation of these events (see e.g., Barnes et al. 2016; Zhu et al. 2021; Korobkin et al. 2021; Fryer et al. 2023, 2024).

One principal component of modeling the kilonova light curve is the spectral emission profiles of radioactive decay products. The large number of nuclear $\beta$-decays which occur during the course of the $r$-process generates an enormous number of emitted particles: electrons ($e^{-}$), neutrinos $(\bar{\nu}_{e})$, $\gamma$-rays ($\gamma$), and neutrons, in addition to reactions that produce alpha particles, and fission fragments. These particles power the associated radioactive heating, and their energy distribution is therefore vital to understanding the details of their thermalization and contribution to the observed light curve (Barnes et al., 2016; Brethauer et al., 2024; Bulla, 2023). The modeling of these emission profiles is itself a challenge, as it not only requires the use of a sophisticated nuclear reaction network, but also appropriate emission spectra for thousands of nuclear species. As a result, kilonova models generally calculate the energy generation rate for each channel via the Q-value of the decay and then make simple assumptions about how that energy is distributed between particles (see e.g. Korobkin et al. (2021)), although some papers have attempted more detailed models, such as Shingles et al. (2023), which have incorporated detailed $\gamma$-ray emission spectra, albeit only those experimentally measured.

Recent work (Mumpower et al., 2025a) has provided a general purpose framework for calculating the emission spectra for $\beta$-decay, as well as tabulated results for the decay of each nuclear species. In this work, we combine these theoretical spectra calculations with those of a nuclear reaction network to simulate state-of-the-art emission spectra for $r$-process nucleosynthesis. These emission spectra are consistently calculated for all nuclei which undergo $\beta$-decay for all particles: electrons, neutrinos, $\gamma$-rays, and neutrons, and their inclusion in future radiative transport models will enable a more precise interpretation of observations. We organize our paper as follows: In $\S\ref{sec:method}$, we describe the methodology for the calculation of these emission spectra. In $\S\ref{sec:result}$, we present the obtained emission spectra, highlighting the unique features that are observed, and discuss the potential implications for the modeling of kilonova light curves. In $\S\ref{sec:observe}$, we estimate the potential for direct observation of these emission spectra. We summarize and conclude in $\S\ref{sec:conclude}$.

II Methodology

The nuclear reaction for a species undergoing $\beta$-decay can be written as:

[TABLE]

where A and Z are the mass and atomic number of the decaying nucleus. In addition, nuclei often decay into an excited state of the daughter nuclei, leading to the emission of large numbers of $\gamma$-rays during the de-excitation of this daughter nucleus. Excited states of the daughter can spontaneously emit neutrons, leading to $\beta$-delayed neutron emission:

[TABLE]

where $C$ is the number of neutrons that are emitted. Because the energy release for $\beta$-decays of $r$-process nuclei is normally large ($\sim$ MeV), the energies of the emitted particles can also be substantial. As a consequence, these emitted particles can have significant effects on the modeling of the kilonova, and it is therefore important to understand the distribution of these emitted particles as a function of energy and time.

We define the emission spectrum for species $(i)$ as $S^{(i)}(E,t)$ giving the number of emitted particles per unit time per unit energy per unit mass. This convention is adopted because spectra are often binned in energy, therefore the value of $S^{(i)}(E,t)$ is independent of the bin width. The total number of particles emitted per unit time per unit mass $N^{(i)}(t)$ can be calculated as:

[TABLE]

where $S_{j}$ and $w_{j}$ are the emission spectrum and bin width for the $j^{th}$ energy bin. We can calculate the total emission flux by assuming the emitted particles in a given bin are uniformly distributed across the energy range:

[TABLE]

where $\Phi^{(i)}(t)$ is the emission flux per unit mass, which has units of energy per unit time per unit mass, $w_{j}$ and $S_{j}$ are as defined above, and $E_{j}$ is the average energy of the bin. We can also calculate the average energy of emitted particles as a function of time by dividing the total emission flux by the number of particles:

[TABLE]

To calculate this emission spectra from $r$-process nucleosynthesis, we shall combine the results from detailed spectra calculations for the $\beta$-decay of individual nuclei with a reaction network which quantifies the number of $\beta$-decays that occur for each nucleus. We therefore factorize the emission spectra as follows:

[TABLE]

where $S^{(i)}_{{}^{A}_{Z}X}(E)$ is the spectral emission for species (i) of the nucleus ${{}^{A}_{Z}X}$ per unit decay as a function of energy (number per unit energy per decay), and $F_{{}^{A}_{Z}X}(t)$ is the reaction flow for the nucleus ${}^{A}_{Z}X$ as a function of time (decays per second per unit mass), and the sum is over every nucleus which decays in the network. We note the importance of the use of reaction flow to characterize the number of decays per unit mass—while this quantity is in principle just the abundance of the species as a function of time multiplied by its decay rate, in practice, there are many species where the abundance is small but the overall reaction flow and spectral contribution are significant. For a prime example of this, we highlight Tl-208, with a half-life of 3 minutes. It is fed through a long decay chain by Ra-228, with a half-life of 5.75 years, which produces significant spectral contributions even though it is never of significant abundance at such timescales (see Vassh et al. (2024) for more details).

We summarize our methodology for the calculation of each of these two quantities below.

II.1 Individual Emission Spectra

We use the tabulated $\beta$-decay spectra from Mumpower et al. (2025a), which provides the electron, $\gamma$-ray, neutrino, and neutron spectra for neutron-rich nuclei with A $>$ 8. The emission spectra are calculated using a statistical, multi-phase approach. he calculation first obtains a $\beta$-strength function, giving the distribution of excited daughter nucleus states which the decaying parent initially transforms into. If the nucleus has a level structure and $\beta$-decay scheme which are experimentally known, it is used, otherwise the Quasi-particle Random Phase Approximation (e.g., Krumlinde and Möller (1984); Möller et al. (1997, 2003a); Marketin et al. (2016)) is used to estimate the $\beta$-strength function. From this, the emission spectra of the emitted electron and anti-neutrino are calculated under the assumption that the total available decay energy is divided between them (zero kinetic energy for the daughter nucleus). The shape of the spectra is determined from statistical considerations, where the density of states can be written as:

[TABLE]

where $p_{e}$ and $p_{\nu}$ are the momenta of the electron and antineutrino, respectively, $F$ is the Fermi function for $\beta$-decay, which accounts for distortion in the electron wave function due to the interaction between the nucleus and the electron, and $f$ is the Fermi-Dirac distribution, which accounts for Pauli blocking in the ambient electrons. From the shape of the density of states, the total decay energy can be used to normalize the electron spectrum, and via conservation of energy, the neutrino spectrum. This result is then summed and appropriately weighted over the $\beta$-strength function for possible daughter nucleus states to get the total emission spectra.

The $\gamma$-ray and neutron emissions are calculated from a Hauser-Feshbach statistical approach. For each excited daughter nucleus state which can be accessed from $\beta$-decay, the relative branching ratios of $\gamma$-ray emission (to all possible lower energy states) and neutron emission (to all possible states of the new daughter nucleus) are calculated, creating a web of possible decay paths. These decay paths chain into additional daughter nuclei when $\beta$-delayed neutron emission occurs. From these pathways, the system is treated statistically as a weighted sum over all possible paths, and the emission spectra is the corresponding weighted sum of the spectral emissions of the individual path components. For more details of the mathematical treatment of the spectra calculation, we refer the reader to Mumpower et al. (2025a) and the references therein.

The frame work used above assumes that beta decays proceed via “allowed” transitions. First-forbidden transitions are allowed in the calculation of the beta-strength; however, the shape of the spectra is assumed to be the same as allowed transitions. In practice, the assumed shape of the spectra for these transitions would include an additional, energy dependent shape factor which accounts for the orbital angular momentum carried by the leptons. In addition, higher-order transitions are neglected, and if included, would allow different levels to be populated, altering the emission spectra. It is uncertain how much the limitations of our approach will affect the resulting emission spectra. It is suggested that most light-to-mid mass nuclei undergo allowed transitions (see e.g. Möller et al. (2003b)); however, it is unclear how much first-forbidden transitions contribute to heavier elements (see e.g. Marketin et al. (2016); Ney et al. (2020); Ravlić et al. (2025)). We emphasize the importance of further development in the calculation of nuclear structure of heavy and superheavy elements in order to accurately quantify the effects on emission spectra.

II.2 Nuclear Reaction Flow

We simulate nucleosynthesis with version 1.6.0 of the Portable Routines for Integrated nucleoSynthesis Modeling (PRISM) reaction network (Sprouse et al., 2021). PRISM uses a robust adaptive timestepping algorithm which automatically reduces the timestep if the abundance of any species or the conditions of the environment change too quickly, which fully resolves all species in the network numerically, even those with low abundance, which allows for accurate calculation of the reaction flow. The nuclear input to PRISM is based on the 2012 version of the Finite Range Droplet Model (Möller et al., 2012, 2016). Radiative capture and fission rates are calculated with the CoH3 statistical Hauser–Feshbach code (Kawano, 2019, 2021). $\beta$-decay rates, including delayed neutron emission, are calculated assuming statistical de-excitation from excited states (Mumpower et al., 2016, 2018). We highlight that these $\beta$-decay rates are calculated from the same assumptions as in the calculation of the spectra above, such that the probability for a nucleus to undergo delayed neutron emission in the reaction network is the same as in the calculation of the spectra. The remaining reaction rates (e.g. alpha decay and other less substantive reaction types for the $r$-process) are obtained from the REACLIB database (Cyburt et al., 2010). Nuclear fission is handled as in Vassh et al. (2019).

To simulate the composition of r-process material we use the trajectory (b) of Mumpower et al. (2025b), which produces a robust $r$-process with substantial actinide production. Fig. 1 shows the final abundances of this trajectory, which are a good match to the solar $r$-process residuals above the second peak of the $r$-process (A $\sim$ 130).

This trajectory is obtained through modeling the cocoon of a $\gamma$-ray burst, which has been suggested to be a site of the $r$-process due to photo-hadronic interactions in the jet head (Mumpower et al., 2025b). Our choice of this trajectory as opposed to many others from more conventional $r$-process sites which can also produce a robust $r$-process (e.g., Sprouse et al. (2024)) is due to the accuracy with which the solar residuals are reproduced by this trajectory; however, we do not expect that the results which we shall discuss below will change significantly with a more conventional choice of trajectory. The cocoon is modeled with the density profile of:

[TABLE]

where $\xi$ = 2, $\tau_{1}$, $\tau_{2}$ are characteristic timescales, for which we use $\tau_{1}=\tau_{2}=3.5\times 10^{-2}$ s, and $\rho_{0}$ is the initial density, for which we use $3.2\times 10^{4}$ g/cm3. The temperature is assumed to evolve as an adiabatic gas:

[TABLE]

where $\gamma=4/3$ (radiation dominated), and we take $T_{0}=2$ GK. We take the initial electron fraction to be $Y_{e}$= 0.034. For more details of the physical motivation and assumptions behind this trajectory, we refer the reader to (Mumpower et al., 2025b).

III Results

In Fig. 2, we show the particle emission spectra for electrons, $\gamma$-rays, neutrinos, and neutrons as a function of energy at several different times. We note the following conclusions.

The electron and neutrino spectra are almost identical. In the first second of emission, their emission profiles are significant to very high energies ($\sim$ 10-20 MeV). For the first 10- /100 seconds, this trend continues, with significant production up to $\sim$ 5 MeV with an average energy of $\sim 1$ MeV, and then a gradual but persistent falloff to $\sim$ 20 MeV. We note that interesting features in the distribution are consistently present, but not constant as a function of time, such as the bump around $\sim$10 MeV at 100 seconds. On longer time scales, the spectra remains broadly distributed in energy, though the distribution generally continues to decrease in the number of higher-than-average energy particles. We attribute the smooth high energy tail in the spectrum which is seen up to long time scales to the decays of the byproducts of fission, which produce a consistent shower of a wide distribution of particles to long time scales. In general, the neutrino spectrum is slightly higher average energy ($\sim 15\%$) than that of the electron. The consistency in this trend that is observed as well as the similarity of the shape of the neutrino and electron spectra is likely a consequence of the statistical approach taken to partition the available energy between the neutrino and electron.

The $\gamma$-ray spectrum is characteristically different from the electron and neutrino spectra. Initially, the distribution is bimodal in nature, with a strong peak at low energy ($\sim$ 1 MeV), then a sharp fall-off, followed by a flattening to $\sim$ 15-20 MeV. This rapidly morphs into a unimodal distribution by $\sim$ 1 second. Initially, the distribution is generally more smooth, with not many specific spectral lines to be identified. This is attributed to the lack of available nuclear data for the many short-lived nuclei far from stability which decay at early times. This lack of data means that there are not specific spectral features, and the $\gamma$-ray spectra is determined from theoretical calculations. As time increases ($\gtrapprox$ 1 hour), we begin to see the emergence of the specific spectral lines of longer-lived nuclei. These come to dominate the spectra by $\sim$ 1 day. These spectral lines can be easily paired with specific nuclei, for instance with the decay of 208Tl as proposed by Vassh et al. (2024), which is responsible for the 2.6 MeV line, and 140La, which is responsible for the 2.0, 3.1, and 3.3 MeV lines seen at 1 week, although radiation transport may be necessary to determine the observability of the lines at earlier times. We expect that additional $\gamma$-rays are emitted during prompt fission emission, and will provide a larger high energy tail to the distribution (see e.g. Wang et al. (2020)), but this is not included in the calculation, which only considers the $\beta$-decay spectra.

The neutron spectrum is the most featureless of the spectra. Initially, this spectrum is much more significant than that of the $\gamma$-ray spectra, but shrinks to become comparable to the $\gamma$-ray spectra by $\sim$10-100 s and then significantly smaller by later times. We attribute this to the fact that while neutron capture is still occurring ( $\lessapprox$ 100 s), the nuclei which are created are farther from equilibrium and therefore have a much larger number of $\beta$-delayed neutron emission. Indeed, as time increases, the $\gamma$-ray spectrum grows to match and surpass the neutron spectrum (at around $\sim$ 1 second); however, the two remain in close proximity until neutron capture ends ($\sim$ 100 seconds). Beyond 100 seconds, neutrons emitted can be attributed to decays of the byproducts of fission, and like the $\gamma$-rays, it is expected that additional high-energy neutrons are emitted during prompt fission emission, but this is not included in the calculation.

The solid lines of Fig. 3a indicate the average particle energy (Eq. 5) for each particle as a function of time. In general, we note the average energies of the electron and neutrino are high for the first several seconds: $\sim 5$ MeV until $~1$ second, then 1.5-2.5 MeV until $\sim 10$ seconds, not dropping below 1 MeV until neutron capture ends ($\sim 100$ seconds). In contrast, the average energies of neutrons and $\gamma$-rays remain lower the entire time ($\lessapprox$ 1 MeV). However, for all species, the average energy is substantially higher than what would be expected from thermal equilibrium (dotted lines). To further illustrate the highly non-thermal nature of the particle spectra, we compare in Fig. 3b the spectral distribution at $t=10$ seconds to what would be expected from a thermal distribution. The difference is striking, showing a much harder spectrum than if thermal equilibrium was assumed.

Fig. 4a gives the total energy flux (Eq. 4) as a function of time, with the relative contributions from each channel shown in 4b. We note that the distribution of available energy is neither equal between channels nor constant in time. The largest fraction of energy goes to the neutrinos (between $\sim 40\%$ and $50\%$), while the rest is mainly divided between the $\gamma$-ray and electron emission. At early times, ($\lessapprox 1$ s), the contribution to $\gamma$-rays is small, and this is the only time period in which the contribution to neutrons is significant, though most of the difference in $\gamma$-rays is allocated to an increased electron contribution.

As the energy distributions of the emission spectra are an important input for the associated kilonova modeling, both the non-thermal distributions and the unequal allocation of available energy between channels may have significant impact. Full radiative transport is necessary to estimate the timescale for the emission spectra to become thermalized and to understand the effects on the resultant light curve. We emphasize that our calculations have no radiative transport, but are intended to provide self-consistent modeling of particle spectra which can be incorporated into the corresponding transport codes to enable improved understanding of kilonova light curves.

The particle emission spectra underlying the results presented are publicly available via Zenodo (Gross and others, 2025).

IV Observability of Signals

While these spectral emission profiles are an important component of the modeling of kilonova with significant impact on the final light curve, we can also consider the possibility of the direct observation of these emission profiles. Emitted electrons and neutrons will interact with the surrounding medium and will not be directly observable. However, emitted neutrinos will escape, and if they are produced in sufficient numbers at early times, they may be observable. $\gamma$-rays are initially opaque to observations, but at later times, the opacity will lower and some $\gamma$-rays will escape. The timescale for which $\gamma$-rays in the ejecta become less opaque is estimated be of order hours to days (see e.g. Guttman et al. (2024), which estimate the a timescale of $\sim$ 1 day for a GW180817-like ejecta). Full transport calculations are required to model the ejecta to determine the expected spectral signatures at early times, but at later times, if the emitted $\gamma$-rays are produced in sufficient number, specific emission lines may be directly observable. In this section, we provide an estimate of the neutrino and $\gamma$-ray flux resulting from $r$-process emission.

We assume that the emitted particles are a uniform, point-like emission. We can calculate the total incident flux from a total ejecta mass M at a distance D from the source:

[TABLE]

where $\Phi^{(i)}(E,t)$ is the total flux of species (i) rays emitted per second per unit energy (per unit abundance), similar to that in Eq. 4, $M_{n}$ is the average nucleon mass ($1.67\times 10^{-24}$ g), $M_{ejecta}$ is the assumed total mass which undergoes nucleosynthesis, and $D$ is the distance to the source. We nominally take $D$ = 15 kpc and $M_{ejecta}=0.01M_{\odot}$ as conservative but reasonable values for an event which occurs in our galaxy.

IV.1 Neutrinos

We display our results for the incident neutrino flux in Fig. 5a at varying times over the first 100 seconds. To estimate the possibility of detection, we consider only the first 10 seconds of emission, which constitutes 95% of the total neutrino emission. Following the approach of An et al. (2023), we consider inverse $\beta$-decay as the detection channel. The number of events in a Cherenkov detector can be estimated as (e.g. An et al. (2023)):

[TABLE]

where $\bar{\sigma}_{\text{IBD}}\approx 9.5\times 10^{-42}\ cm^{2}\times[\bar{E}^{\nu}/(10\text{ MeV})]^{2}$, and $N_{p}\sim 1.3\times 10^{34}\times[M_{T}/(200\text{ kton})]$, $\bar{E}^{\nu}$ is the average neutrino energy, $M_{T}$ is the detector size, and $F^{\nu}_{total}(t)$ is the total number of incident neutrinos by time t, which we estimate from the integration of Eq. 10 to be $F^{\nu}_{total}(10\ s)=1.54\times 10^{8}$ cm*-2*. We adopt $\bar{E}^{\nu}$ = 3 MeV and $M_{T}$ = 260 kton, which corresponds to the approximate size of Hyper-Kamiokande. With these values, we find $\mathcal{N}\sim 2.2$ observed events. Therefore, $r$-process nucleosynthesis events in our galaxy may be observable, though we emphasize that this only an estimate which ignores many important components of a thorough analysis, such as oscillations, detector response, etc.

IV.2 $\gamma$-Rays

The $\gamma$-rays will undergo spectral broadening due to the expansion of the ejecta. Assuming a spherical mass distribution for the ejecta with homologous expansion at velocity $v$, the resultant distribution of a spectral line of energy $E_{0}$ can be written as (See §A for a derivation):

[TABLE]

where $N_{total}$ is the number of particles emitted per unit time, and the range of possible energies E (for a line of energy $E_{0}$) are $E_{0}/(1+v)\leq E\leq E_{0}/(1-v)$, where $v$ is in units of $c$. Taking a nominal value of $v=0.1\ c$, and concatenating the above spectral broadening with Eq. 10, we calculate and display the incident $\gamma$-ray flux for an ejecta of mass $0.01M_{\odot}$ at distance 15 kpc in Fig. 5b. We compare this to the predicted observational capabilities of existing and future $\gamma$-ray detectors. Our results suggest that $\gamma$-rays from events within our galaxy should be observable for weeks to months. We note that this is a conservative estimate: the ejecta mass was taken to be only $0.01M_{\odot}$, while some models suggest that the ejection mass could be as high as $0.1M_{\odot}$. Our results are comparable with those of other studies which perform detailed radiation transport (see e.g. Jacobi et al. (2025), Bernuzzi et al. (2025)), though we highlight the difference in focus —while these studies focus on the observability at early times out to large distances (40 Mpc), we wish to highlight the direct observability of the emission profiles for more modest distances (15 kpc) at later times after the opacity has lowered.

V Conclusions

We have presented a first-principles calculation of the emission spectra of electrons, $\gamma$-rays, neutrinos, and neutrons as a function of energy and time for $r$-process nucleosynthesis. The spectra that we calculate have several unique properties: they are decisively non-thermal, with average energies far higher ($\gtrapprox$ 1 MeV) than would be expected from a thermal distribution ($\sim$ keV). We also find that the partitioning of energy is not even between channels as a function of time. For the majority of the evolution, the largest fraction of the energy is partitioned to neutrinos, with smaller but significant amounts going to electrons and $\gamma$-rays, and very little to neutrons, except at early times. Our microscopic calculations provide the spectral detail necessary for higher fidelity kilonova modeling, and their inclusion in future radiative transport models will enable a more precise interpretation of observations.

We also estimate the observability of the produced $\gamma$-ray and neutrino signal, finding that $r$-process activity resultant from GRBs may be directly observable within our galaxy for significant nucleosynthesis $(M\gtrapprox 0.01M_{\odot}$). The possibility of direct observation offers an exciting opportunity to bypass the degeneracies of kilonova modeling and probe nucleosynthesis through other messengers. In particular, spectral features of the $\gamma$-rays can be tied directly to the specific nuclei (such as the 2.6 MeV line produced by the $\beta$-decay of 208Tl). More work is necessary to identify and characterize these observational features in a variety of nucleosynthetic scenarios.

LANL is operated by Triad National Security, LLC, for the National Nuclear Security Administration of U.S. Department of Energy (Contract No. 89233218CNA000001). Research presented in this article was supported by the Laboratory Directed Research and Development program of Los Alamos National Laboratory under project number 20230052ER.

Appendix A Doppler Broadening of Line Profiles for Ejecta

We model our ejecta as a spherically symmetric mass distribution which is homologously expanding. The ejecta is emitting photons at a rate $\dot{\epsilon}$ per unit mass per second with wavelength $\lambda$. Homology implies that the velocity profile is proportional to the radius, and we can use this to write the mass profile in terms of the velocity:

[TABLE]

where $m_{ex}$ is the total ejection mass and $v_{ex}$ is the velocity of the edge of the distribution.

The Doppler shift due to expansion can be written as:

[TABLE]

where $v_{z}$ is in units of $c$ and gives the velocity component along the line of sight. We note that since we have homologous expansion, the value of $v_{z}$ is solely based on the distance from the center of the source along this line of sight. Therefore, we consider the circular slice which is centered along the line of sight with velocity component $v_{z}$, and can calculate the differential mass of this slice:

[TABLE]

The total photon emission can then be expressed as:

[TABLE]

Noting that we have $dv_{z}/d\lambda=-1/\lambda_{0}$, and using Eq. A2, we find:

[TABLE]

where the bounds on this distribution are $\lambda_{0}(1-v)\leq\lambda\leq\lambda_{0}(1+v)$. Translating this into the energy regime, we find:

[TABLE]

where $N_{total}=m_{ex}\dot{\epsilon}$ is the total emission rate for the entire ejecta, and the bounds on the distribution are now $E_{0}/(1+v)\leq E\leq E_{0}/(1-v)$.