Gamma-rays from kilonova: a potential probe of r-process nucleosynthesis

TL;DR

This paper explores gamma-ray emissions from neutron star merger remnants as a direct method to measure r-process nucleosynthesis yields, using simulations to predict spectra and assess detectability with future gamma-ray observatories.

Contribution

It introduces detailed 3D radiative transport simulations of gamma-ray spectra from kilonova remnants, highlighting their potential to reveal nucleosynthesis yields and distinguish fission models.

Findings

Gamma-ray spectra can serve as fingerprints of ejecta properties.

Detectability of gamma rays is feasible with next-generation detectors.

Fission fragment signatures can differentiate nucleosynthesis models.

Abstract

The mergers of compact binaries with at least one neutron star component are the potential leading sites of the production and ejection of $r$-process elements. Discoveries of galactic binary pulsars, short gamma-ray bursts, and gravitational-wave detections have all been constraining the rate of these events, while the gravitational wave plus broadband electromagnetic coverage of binary neutron star merger (GW170817) has also placed constraints on the properties (mass and composition) of the merger ejecta. But uncertainties and ambiguities in modeling the optical and infrared emission make it difficult to definitively measure the distribution of heavy isotopes in these mergers. In contrast, gamma rays emitted in the decay of these neutron-rich ejecta may provide a more direct measurement of the yields. We calculate the gamma production in remnants of neutron star mergers, considering two epochs: a kilonova epoch, lasting about two weeks, and a much later epoch of tens and hundreds of thousands of years after the merger. For the kilonova epoch, when the expanding ejecta is still only partially transparent to gamma radiation, we use 3D radiative transport simulations to produce the spectra. We show that the gamma-ray spectra associated with beta- and alpha-decay provide a fingerprint of the ejecta properties and, for a sufficiently nearby remnant, may be detectable, even for old remnants. We compare our gamma spectra with the potential detection limits of next generation detectors, including the Lunar Occultation Explorer (LOX), the All-sky Medium Energy Gamma-ray Observatory (AMEGO), and the Compton Spectrometer and Imager (COSI). We show that fission models can be discriminated via the presence of short-lived fission fragments in the remnant spectra.

Tables

Table 1. Table 1: Models Summary.

	Weak	Strong	Fission
Model	$r$-process	$r$-process	Model	Shape
	(spherical)	(toroidal)
S1	$Y_e=0.4$	—	—	Spherical
S2	$Y_e=0.3$	—	—	Spherical
As	—	$Y_e=0.05$	Symmetric	Torus
Ak	—	$Y_e=0.05$	KT∗	Torus
S1As	$Y_e=0.4$	$Y_e=0.05$	Symmetric	Sphere+torus
S1Ak	$Y_e=0.4$	$Y_e=0.05$	KT	Sphere+torus
S2As	$Y_e=0.3$	$Y_e=0.05$	Symmetric	Sphere+torus
S2Ak	$Y_e=0.3$	$Y_e=0.05$	KT	Sphere+torus

Table 2. Table 2: Some of the Isotopes with Bright Lines which Produce the Spectral Peaks Visible in Figure  4 .

Models	Time Range	Line Energy [keV]	Isotope	$T_{1/2}$
S1	1 hr–1 day	1384	${}_{38}{}^{92}\text{Sr}$	2.611 hr
	2 hr–2 days	934	${}_{39}{}^{92}\text{Y}$	3.54 hr
	12 hr–2 days	658	${}_{41}{}^{97}\text{Nb}$∗	72.1 minutes
	$>8$days	756	${}_{40}{}^{95}\text{Zr}$	64.0 days
S2	6 hr–4 days	743	${}_{51}{}^{128}\text{Sb}$	9.05 hr
	6 hr–4 days	754	${}_{51}{}^{128}\text{Sb}$	9.05 hr
	$>4$days	364	${}_{53}{}^{131}\text{I}$	8.02 days
	$>4$days	80.2	${}_{53}{}^{131}\text{I}$	8.02 days
	$>4$days	29.8	${}_{53}{}^{131}\text{I}$	8.02 days
	$>4$days	2002	${}_{50}{}^{125}\text{Sn}$	9.64 days
As, Ak	12 hr–1 day	77.4	${}_{84}{}^{197}\text{Pt}$	19.9 hr
	$>2$days	81.0	${}_{54}{}^{133}\text{Xe}$	5.25 days
As	hr	847	${}_{53}{}^{134}\text{I}$	52.5 minutes
	hr	884	${}_{53}{}^{134}\text{I}$	52.5 minutes
	12 hr–2 days	249	${}_{54}{}^{135}\text{Xe}$	9.14 hr
Ak	$>10$days	145	${}_{58}{}^{141}\text{Ce}$	32.5 days
	$>10$days	1596	${}_{57}{}^{140}\text{La}$∗	1.67 days

Table 3. Table 3: Possible Influential γ 𝛾 \gamma -ray Emitters and their Long-lived Populating Ancestors as Found by Network Calculations on a 10 ky Observational Timescale.

Isotope	$T_{1/2}$	Mass Range [M⊙]	Ancestor(s)	$T_{1/2}$	Ancestor Mass Range [M⊙]	Line Energy [keV]	Flux [ph s-1 cm-2]
${}_{95}{}^{241}\text{Am}$	432.6 yr	(1 – 10)$\times {10}^{-9}$	${}_{96}{}^{245}\text{Cm}$	8423 yr	$\left(2-20\right)\times {10}^{-8}$	59.5409	$\left(2-20\right)\times {10}^{-8}$
${}_{95}{}^{243}\text{Am}$	7364 yr	(1 – 10)$\times {10}^{-8}$	Self			74.66	$\left(2-30\right)\times {10}^{-8}$
${}_{95}{}^{246}\text{Am}$	39 minutes	(1 – 10)$\times {10}^{-17}$	${}_{96}{}^{250}\text{Cm}$	8300 yr	$\left(8-80\right)\times {10}^{-9}$	679.2	$\left(3-30\right)\times {10}^{-9}$
						756	$\left(7-70\right)\times {10}^{-10}$
${}_{83}{}^{213}\text{Bi}$	45.59 minutes	$\left(5-60\right)\times {10}^{-17}$	${}_{90}{}^{229}\text{Th}$	7880 yr	$\left(5-50\right)\times {10}^{-9}$	440.45	$\left(4-40\right)\times {10}^{-9}$
${}_{83}{}^{214}\text{Bi}$	19.9 minutes	(1 – 10)$\times {10}^{-17}$	${}_{88}{}^{226}\text{Ra}$	1600 yr	$\left(5-50\right)\times {10}^{-10}$	609.32	$\left(3-40\right)\times {10}^{-9}$
			${}_{90}{}^{230}\text{Th}$	75.4 ky	$\left(2-20\right)\times {10}^{-8}$	1120.294	(1 – 10)$\times {10}^{-9}$
						1764.491	(1 – 10)$\times {10}^{-9}$
${}_{97}{}^{250}\text{Bk}$	3.21 hr	$\left(3-30\right)\times {10}^{-17}$	${}_{96}{}^{250}\text{Cm}$	8300 yr	$\left(8-80\right)\times {10}^{-9}$	1028.654	(1 – 10)$\times {10}^{-10}$
${}_{96}{}^{245}\text{Cm}$	8423 yr	$\left(2-20\right)\times {10}^{-8}$	Self			99.5232	(1 – 10)$\times {10}^{-8}$
						103.741	$\left(2-20\right)\times {10}^{-8}$
						117.2322	$\left(4-40\right)\times {10}^{-9}$
						175.01	$\left(5-50\right)\times {10}^{-9}$
${}_{53}{}^{134}\text{I}$	52.5 minutes		K–T fission			847.025
						884.09
						1072.55
${}_{57}{}^{140}\text{La}$	1.68 days		K–T fission			487.021
						1596.21
${}_{93}{}^{239}\text{Np}$	2.36 days	(1 – 10)$\times {10}^{-14}$	${}_{95}{}^{243}\text{Am}$	7364 yr	(1 – 10)$\times {10}^{-8}$	99.5232	$\left(4-50\right)\times {10}^{-9}$
						103.741	$\left(7-80\right)\times {10}^{-9}$
						106.123	$\left(8-90\right)\times {10}^{-9}$
						277.599	$\left(5-50\right)\times {10}^{-9}$
${}_{82}{}^{214}\text{Pb}$	27.06 minutes	$\left(2-20\right)\times {10}^{-17}$	${}_{88}{}^{226}\text{Ra}$	1600 yr	$\left(5-50\right)\times {10}^{-10}$	241.995	$\left(5-60\right)\times {10}^{-10}$
			${}_{90}{}^{230}\text{Th}$	75.4 ky	$\left(2-20\right)\times {10}^{-8}$	295.2228	(1 – 10)$\times {10}^{-9}$
						351.9321	$\left(3-30\right)\times {10}^{-9}$
${}_{51}{}^{125}\text{Sb}$	2.76 yr	${10}^{-18}$ – ${10}^{-11}$	Symm fission			427.874	$8\times {10}^{-15}$ – ${10}^{-7}$
						463.365	$3\times {10}^{-15}$ – $4\times {10}^{-8}$
						600.597	${10}^{-15}$ – ${10}^{-7}$
						635.95	${10}^{-15}$ – ${10}^{-8}$
${}_{51}{}^{126}\text{Sb}$	12.35 days	$\left(.3-30\right)\times {10}^{-14}$	${}_{50}{}^{126}\text{Sn}$	230 ky	$\left(.2-20\right)\times {10}^{-7}$	414.7	$\left(.2-20\right)\times {10}^{-8}$
						666.5	$\left(.2-20\right)\times {10}^{-8}$
						695.0	$\left(.2-20\right)\times {10}^{-8}$
						720.7	(.1 – 10)$\times {10}^{-8}$
${}_{50}{}^{125}\text{Sn}$	9.64 days		Symm fission			822.48
						915.55
						1067.1
						1089.15
${}_{50}{}^{126}\text{Sn}$	230 ky	$\left(.2-20\right)\times {10}^{-7}$	K–T, Self			64.281	$\left(.2-20\right)\times {10}^{-9}$
						86.938	$\left(.2-20\right)\times {10}^{-9}$
						87.567	$\left(.7-70\right)\times {10}^{-9}$

Table 4. Table 4: Same as Table 3 , Except with a 100 ky Observational Timescale.

Isotope	$T_{1/2}$	Mass Range [M⊙]	Ancestor(s)	$T_{1/2}$	Ancestor Mass Range [M⊙]	Line Energy [keV]	Flux [ph s-1 cm-2]
${}_{95}{}^{243}\text{Am}$	7364 yr	(1 – 10)$\times {10}^{-9}$	${}_{96}{}^{247}\text{Cm}$	15.6 My	(4 – 40)$\times {10}^{-9}$	74.66	(3 – 30)$\times {10}^{-9}$
${}_{83}{}^{213}\text{Bi}$	45.59 minutes	(7 – 70)$\times {10}^{-18}$	${}_{92}{}^{233}\text{U}$	159.2 ky	(1 – 10)$\times {10}^{-8}$	440.45	(5 – 50)$\times {10}^{-10}$
${}_{83}{}^{214}\text{Bi}$	19.9 minutes	(6 – 60)$\times {10}^{-18}$				609.32	(2 – 20)$\times {10}^{-9}$
			${}_{90}{}^{230}\text{Th}$	75.4 ky	(1 – 10)$\times {10}^{-8}$	1120.294	(5 – 60)$\times {10}^{-10}$
			${}_{92}{}^{234}\text{U}$	245.5 ky	(2 – 20)$\times {10}^{-8}$	1238.122	(2 – 20)$\times {10}^{-10}$
						1764.491	(5 – 60)$\times {10}^{-10}$
${}_{93}{}^{239}\text{Np}$	2.36 days	(1 – 10)$\times {10}^{-15}$				99.5232	(6 – 60)$\times {10}^{-10}$
			${}_{95}{}^{243}\text{Am}$	7364 yr		103.741	(6 – 60)$\times {10}^{-10}$
			${}_{96}{}^{247}\text{Cm}$	15.6 My	(2 – 20)$\times {10}^{-8}$	106.123	(6 – 60)$\times {10}^{-10}$
						277.599	(6 – 70)$\times {10}^{-10}$
${}_{91}{}^{233}\text{Pa}$	26.98 days	(3 – 30)$\times {10}^{-15}$	${}_{93}{}^{237}\text{Np}$	2.144 My	(8 – 90)$\times {10}^{-8}$	300.129	(6 – 70)$\times {10}^{-11}$
						311.904	(4 – 40)$\times {10}^{-10}$
						340.476	(4 – 50)$\times {10}^{-11}$
${}_{82}{}^{214}\text{Pb}$	27.06 minutes	(8 – 80)$\times {10}^{-18}$	${}_{90}{}^{230}\text{Th}$	75.4 ky	(1 – 10)$\times {10}^{-8}$	241.995	(3 – 30)$\times {10}^{-10}$
			${}_{92}{}^{234}\text{U}$	245.5 ky	(2 – 20)$\times {10}^{-8}$	295.2228	(7 – 70)$\times {10}^{-10}$
						351.9321	(1 – 10)$\times {10}^{-9}$
${}_{51}{}^{126}\text{Sb}$	12.35 days	(.2 – 20)$\times {10}^{-14}$	${}_{50}{}^{126}\text{Sn}$	230 ky	(.1 – 10)$\times {10}^{-7}$	414.7	(.1 – 10)$\times {10}^{-8}$
						666.5	(.2 – 10)$\times {10}^{-8}$
						695.0	(.2 – 10)$\times {10}^{-8}$
						720.7	(.8 – 80)$\times {10}^{-9}$
${}_{50}{}^{126}\text{Sn}$	230 ky	(.1 – 10)$\times {10}^{-7}$	K–T, Self			64.281	(.1 – 10)$\times {10}^{-9}$
						86.938	(.1 – 10)$\times {10}^{-9}$
						87.567	(.6 – 50)$\times {10}^{-9}$

Equations

$$S(E,t)=N_A\sum _{\mathrm{iso}}{\lambda }_{\mathrm{iso}}{Y}_{\mathrm{iso}}(t)\left({\mathrm{RP}}_{\mathrm{iso}}(E)+\sum _{\gamma (\mathrm{iso})}{I}_{\gamma }^{\mathrm{iso}}{E}_{\gamma }^{\mathrm{iso}}\delta (E-E_{\mathrm{iso}})\right),$$

$${\sigma }_{\mathrm{PE}}=n_{\mathrm{abs}}{\sigma }_{\mathrm{abs}}={\rho }_{\mathrm{eje}}/(m_{\mathrm{proton}}\bar{a}){\sigma }_{\mathrm{Fe}}$$

$$r_{\mathrm{shock}}=(E_{\exp }/{\rho }_{\mathrm{CSM}}{)}^{2/5}{t}^{2/5}$$

$$v_{\mathrm{shock}}=d{r}_{\mathrm{shock}}/dt\approx (E_{\exp }/{\rho }_{\mathrm{CSM}}{)}^{2/5}{t}^{-3/5}.$$

$$v_{\mathrm{shock}}=v_{\mathrm{ejecta}}{m}_{\mathrm{ejecta}}/(m_{\mathrm{ejecta}}+4\pi r_{\mathrm{shock}}^3{\rho }_{\mathrm{CSM}})$$

$$r_{\mathrm{shock}}=(v_{\mathrm{ejecta}}{m}_{\mathrm{ejecta}}/\pi {\rho }_{\mathrm{CSM}}{)}^{1/4}{t}^{1/4}.$$

Full text

Neutron stars (1108) — Explosive nucleosynthesis (503) — R-process (1324) — Gamma-ray sources (633) — Gamma-ray lines (631)

††software: Maverick (Hungerford et al., 2003, 2005), Numpy (van der Walt et al., 2011), PRISM (Mumpower et al., 2017, 2018)

1 Introduction

Over four decades ago, Lattimer & Schramm (1974) proposed that rapid decompression of neutron-rich matter from a tidally disrupted neutron star could account for the $r$-process production of the universe. Proving this point requires demonstrating that the rate of neutron star mergers is sufficiently high and that the cumulative nucleosynthetic yield is plentiful, given the merger rate, and furthermore, produces the solar-like distribution in proper agreement with $r$-process enriched metal-poor stars (Sneden et al., 1996; Beers & Christlieb, 2005; Hansen et al., 2018; Ji & Frebel, 2018). Rates of these mergers from theoretical (e.g. Fryer et al., 1999b; Dominik et al., 2012) and observed binary pulsars (e.g. Kalogera et al., 2004; Chen & Holz, 2013), and gamma-ray bursts (GRBs; e.g. Paul, 2018) span a wide range, arguing that they produce between $<1$% and 100% of the $r$-process (Côté et al., 2017). Theoretical rates are uncertain because binary population synthesis models suffer from large uncertainties in stellar evolution (e.g. stellar radii and shell sizes), binary evolution (e.g. common envelope evolution and mass transfer) and supernova (e.g. neutron star kicks) properties. Observations, on the other hand, are prone to bias (e.g. determining the completeness of the observed sample). The gravitational-wave detection of GW170817 provided an independent observational constraint, arguing for a sufficiently high rate that, with yields currently predicted by simulations, mergers could produce most, if not all, of the $r$-process elements (Côté et al., 2018; Rosswog et al., 2018). While ongoing gravitational-wave detections are refining these rate estimates, studies from the perspective of galactic chemical evolution indicate that several $r$-process sites were operating in the early universe (Hotokezaka et al., 2018; Côté et al., 2019; Simonetti et al., 2019).

With the merger rate increasingly constrained, the viability of mergers as an $r$-process source depends more upon the uncertainties in the amount and composition of the merger ejecta. The ejecta from the merger occurs during the initial tidal disruption, as well as at late times, as the debris accretes onto the merged core (Dessart et al., 2009; Perego et al., 2014; Martin et al., 2015; Siegel & Metzger, 2017). Theory predicts a range of ejecta masses about ${10^{-3}-10^{-2}M_{\odot}}$ (Korobkin et al., 2012; Bauswein et al., 2013; Hotokezaka et al., 2013; Endrizzi et al., 2016; Radice et al., 2016; Sekiguchi et al., 2016). While the tidal (or dynamical) ejecta is believed to be neutron rich, and hence has been argued to produce a composition that is "robust" in $r$-process elements, the neutron fraction can be reset by neutrinos, producing everything from iron peak elements to the heavy $r$-process (Fernández & Metzger, 2013; Wanajo et al., 2014; Fernández et al., 2015; Just et al., 2015). To truly understand the yields from neutron star mergers, we must understand both the ejecta composition and their amount.

Optical, ultraviolet and infrared electromagnetic counterparts of neutron star mergers provide one venue for inferring the nature of these ejecta (Li & Paczyński, 1998; Piran, 2005; Metzger et al., 2010). Specifically, astronomers argued for both “red” (produced from ejecta with heavy $r$-process) and “blue” (ejecta with atomic masses only up to and including the second $r$-process peak) components (Metzger & Berger, 2012). The Lanthanides synthesized as part of the heavy $r$-process have many lines in the ultraviolet, optical, and near-infrared wavelength bands, driving the emission to the mid-infrared. These ejecta produce the "red" component in the emission seen in many calculations (e.g. Kasen et al., 2013; Tanaka & Hotokezaka, 2013; Fontes et al., 2015). If the late-time ejecta is less neutron rich to the point that there are insufficient neutrons to produce the heavy $r$-process elements, the ejecta will generate a bright, short-lived blue transient (e.g. Metzger & Berger, 2012; Barnes & Kasen, 2013; Perego et al., 2014; Wollaeger et al., 2018).

Prior to GW170817, astronomers had to make a series of assumptions to probe the ejecta properties of neutron star mergers. First, they established a connection between short GRBs and neutron star mergers by observing that offset distributions (Fong & Berger, 2013) of short GRBs match predictions of neutron star populations (Bloom et al., 1999; Fryer et al., 1999b). They then assumed that deviations in the power-law decay of GRB afterglows could arise from the emergence of radioactive emission from the ejecta. A number of kilonova candidates were identified (Perley et al., 2009; Tanvir et al., 2013; Fong et al., 2014; Jin et al., 2015, 2016; Lamb et al., 2019). However, observing such components is difficult, because the kilonova light-curve signal must be separated from much brighter background of the GRB afterglow, and shocks in the afterglow may produce bumps in the optical/infrared that can be mistaken for kilonova light (Kasliwal et al., 2017). If the infrared excess has a corresponding X-ray flare, it is more likely to be caused by shock interactions with the inhomogeneities in the circumstellar medium rather than powered by the ejecta radioactivity. With GW170817, the ejecta emission –kilonova– was observed unambiguously for the first time, providing a first direct probe of this phenomenon. The combined strong blue and red components of this merger seemed to fit the models predicted for both dynamical/tidal and late-time wind/disk ejecta, allowing to infer the masses of individual components.

But recent analysis of the GW170817 kilonova spectra has made it clear that uncertainties in the model would make it difficult to make concrete claims about the amount and composition of the ejecta. Overviews of the analyses from different groups show a broad range of inferred ejecta masses (Côté et al., 2018; Ji et al., 2019). Much of the uncertainty in light-curve calculation comes from the modeling of opacities and their incorporation into transport codes (Kasen et al., 2013; Tanaka & Hotokezaka, 2013; Fontes et al., 2017). The methods used to calculate the opacities, the number of levels (and lines) considered, and the methods to combine these opacities in an expanding medium all can affect the light curve (Fontes et al., 2019). However, the uncertainties in the ejecta properties (density, velocity, and composition distributions) and morphology produce even larger uncertainties (Grossman et al., 2014; Wollaeger et al., 2018). Thus, even with the pristine data from GW170817, it is difficult to determine the ejecta masses to better than an order of magnitude. Other effects also muddle the interpretation and analysis of the kilonova emission. For example, the flux (especially in the optical and ultraviolet) can vary dramatically with the viewing angle (see, e.g. Fernández et al., 2017; Wollaeger et al., 2018). All of these studies assumed that radioactive decay powers the emission, but additional energy sources (a pulsar or emission from accretion onto the compact remnant) can also impact the light curve (Wollaeger et al., 2019).

With all of these uncertainties, it is difficult to estimate accurate ejecta masses based solely on the broadband light curves. Obtaining detailed abundances is even more challenging. It is possible that spectral features can provide evidence of the composition, and there are hints that the GW170817 must have ejected at least some light $r$-process elements (Pian et al., 2017; Watson et al., 2019), but obtaining detailed yields requires full, time-dependent and out-of-equilibrium opacities. It is also possible to constrain ejecta masses via radioactive heating (Piran et al., 2014; Rosswog et al., 2018), but this approach is only partially successful, as heating is expected to look similar for many initial conditions (Lippuner & Roberts, 2015; Eichler et al., 2019).

As with 56Ni yields in thermonuclear supernovae (Churazov et al., 2014; The & Burrows, 2014) and 56Ni and 44Ti yields in core-collapse supernovae (Hungerford et al., 2005; Grefenstette et al., 2014, 2017), a more direct measurement of the yields can be obtained by observing the photons from the decay of radioactive nuclei in the ejecta. In this paper, we study the potential of measurements of decay photons to probe the nucleosynthetic yields and nuclear physics in neutron star mergers. We focus our efforts on the study of $\gamma$-rays produced by the nuclear decay of neutron-rich nuclei. In a pioneering work on this subject, Hotokezaka et al. (2016) calculated the $\gamma$-ray signal from kilonova ejecta and found that it would be detectable out to $\sim 3-10\,{\rm Mpc}$ with current detectors. However, their work was done without modeling $\gamma$-ray transport, which can significantly redistribute emission to lower energies, impairing detectability. Recently, Li (2019) constructed a semianalytic model of the radioactive $\gamma$-ray emission from kilonovae powered by nuclear decays. An earlier study by Janiuk (2014) suggested detection of X-ray emission from iron-group isotopes synthesized in central engines of GRBs.

Although the main peak flux of $\gamma$-rays happens at early times, the emission continues for more than a hundred thousand years after the merger. Therefore, it is possible that there is a nearby kilonova remnant that can be observed. In a complementary study, Wu et al. (2019) consider prospects of finding such neutron star merger remnants in the Milky Way galaxy. In earlier work, Qian et al. (1998) concluded that sensitive $\gamma$-ray detector can observe lines from a few long-lived heavy radioactive isotopes decaying in supernova remnants, in particular $\prescript{125}{51}{\text{Sb}}$, $\prescript{126}{50}{\text{Sn}}$, $\prescript{137}{55}{\text{Cs}}$, $\prescript{144}{58}{\text{Ce}}$, $\prescript{155}{63}{\text{Eu}}$ and $\prescript{194}{76}{\text{Os}}$. Subsequently, Ripley et al. (2014) investigated search prospects for both supernova and neutron star merger remnants within the Galactic plane using the NuGRID and LOFT X-ray observatories. It was found that $>10^{2}$ overabundance is required to detect the lines of the most promising isotope, $\prescript{126}{50}{\text{Sn}}$. Fuller et al. (2019) argued that thermal positron production at the initial stage of kilonova explosion could generate strong 511 keV annihilation line signature which might help with identifying such remnants. We further explore possible $\gamma$-ray emission from the remnants, using detailed $r$-process nucleosynthesis calculations and models for ejecta deceleration in the interstellar medium.

Section 2 introduces our method, including the ejecta morphologies, detailed nucleosynthesis models, and $\gamma$-ray source calculation. Section 3 discusses early-time $\gamma$-ray signatures of kilonova, following the fully 3D transport of the emitted $\gamma$-rays through the ejecta. In section 4, we calculate the properties of neutron star merger remnants. We conclude with a comparison with upcoming $\gamma$-ray missions.

2 Gamma Rays from $r$-process Yields

The neutron-rich ejecta from neutron star mergers are expected to produce a wide range of elements from the iron peak to third $r$-process peak and beyond. Gamma-ray signatures would therefore depend not only on the neutron richness, but also on thermodynamic history and the morphology of the ejecta, which affect this history. The ejecta neutron richness ranges from extremely high in the dynamical part produced in the process of tidal disruption of the neutron stars, to the medium-richness outflows from the accretion disk (Janiuk, 2014; Siegel & Metzger, 2017; Miller et al., 2019) up to the much more symmetric ejecta in the outflows from central merger product (Perego et al., 2014; Martin et al., 2015).

The extent of corresponding variability of the $\gamma$-ray production is demonstrated in Figure 1 which shows it as a function of time for different thermodynamic conditions and nuclear mass models. Here, the $\gamma$-radiation rate $\varepsilon_{\gamma}(t)$ is normalized to the power-law decay fit (Metzger & Berger, 2012; Korobkin et al., 2012) $\varepsilon_{0}(t)=2\times 10^{10}\,{\rm erg}\,{\rm g}^{-1}{\rm s}^{-1}\;t_{d}^{-1.3}$, where the time $t_{d}$ is measured in days. Such a plot allows forthe easy identification of epochs of active $\gamma$-ray emission for a given model.

In Figure 1, variability due to nuclear mass model is shown with a color band, while different colors represent different initial hydrodynamic conditions (see details below). The two top panels show composition outcomes from moderately neutron-rich ejecta ($Y_{e}=0.4$ and $Y_{e}=0.3$), and exhibit several orders of sensitivity to the hydrodynamics as opposed to only about one order of magnitude sensitivity to nuclear mass model. On the other hand, models starting with extremely neutron-rich conditions ($Y_{e}=0.05$) exhibit very little variation in $\gamma$-ray production (bottom panel). This is because nucleosynthesis in this regime is robust and governed by fission recycling much more than by hydrodynamics (Korobkin et al., 2012; Holmbeck et al., 2019).

The power law with fractional power index emerges from a multitude of $\beta^{-}$-decaying isotopes (Metzger & Berger, 2012; Hotokezaka et al., 2017). At late times, it breaks into individual peaks produced by individual radioactive decay chains. A chain generating peak at time $t$ starts with an isotope having mean lifetime $\tau\approx t$ which does not necessarily produce $\gamma$-rays: rather, a different $\beta^{-}$-decaying isotope downstream with much shorter lifetime may be responsible. Decay chains \prescript{126}{50}{\text{Sn}}$$\to$$\prescript{126}{51}{\text{Sb}} and \prescript{237}{93}{\text{Np}}$$\to\dots$$\prescript{213}{83}{\text{Bi}} shown schematically at the bottom panel of Fig. 1 give examples of such scenario.

To keep parameter space manageable, we only explore a few thermodynamic trajectories representative of different ejecta types and nucleosynthesis models. We adopt a two-component model motivated by numerical simulations of neutron star mergers (Rosswog et al., 2014) and used in our two-dimensional study of kilonova light curves (Wollaeger et al., 2018). As in Wollaeger et al. (2018), the two components are neutron-rich "dynamical ejecta" and lighter $r$-process-producing “wind”.

The morphology of the dynamical ejecta is derived from model A in Rosswog et al. (2014, see their Table 1), which was computed in the neutron star merger simulation and followed by the subsequent expansion of the ejecta up to homology. We rescale its mass for the best fit to the GW170817 kilonova (as in our models for this event presented in Evans et al., 2017; Tanvir et al., 2017; Troja et al., 2017). For the secondary, less neutron-rich and wind-like outflow, we pick an analytic spherically symmetric background solution as introduced in Wollaeger et al. (2018). Dynamical ejecta is rescaled to have mass ${m_{\rm dyn}=0.0065M_{\odot}}$ and median expansion velocity ${v_{\rm dyn}=0.2c}$, while the wind outflow is heavier and slower: ${m_{\rm wind}=0.03M_{\odot}}$, ${v_{\rm wind}=0.08c}$ (Tanvir et al., 2017). The morphologies of the two components in our models are depicted in Figure 2.

To produce the nucleosynthetic composition for our model components, we use the Portable Routines for Integrated nucleoSynthesis Modeling (PRISM) reaction network, most recently used in Côté et al. (2018), Vassh et al. (2019), Sprouse et al. (2019). This network uses state-of-the-art nuclear physics inputs (e.g. Mumpower et al., 2017, 2018; Möller et al., 2019), including a consistent treatment of capture rates as well as neutron-induced and $\beta$-delayed fission using the theoretical framework of Kawano et al. (2008, 2016), Mumpower et al. (2016). This framework combines together various statistical nuclear model inputs such as nuclear level densities, $\gamma$-ray strength functions and optical potentials to produce well-tested predictions (Spyrou et al., 2016; Yokoyama et al., 2019) for nucleosynthesis calculation. Variations in nuclear binding energies proceed as in Mumpower et al. (2015) with the statistical model inputs held fixed.

The time evolution of the abundances $Y_{\rm iso}(t)$ is used to calculate the detailed $\gamma$-ray source. The source represents finely binned spectrum, based on the total spectrum $S(E,t)$, which in turn is computed using abundances of the decaying isotopes, their known $\gamma$-radiation lines and the spectrum of the continuum component (if present), ${\rm RP}_{\rm iso}(E)$:

[TABLE]

where the first and second sums are over all decaying isotopes and $\gamma$-radiation lines for each isotope, respectively. Each $\gamma$-radiation line is characterized by energy, $E^{\rm iso}_{\gamma}$, and absolute intensity, $I^{\rm iso}_{\gamma}$, per single decay. The spectrum ${\rm RP}_{\rm iso}(E)$ is normalized to unity, ${\int{\rm RP}_{\rm iso}(E)dE=1}$. Here, we use recent data, provided by the Evaluated Nuclear Reaction Data Library ENDF/B-VIII.0 111https://www-nds.iaea.org/public/download-endf/ENDF-B-VIII.0/ library (Brown et al., 2018). Finally, $\lambda_{\rm iso}$ is the decay rate of the isotope, and $N_{A}$ is Avogadro’s number. Total $\gamma$-ray energy production $\varepsilon_{\gamma}(t)$ is easily obtained by integrating the spectrum $S(E,t)$ over energy.

While Figure 1 was computed with a large number of nuclear mass models (25 models), three hydrodynamics conditions for each model and two fission prescriptions for the neutron-rich case, in the rest of the paper we focus on just four representative yield distributions. The nucleosynthesis is computed with parameterized trajectories (an exponential plus power-law decay described in Lippuner & Roberts, 2015) and self-consistent nuclear reheating. Composition of dynamical ejecta is calculated assuming initial entropy ${s=10\,k_{\rm B}/{\rm baryon}}$, electron fraction ${Y_{e}=0.05}$, and expansion timescale ${\tau=10\ {\rm ms}}$. We explore sensitivity to nuclear physics by using two different fission prescriptions: a symmetric splitting (following Mumpower et al., 2018), and the fission fragment distributions of Kodama & Takahashi (1975). Composition for both of the "winds" is computed with expansion timescale ${\tau=100\ {\rm ms}}$, initial entropy ${s=50\,k_{\rm B}/{\rm baryon}}$, and electron fractions ${Y_{e}=0.4}$ for “wind 1” and ${Y_{e}=0.3}$ for “wind 2”. This is similar to the basic compositions used in Wollaeger et al. (2018).

The resulting yield distributions one day after the merger are shown in Figure 3. Two distributions for the low-$Y_{e}$ dynamical ejecta (red, As and black, Ak) represent strong $r$-process between the second and third peak, computed with two different types of fission model, as previously described. The medium-$Y_{e}$ wind component (green, S2) spans the range from first to the second $r$-process peaks, while the high-$Y_{e}$ component (blue, S1) only produces the first $r$-process peak. These four uniform-composition models are selected to represent dominant peak contribution. Models S1 and S2 have spherically symmetric morphology (“S”) and correspond to the yields with $Y_{e}=0.4,0.3$, respectively. Models As and Ak have morphology of model "A" from Rosswog et al. (2014) and correspond to the strong $r$-process production with symmetric split and Kodama–Takahashi fission models, respectively. Superimposing these models, we additionally construct four two-component models. Our models are summarized in Table 1.

3 Gamma Rays from Kilonova Transients

The $\gamma$-ray emission is strongest at early times (first 10 days) when it emerges from the expanding ejecta. This is the so-called kilonova epoch. For a nearby event, the gravitational-wave and follow-up electromagnetic detections of this event will provide exact localization, allowing dedicated $\gamma$-ray follow-up of the kilonova. Initially, most of the emitted gammas are trapped in the flow and the escape of this emission requires transport calculations. The transport necessarily distorts $\gamma$-ray source spectrum: it broadens every line and absorbs or redistributes energy. We use the $\gamma$-ray production spectra for four representative yields calculated in the previous Section 2 to source $\gamma$-rays in the transport code, and perform 3D transport simulations on a material background of the two morphologies as described above.

For this optically thick transport regime, we use the Monte Carlo $\gamma$-ray transport code Maverick described in Hungerford et al. (2003, 2005). In the context of 56Ni decay in thermonuclear supernovae, this code has been verified in a code comparison effort against most major codes in the community (Milne et al., 2004). Maverick assumes the material properties (density and composition) are in steady state for each time slice. Average escape time of gamma-ray packets is $<25$% of the age of the explosion for all time slices considered, so this steady-state assumption is reasonable. The ejecta is followed assuming a homologous expansion and then mapped into a 3D ($50^{3}$) grid for the transport.

We assume that the source spectrum is proportional to the mass in each zone. We use luminosity-weighted packets, so the number of Monte Carlo packets in each zone is also proportional to the mass. The packets sample the energy of the $\gamma$-rays based on our emission spectrum and are binned into 2000 energy groups ranging from 5 keV to 20 MeV.

The $\gamma$-ray opacity includes components from Compton scattering, photoelectric absorption, and pair production absorption. It is dominated by Compton scattering above roughly 100-300 keV. Photoelectric absorption becomes important below 100-300 keV, depending on the $Z$ of the absorbing material. Compton scattering is treated by sampling the outgoing photon properties (energy and angle) from the complete Klein–Nishina scattering kernel in the free electron limit. The electron density in each zone is contributed by electrons from the wind component as well as electrons from the dynamic ejecta component.

The photoelectric absorption opacity ($\sigma_{\rm PE}$) is represented as an effective absorption as follows:

[TABLE]

where the number of absorbers ($n_{\rm abs}$) is set to the density of the ejecta ($\rho_{\rm ejecta}$) divided by the average atomic mass ($\bar{A}$) and the proton mass ($m_{\rm proton}$). Here, the ejecta can include both wind and dynamical ejecta components. The cross section per absorber is taken to be the relatively well-known cross section of iron ($\sigma_{\rm Fe}$). This simplifying assumption for the cross section can lead to errors in our opacity, especially below 100 keV, where it dominates the opacity as the photoelectric cross section scales as roughly the proton fraction to the fourth power, but it provides a rough estimate for the opacity. However, features below 100 keV should be taken with some caution.

With this physics, we use Maverick to calculate the escape fraction and energy of the Monte Carlo packets. These packets are tallied into a spectrum that has 250 logarithmically spaced energy bins from 3 keV to 20 MeV. Figure 4 shows the resulting spectra for both one- and two-component models. There are distinct differences in the $\gamma$-ray signal between all of our models in the first few hours, which persist to late times.

All models show more high-energy spectrum at early times. Models Ak and As lose about one order of magnitude in brightness between the first hour and the first day, while S1 and S2 gain approximately the same amount, catching up and becoming dominant emitters compared to As and Ak around the spectral peak at 1 MeV. Ak model (as well as S2Ak and S1AK) has a distinct enhancement at the first two hours in the range of high $\gamma$-ray energies $>5$ MeV. At early times, the line broadening is noticeably blue-shifted due to the photosphere approaching the observer, while at later times (after one day or so) the broadened lines become much more symmetric. Dynamical ejecta models show much less features; this is not so much due to the morphology expanding twice as fast on average as because there are many isotopes contributing and blending to form a pseudo-continuous spectra. Nevertheless, lines from certain radioactive nuclides such as $\prescript{197}{84}{\text{Pt}}$ and $\prescript{133}{54}{\text{Xe}}$ can be clearly identified. Isotopes $\prescript{140}{57}{\text{La}}$ and $\prescript{141}{58}{\text{Ce}}$ are only prominent for the Ak model, while $\prescript{135}{54}{\text{Xe}}$ and $\prescript{134}{53}{\text{I}}$ emerge for the As model, for which a different fission prescription was used. This gives a hint that potentially a correct fission model can be decided from the observation. Table 2 lists the properties of the most prominent $\gamma$-ray producing isotopes for each of the models.

An important feature distinguishing kilonova transients in $\gamma$-rays is that the lines of individual decaying nuclides become prominent only on the timescale comparable to their mean lifetime. This is unlike the optical or infrared signal, which is affected by the entire yield at all times. Bright $\gamma$-ray emitting isotopes take turns to emerge in the spectrum, allowing the potential to trace evolution of the composition in real time. However, this effect is mitigated by the long integration time, even for the most sensitive detectors.

In summary, $\gamma$-ray observations will be able to determine whether the ejecta originates from the electron poor or electron rich initial conditions. However, the differences between fission models Ak versus As are very small and will be difficult to detect. In models with mixed electron fractions and multiple components, it will be difficult to determine the exact yield (only that the material is mixed and not dominated by a low- or high-electron fraction abundance. After 10 days, the emission has dropped by two orders of magnitude, becoming increasingly difficult to detect.

4 Gamma Rays from Kilonova Remnants

The detection of old kilonova remnants provides an alternate observational prospect to constraining the nucleosynthesis in neutron star mergers. Although the rate of neutron star mergers is about three orders of magnitude lower that that of supernovae, given the fact that a few hundred supernova remnants have been discovered, it is not unreasonable to assume that kilonova remnants younger than 100 kyr can be found in our neighborhood of the Milky Way (Wu et al., 2019). If a relatively young remnant exists close to the Earth, we may be able to detect it and probe the yields of the merger. The $\gamma$-ray spectrum of a kilonova remnant would consist of multiple lines generated by long-lived residual nuclides from the $r$-process, providing unique perspective on its nuclear physics. As the remnant decelerates, line broadening is less important (Piran et al., 2013), producing individually identifiable lines of specific radionuclides. This can be particularly helpful for discriminating between various $r$-process scenarios. This is true even for very dilute interstellar medium in the Galactic halo. In this section, we study both the remnant evolution to determine velocities and spatial sizes of kilonova remnants and the expected $\gamma$-ray signals, comparing the results from two fission models.

4.1 Kilonova Remnant Evolution and Properties

An explosive remnant (whether it be a supernova or kilonova) passes through four evolutionary phases: free expansion, Sedov–Taylor, snowplow and merger with interstellar medium. The free expansion phase is assumed to last until the ejecta sweeps up mass comparable to itself. During this phase, the expectation is that the ejecta is expanding without decelerating. The velocity of the shock ($v_{\rm shock}$) is a constant and the radius of the shock ($r_{\rm shock}$) increases with time ($t$) linearly.

When radiative cooling is slow compared with the shock evolution, the Sedov–Taylor similarity solution (Taylor, 1941, 1950; Sedov, 1946) is used to model the shock evolution. This similarity solution can be derived through simple dimensional analysis: $[E_{\rm exp}/\rho_{\rm CSM}]$ have units of $({\rm g}\;{\rm cm}^{2}\;{\rm s}^{-2})/({\rm g}\;{\rm cm}^{-3})={\rm cm}^{5}\;{\rm s}^{-2}$. With these units, we can derive the shock position,

[TABLE]

where $E_{\rm exp}$ is the explosion energy, $\rho_{\rm CSM}$ is the density of the circumstellar medium which is, for massive stars, the stellar wind, and for neutron stars, the interstellar medium. For a blast wave moving through a constant-density medium, the radius increases as time to the $2/5$ power. The corresponding shock velocity ($\rm v_{\rm shock}$) is

[TABLE]

This phase continues until radiative cooling becomes faster than the evolution of the shock. At this point, the shock evolves through a snowplow phase where the evolution is dictated by momentum conservation. In this phase, the remnant velocity ($v_{\rm shock}$) is

[TABLE]

where $m_{\rm ejecta}$ is the ejecta mass and $v_{\rm ejecta}$ is the ejecta velocity. At late times, the ejecta mass can be neglected in the denominator and the radius as a function of time is

[TABLE]

To determine how well these simple analytic estimates match the properties of the remnant, we have modeled the ejecta expansion numerically in 1D to late times. For the purposes of this study, two properties are most crucial: the velocity distribution of the radioactive ejecta to get line broadening and the extent of the remnant. Our numerical model uses a 1D Lagrangian hydrodynamics code initially designed for supernovae (Fryer et al., 1999a), but modified (using a simple $\gamma=5/3$ equation of state) to follow the ejecta out to large distances. With this code, we calculate several models with varying ejecta masses, velocities and densities of the circumstellar medium. Because of the kicks imparted on neutron stars at birth, the merger can happen far off of the Galactic plane where the density of the surrounding medium is low, spanning a large range: ${10^{-4}-10^{2}\;{\rm cm^{-3}}}$ (Wiggins et al., 2018).

Before we discuss the full suite of results, consider the evolution of the explosion better. Figure 5 shows the velocity profile of a shock from a kilonova explosion with $0.01\,M_{\odot}$ of ejecta, $6\times 10^{49}$ erg of energy, and an interstellar medium (ISM) density of roughly $0.001\,cm^{-3}$ for times ranging from 10 days to 500 yr. It takes over 100 yr for the shock to sweep up a mass equal to the ejecta mass fully transition to the Sedov–Taylor phase where the velocity decreases with the radius to the 3/2 power. Note that there is a transition region where the shock decelerates but not as quickly as expected with Sedov–Taylor.

We have constructed models of the kilonova remnant, coupling the four phases of the remnant evolution to determine both the remnant size and velocity (Fig. 6) as a function of time. Within a factor of two or so, the late-time properties of these remnants ($>10^{4}$ yr) are not very different from supernovae. Although the velocities are higher in kilonovae, the lower ejecta masses mean that the kilonova remnant decelerates faster than normal supernovae. We also expect radiative cooling to dominate sooner when the kilonova remnant transitions from the Sedov to the snowplow phase at earlier times, leading to more rapid deceleration after roughly 10,000 yr. At $10^{4}$ yr, the kilonova remnant is expanding at between one hundred and a few thousand ${\rm km\,s^{-1}}$ and at $10^{5}$ yr may already have decelerated to the sound speed of the ISM (tens of ${\rm km\,s^{-1}}$) and still expanding at one hundred ${\rm km\,s^{-1}}$.

An interesting feature of kilonova remnants is the rapid evolution to the Sedov phase. Whereas supernova remnants are free streaming for the first 100–1,000 yr (depending on the density of the interstellar medium), kilonovae enter this phase between 0.25 and 100 yr. During the Sedov phase, a reverse shock is produced that heats the ejecta, driving strong radio emission, detectable sometimes within a year or a few years from the outburst (Piran et al., 2013).

4.2 Remnant Gamma Rays

It is estimated that a few neutron star merger remnants exist in our neighborhood of the Milky Way with ages below a 100 ky (Ripley et al., 2014; Wu et al., 2019). As follows from the previous section, for remnants with 10-100 kyr ages, the ejecta velocities are likely to lie between 100 and $3,000{\rm km\,s^{-1}}$ and the remnant size lies between 5 and 300 pc. A remnant 3,kpc away from the Earth will have angular size $0.3^{\circ}-6^{\circ}$. In this section, we review the expected $\gamma$-ray spectra from these remnants as a function of composition. In the same way as atomic spectra can be used to infer composition, these decay spectra can be used as fingerprints of the yields. At these ages, velocity broadening will not significantly alter the line signals.

For neutron-poor ejecta with a composition peaking near the first $r$-process peak, many of the isotopes have already decayed by 10–100 kyr. But a few isotopes, 99Tc, 126Sn, 126Sb, and 129I, contribute to the $\gamma$-ray spectra with energy spanning from roughly 30 eV to a few MeV. Figure 8 shows the $\gamma$-rays for both $Y_{e}=0.3$ and $Y_{e}=0.4$ ejecta. The decay timescales for these isotopes are long (more than 100 kyr) and the signal at 10 kyr is not so different than the signal at 100 kyr.

If we instead focus on the neutron-rich ejecta, both the spectra and the physics are more involved, as complex decay pathways may arise leading to nonintuitive $\gamma$-ray emitters. Figure 7 shows one such example where 237Np, the long-lived ancestor with $T_{1/2}\sim 2.1\times 10^{6}$ yr, decays into $\gamma$-ray producing 213Bi, a nucleus whose half-life is roughly 45 minutes.

At sufficiently high neutron richness, the yields are less sensitive to the exact neutron fraction, but more sensitive to the nuclear physics such as the fission model, reinforcing the need for improved nuclear physics modeling for the $r$-process (Horowitz et al., 2019). Figure 9 shows the spectra at 10 and 100 kyr for our dynamical ejecta with two different fission models. As with the atomic spectra, there is a forest of decay lines. Because of the forest of lines, velocity broadening can merge lines and we include plots at the low and high end (100, 3000 km s*-1*) of our remnant velocities. With expected energy resolutions, it may still be possible to distinguish between the yields of different nuclear physics models.

Network calculations of neutron-rich ejecta suggest 126Sb, 128Sb, 214Bi, 214Pb, 243Am, 246Am, 245Cm, and 250Bk are the dominant isotopes contributing to the spectra on the observational timescale of 10 and 100 kyr. These isotopes are the result of decays from long-lived ancestors that set the observational timescale. We summarize possible influential $\gamma$-ray emitters and their long-lived ancestors in Tables 3 and 4. Some of these isotopes have been studied as potential indicators of the $r$-process in previous works Qian et al. (1998); Ripley et al. (2014); Wu et al. (2019) as they can be found by surveying the half-lives in the nuclear chart. Other isotopes emerge as prominent emitters despite having short half-lives as decay products downstream of long-lived populating ancestors. This demonstrates how comprehensive treatment of nucleosynthesis is significant when catching $\gamma$-ray emitters.

Fission of actinides produces several short-lived isotopes which leave distinct mark on the $\gamma$-ray spectrum. Table 3 lists two such isotopes, $\prescript{134}{53}{\text{I}}$ and $\prescript{140}{57}{\text{La}}$ being signature for the symmetric-split fission model, and two other, $\prescript{125}{51}{\text{Sb}}$ and $\prescript{125}{50}{\text{Sn}}$ as signature isotopes for the Kodama–Takahashi fission. Both Tables 3 and 4 list $\prescript{126}{50}{\text{Sn}}$ as an isotope which can be either synthesized from the very beginning, or as a signature fission product of the Kodama–Takahashi fission prescription. If a remnant were to be discovered nearby, this hints to the exciting possibility to observe $\gamma$-ray lines of short-lived fragments of actinide fission, such as those actively studied in experiment, e.g., the CAlifornium Rare Isotope Breeder Upgrade (CARIBU) facility at Argonne National Laboratory (Marley et al., 2013; Van Schelt et al., 2013). Observing such lines will allow to reason about models of nuclear fission.

Tables 3 and 4 show that for a remnant at distance of 3 kpc away none of the isotopes produce line flux in excess of $10^{-6}$ photons s*-1* cm*-2*. Sensitivity of at least $10^{-7}$ photons s*-1* cm*-2* is required to detect such a remnant. For a remnant that is 100 kyr old, the line flux is an order of magnitude less. Ripley et al. (2014) suggested blind search for neutron star merger remnants within the Galactic plane. However, an older remnant is possible in the halo, away from the Galactic plane. At 100 kyr in the halo with low density of interstellar medium the remnant would have the radius of $60-200$ pc (see Figure 6) and the angular size of $1^{\circ}-6^{\circ}$ (see also Wu et al., 2019).

5 Discussion: Detectability Prospects

To assess the detectability of kilonova remnants, we compare our results to a number of existing detectors as well as three different proposed satellite missions: Lunar Occultation Explorer (LOX), the Compton Spectrometer and Imager (COSI), and the All-sky Medium Energy Gamma-ray Observatory (AMEGO). Each of these missions has different strengths and weaknesses in observing kilonovae and their remnants and we review each of them here. The most recent proposal for COSI is a COSI-SMEX mission and its sensitivity should lie between the COSI-X and GRX proposals (A. Zoglauer, private communication). We take the latest sensitivity curves from the AMEGO (Rando, 2017) and LOX (R. Miller, private communcation) collaborations. The LOX satellite is focused on the 0.1-few MeV range and its predicted sensitivities in this range are nearly two orders of magnitude lower than AMEGO, but it has only roughly 10% energy resolution.

Figure 10 shows the transient signal, integrated over 1 Ms starting from 1 hour from the moment of merger for our ejecta models at 3 Mpc. There are roughly 100 galaxies within 3 Mpc and it is likely that the transient will be localized quickly, allowing a nearly instantaneous observation. Integrating over this period provides a reasonable estimate of the observed flux for these transients. At the 10-100 keV range, NuSTAR might be able to detect the signal from some mergers with 3 Mpc. For COSI-SMEX and AMEGO, telescopes that can be pointed, we assume a steady 1 Ms observation. Even assuming continuous observation by COSI-SMEX or AMEGO, such an event will be difficult to detect. The LOX satellite should be able to detect a merger at 3 Mpc, but a merger at 10 Mpc will be just at the observing threshold. We have also included the sensitivity of the balloon-based concentrator concept (Shirazi et al., 2018).

Probably more exciting is the possibility of a nearby, old kilonova remnant. Figure 11 shows the detectability of a 10 kyr remnant 3 kpc from the Earth with sensitivities assuming 1 year of directed time. With only a handful of merger remnants younger than 100 kyr in the entire Milky Way, the odds of remnant this close to the Earth is less than 1%. A nearby remnant will likely be in a denser interstellar medium (e.g. $1\,{\rm cm}^{-3}$), slowly expanding ($\sim 500{\rm km}\,{\rm s}^{-1}$) and with a small extent ($1\,{\rm pc}\approx 0.5^{\circ}-1^{\circ}$). Since these remnants vary slowly with time, sky surveys can be mined to look for this data (to achieve 1 yr of directed time will require multiple years of telescope runtime). With the high-energy resolution of COSI-SMEX it would be possible to distinguish the individual features in this signal, but would require a remnant less distant than 3 kpc. If we do not have to correct for the fact that nearby remnants can not be treated as point sources, LOX will be able to observe these objects up to 10 kpc and identify some of the largest features enough to distinguish between our two fission rate results. Sensitivity and energy resolutions of AMEGO (Kierans, 2019) lies in between COSI-SMEX and LOX.

6 Summary

We have studied the $\gamma$-rays that arise primarily from the $\beta^{-}$ and $\alpha$-decays of radioactive isotopes in the kilonova ejecta of neutron star mergers222More work can be done to better understand the $\gamma$-ray signal and we have not included all sources of gamma-rays (e.g. fission-induced $\gamma$-ray emission).. Comprehensive nucleosynthesis network modeling was used to generate the yields. We compared the signatures of a few representative compositions reflecting different types of neutron star merger ejecta and variations in nuclear physics. We studied both the transient kilonova phase that requires a new merger event, and an old neutron star merger remnant phase which could be identified in in our galaxy via its gamma emission.

We further analyzed detectability prospects with upcoming telescope proposals for several existing (NuSTAR, Chandra) and future (COSI, AMEGO, LOX) $\gamma$-ray missions. For the kilonova epoch, a neutron star merger event must happen within 10 Mpc to be detectable (see Figure 10). If detected, it may be possible to distinguish broadened lines from very few individual isotopes (see Figure 4 and Table 2) and reason about the ejecta composition. A nearby old merger remnant presents another potential detection possibility (see Figure 11). For the remnant to be detectable with the high-energy resolution of COSI-SMEX it must be located closer than 3 kpc. LOX and AMEGO will be able to observe such remnants up to 10 kpc and identify some of the emitting isotopes. We demonstrated that in such case it is possible to discriminate between our two fission prescriptions (see Tables 3 and 4). Thus, if the spectrum can be measured, short-lived $\gamma$-emitters from fragments of actinide fission will allow to reason about fission models as well as potentially infer the amount of parent actinides.

We stress that a number of physics effects could alter the signals presented in this work. For example, synchrotron radiation may generate a background in the same energy range as our nuclear decay lines and residual thermal positron annihilation can produce a strong 511 keV feature (Fuller et al., 2019). Furthermore, we have not included in our model the $\gamma$-rays from nuclear fission and nuclear isomeric states that are populated in radioactive decays, which may influence the observed spectrum. An example of an isomer that may have observable consequences for kilonovae is the $\prescript{97}{41}{\text{Nb}}$ metastable state at 743 keV which has a 97.9% $\gamma$ branch to the ground state. The de-excitation of this isomer, which is populated by the $\beta^{-}$-decay of $\prescript{97}{40}{\text{Zr}}$, may produce an observable feature near this energy beginning around 12 hours post-merger.

Just as with supernovae, $\gamma$-ray observations from the decay of radioactive nuclei require nearby events with a rate much lower than those achieved with optical and infrared observations. But, as with supernovae, this study, along with the work of Hotokezaka et al. (2016) and Wu et al. (2019), shows the unique potential of $\gamma$-rays to probe the details of nucleosynthesis in neutron star mergers (including nuclear physics), thereby ensuring their importance in understanding these powerful transients.

Acknowledgments

We thank Andreas Zoglauer for providing data on upcoming gamma-ray detectors. O.K. is grateful to Brian Metzger for providing the original motivation for this work. We thank Almudena Arcones, Ryan Wollaeger, Kenta Hotokezaka, Masaomi Tanaka, Xilu Wang, Marius Eichler, Jonah Miller, Stephan Rosswog, Nicole Lloyd-Ronning, and Kei Davis for valuable discussions. We also thank the anonymous referee for constructive suggestions. This work was supported by the US Department of Energy through the Los Alamos National Laboratory. Los Alamos National Laboratory 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 20190021DR. The work on the concentrator detector was supported by NASA grant NNX17AC85G. All LANL calculations were performed on LANL Institutional Computing resources. This work has benefited from support by the National Science Foundation under grant No. PHY-1430152 (JINA Center for the Evolution of the Elements).

.