Gamma-ray Diagnostics of r-process Nucleosynthesis in the Remnants of Galactic Binary Neutron-Star Mergers

We perform a full nuclear-network numerical calculation of the $r$-process nuclei in binary neutron-star mergers (NSMs), with the aim of estimating $\gamma$-ray emissions from the remnants of Galactic NSMs up to $10^6$ years old. The nucleosynthesis calculation of 4,070 nuclei is adopted to provide the elemental composition ratios of nuclei with an electron fraction $Y_{\rm e}$ between 0.10 and 0.45 . The decay processes of 3,237 unstable nuclei are simulated to extract the $\gamma$-ray spectra. As a result, the NSMs have different spectral color in $\gamma$-ray band from various other astronomical objects at less than $10^5$ years old. In addition, we propose a new line-diagnostic method for $Y_{\rm e}$ that uses the line ratios of either $^{137{\rm m}}$Ba/$^{85}$K or $^{243}$Am/$^{60{\rm m}}$Co, which become larger than unity for young and old $r$-process sites, respectively, with a low $Y_{\rm e}$ environment. From an estimation of the distance limit for $\gamma$-ray observations as a function of the age, the high sensitivity in the sub-MeV band, at approximately $10^{-9}$ photons s$^{-1}$ cm$^{-2}$ or $10^{-15}$ erg s$^{-1}$ cm$^{-2}$, is required to cover all the NSM remnants in our Galaxy if we assume that the population of NSMs by \citet{2019ApJ...880...23W}. A $\gamma$-ray survey with sensitivities of $10^{-8}$--$10^{-7}$ photons s$^{-1}$ cm$^{-2}$ or $10^{-14}$--$10^{-13}$ erg s$^{-1}$ cm$^{-2}$ in the 70--4000 keV band is expected to find emissions from at least one NSM remnant under the assumption of NSM rate of 30 Myr$^{-1}$. The feasibility of $\gamma$-ray missions to observe Galactic NSMs are also studied.


Introduction
Elements heavier than Bi exist in our universe, but their origin remains a mystery. Most cosmic isotopes heavier than the iron group are expected to be created by the rapid-neutron capture process, also known as the r-process (Burbidge et al. 1957;Cameron 1957;Cowan et al. 1991;Wanajo & Ishimaru 2006;Arnould et al. 2007;Qian & Wasserburg 2007), but the actual nucleosynthesis sites capable of achieving such neutron-rich environments remain a matter of debate. Before the discovery of binary neutron-star mergers (NSMs) observed as gravitational wave objects like GW170817 (Abbott et al. 2017), NSMs were considered to be more promising as rprocess nucleosynthesis sites than other primary candidates, such as core-collapse supernovae (SNe) because NSMs could achieve more neutron-rich (lower electron fraction Y e ) environments (Lattimer & Schramm 1974;Metzger et al. 2010;Wanajo et al. 2011). The event rate of NSMs is much lower than that of SNe, but the yield of r-process nuclei in one event is expected to be very high (Hotokezaka et al. 2015;Wallner et al. 2015). Observational evidence of the existence of r-process nuclei has already been obtained by infrared observations of kilonovae (also called macronovae or r-process novae) in some short gamma-ray bursts, such as GRB 130603B (Tanvir et al. 2013) and the gravitational wave event GW170817 (Villar et al. 2017). However, the infrared radiation from NSMs is, in principle, the result of indirect emissions from unstable r-process nuclei, and any hint of elements heavier than the lanthanoids is still missing from the infrared information. Given that the nuclear levels of nuclei are in the megaelectronvolt energy range, the gamma-rays from r-process nuclei should be the best probe for searching for r-process sites in the universe.
According to theoretical estimates of the gamma-ray flux from binary NSMs (Hotokezaka et al. 2016), the gamma-ray radiation immediately following a merging event is very dim at about 10 −8 -10 −7 photons s −1 cm −2 keV −1 , even at an extremely close distance d of 3 Mpc. This flux is comparable to or below what the sensitivities of current and near-future megaelectronvolt missions can detect. The precise measurements of photon energies are, in principle, rather difficult in the megaelectronvolt band, where Compton scattering dominates over the photon-absorption process. Therefore, the ability to detect gamma-rays from NSMs by an immediate follow-up observation (a Target-of-opportunity observation; ToO) would be limited by the sensitivity of the gamma-ray instruments. Instead, a non-ToO observation of gamma-rays from long-lived nuclei in NSMs would be an alternative way to survey r-process sites, and this has been proposed by Wu et al. (2019) and Wang et al. (2020). The gamma-ray luminosity from nuclei with long lifetimes, on the order of 10 3 -10 6 yr, becomes much lower than that from short-lived nuclei, but if we limit the survey area within our Galaxy (d  10 kpc), then the gammaray flux in non-ToO observations is expected to become comparable to that required for ToO observations. Therefore, non-ToO observations should provide more sensitive gammaray surveys of NSMs because the exposure time (the accumulation time of signals) is not limited as it is in ToO observations. Another benefit of performing a non-ToO survey is the better identification of gamma-ray lines; we expect the effect of Doppler broadening to be smaller for older NSM remnants than for very young NSMs.
Here, we focus on the non-ToO survey of gamma-rays from r-process nuclei in a possible Galactic NSM remnant. In this paper, we estimate gamma-ray emissions from Galactic NSM remnants in an older age range than in previous work (Hotokezaka et al. 2016;Wang et al. 2020) by using nuclearnetwork numerical calculations with a complete nuclear database. This paper also aims to provide gamma-ray diagnostic methods for NSMs, showing the required sensitivities for future gamma-ray observatories. In our study, we assume that gamma-ray instruments have a wider field-of-view (FOV) than the object size of the NSM remnants, which are larger than early NSMs in a ToO observation. We also assume that the instruments accumulate all of the gamma-ray emissions from the NSM remnants, even though the nuclei may mix with the circumstellar medium (CSM) during the evolution of the remnants.
The rest of this paper is organized as follows. In Section 2, we summarize our environments and procedures for the nuclear-network numerical calculation and show the results for gamma-ray emissions from NSM remnants. In Section 3, we present the gamma-ray diagnostics, which utilize spectral color to identify NSM remnants and provide the line properties for estimating the age t and Y e . In Section 4, we discuss the survey distance and coverage in our Galaxy permitted by the instrument sensitivities, the corresponding limitation of the NSM rate in our Galaxy, and expectations for future missions.

Overview of Numerical Calculation
To estimate the gamma-ray emissions from binary NSM remnants of various ages, we performed a numerical simulation comprising the following three steps: (1) calculation of the mass distribution of r-process nuclei for NSMs at t = 1 yr, (2) calculation of the decay processes of unstable nuclei emitting gamma-rays, and (3) a simple calculation of the radiation transfer of gamma-rays from NSMs.
For the first step, we adopted the nucleosynthesis calculation for around 4070 nuclei performed by Fujimoto et al. (2007), which was cooled using the adiabatic expansion modeled from Freiburghaus et al. (1999) to provide the elemental composition ratios of nuclei for Y e = 0. 10, 0.15, 0.20, 0.25, 0.30, 0.35, 0.40, and 0.45. This estimation assumes that the initial environment has a temperature of 9 × 10 9 K, a radius of 100 km, entropy per baryon of 10 k B , where k B is the Boltzmann constant, and a velocity of 2 × 10 9 cm s −1 , along with the initial abundances of the 4070 nuclei in nuclear statistical equilibrium. As a result, the calculation provides the mass fractions at t = 1 yr evaluated with the nuclear reaction network (network A in Fujimoto et al. 2007), by using Y e = 0.10-0.45 in steps of 0.05. To set up the mass distribution of nuclei for the NSMs at t = 1 yr, we blended the nuclei with the mass fraction using the Y e provided in Wanajo et al. (2014). Specifically, the fractions are 4.54%, 4.85%, 14.6%, 29.7%, 10.3%, 25.1%, 10.5%, and 0.33% for Y e = 0. 10, 0.15, 0.20, 0.25, 0.30, 0.35, 0.40, and 0.45, respectively. Note that this Y e -fraction model by Wanajo et al. (2014) describes a slightly less neutron-rich environment than those by the recent dynamical-ejecta models after the kilonova observations of the gravitational event GW170817, such as four models in Kullmann et al. (2022) under two kinds of equations of state, density dependent 2 (DD2; Hempel & Schaffner-Bielich 2010;Typel et al. 2010), and SFHo (Steiner et al. 2013). In this paper, we adopted the first one by Wanajo et al. (2014) as a pessimistic case for the r-process site, but changing the Y e -fraction models does not change the conclusions from the gamma-ray spectra as tested in Section 3. Figure 1 shows the mass fraction of multiple nuclei at t = 1 yr generated in an NSM case, information that is given in the table of nuclides (neutron number N versus atomic number Z). Using the same data set, Figure 2 summarizes the distribution of nuclei with mass number A at t = 1 yr, showing the contributions of Y e . This plot demonstrates that the environment with lower Y e contributes to the generation of heavier elements.  For the second step, we simulated the decay processes of unstable nuclei, starting from the mass distribution at t = 1 yr calculated in the first step. We used the Decay Data File 2015 (DDF-2015; Katakura & Minato 2016) from the Japanese Evaluated Nuclear Data Library (Katakura 2012), which provides the decay profiles of 3237 nuclei up to Z = 104 (Rf). Here, we applied a correction to the gamma-ray information for 241 Am; this error was reported in our study and was corrected in the next version of the database. The originality of this study lies in the comprehensiveness of nuclei treated in the calculation. In the nuclear-decay calculation, we adopted the α decay, β − -decay, β + -decay, electron capture, isomeric transition, and gamma-decay processes. In our calculation, the internal conversion process is ignored, which emits soft X-rays and makes a negligible contribution to the gamma-ray band. The neutron-and proton-emission processes are also ignored because they contribute to the very early phase, which is beyond the scope of this study. The spontaneous fission process may occur on 257 Es, 256 Cf, 254 Cf, and 250 Cm, but its contribution is negligible. Thus, this process is also excluded from our calculation. In addition, we do not calculate the gamma-ray emission from secondary electrons (electrons from β-decay, δ-rays, and so on) after the decay of unstable nuclei.
To verify the calculations in the second step, we refer to Figure 3, which represents the relative abundances of nuclei at t = 10 3 , 10 6 , and 10 9 yr as a function of A. The distribution of A does not change dramatically after t ∼ 10 3 yr, except for the disappearance of the small dips at the magic numbers. The distribution for t > 10 6 yr becomes roughly consistent with the semi-empirical abundance distribution of cosmic r-process nuclei in Beer et al. (1997), which is close to the solar abundance distribution.
For the outputs of the second step, we obtain the gamma-ray flux F γ,i from the NSM at distance d. Using the nuclear gamma-ray intensity I γ,i of the ith element with the mass number A i , the mass M i , and the half-life T 1/2 , the F γ,i in a small time interval dt is described as Finally, the third step is to calculate the transfer of gammarays through the NSM ejecta. However, we omitted the detailed Monte Carlo calculation of the radiation transfer because the optical depth decreases rapidly after the merger event, by roughly µ t (Li 2019); within the scope of our study at t ? 1 yr, the optical depth is thin and negligible. Therefore, the degradation of the line profiles by Compton scattering is not included in our calculation, which would be dominant in only the very early phase. Note that the detailed calculations of the megaelectronvolt gamma-ray spectra from Galactic NSMs in the initial phase were performed by Wang et al. (2020Wang et al. ( , 2021 and gamma-rays from NSMs per Y e by Chen et al. (2021). In this step, we only apply the bulk Doppler-broadening effect caused by the expansion velocity v(t). The thermal Dopplerbroadening effect is ignored in this calculation because it is two or three orders of magnitude smaller than that from the expansion motion of the heavy elements in the A = 50-200 range. In reality, the line profile from the bulk Doppler effect becomes complicated due to the complex contributions of various velocity components, as has been observed in the X-ray lines from heavy elements in supernova (SN) remnants (Grefenstette et al. 2017;Kasuga et al. 2018). For simplicity, we applied the Gaussian distribution function for the line profile in the calculation. Of the various velocity elements in the remnant, we only applied single Gaussian broadening to the maximum velocity component, which we assume to be the forward shock motion. This was done to simulate the most robust case for considering the gamma-ray sensitivity. In our assumption, v(t) starts from the initial value v(0) = 0.3 c, where c is the speed of light, and evolves at a constant rate during the free expansion phase. During the Sedov-Taylor phase, v(t) evolves as v(t) ∝ t −(3/5) (Taylor 1950), and then as v(t) ∝ t −(0.7) during the pressure-driven snowplow phase (McKee & Figure 3. The relative abundance of 195 Pt in the NSM case at t = 10 3 , 10 6 , and 10 9 yr, shown in green, blue, and thick black, respectively, compared with the semiempirical abundance distribution of Beer et al. (1997).
Ostriker 1977). We assume that the free expansion, Sedov-Taylor, and probability density function phases end at t = 10, 4.7 × 10 4 , and 1.65 × 10 6 yr, respectively. The ages of these phase transitions may change by about one order of magnitude due to differences in density of the CSM, but this modification only affects the Doppler-broadening effect. It becomes negligible when compared with the typical energy resolutions of gamma-ray instruments for ages older than t ∼ 10 3 yr, the range that lies within the scope of this study. Note that even at t = 10 6 yr, v(t) approaches ∼20 km s −1 , which is about double the speed of sound for a typical CSM density of 0.01 cm −3 . The radius becomes ∼100 pc. Finally, we get the gamma-ray spectra for NSMs at t, accumulated from all of the r-process nuclei in the ejecta.

Gamma-Ray Emission and Evolution
From the numerical calculation described in Section 2.1, the gamma-ray spectra from t = 3-10 6 yr, under the assumption that the ejecta mass is M ej = 0.01 M e at d = 10 kpc, are summarized in Figures 4 and 5. As described in Section 1, we assume that all of the emissions from the NSM remnants are observable within the wider FOV; this is assumed to be larger than the object size, which becomes around 10 pc at t = 10 3 yr and expands into around 100 pc at t > 10 6 yr. The spectra contain many nuclear lines broadened by the Doppler effect. They appear to form a continuous spectrum in the early phase, but they become separated at ages older than 10 3 yr. Note that the gamma-ray data without the Doppler-broadening effect (the outputs from the second step of the calculation in Section 2.1) is provided as the numerical model for the XSPEC tool (Arnaud 1996) in the HEAsoft package (Appendix).
To identify the gamma lines in the spectra, we checked the most prominent lines in the gamma-ray spectra generated by a single Y e condition. Table 1 lists the brightest lines shown for the nuclei in each Y e . Roughly speaking, the bright lines seen for objects of a younger age lie in the higher-energy gamma-ray band of the spectrum. Further details of the diagnostics will be discussed in Section 3.

Spectral Color Changes of NSM Remnants
Using the energy spectra of the NSM remnants (Figures 4 and 5), we first checked the properties of the spectral shapes from the hard X-ray to the soft gamma-ray bands. As shown in the normalized spectra plotted in Figure 6, the energy spectra roughly evolve from hard to soft slopes. Gamma-ray emission decreases rapidly leaving the hard X-ray emission in old age, as is indicated in Table 1. This phenomenon of gamma-rays is equivalent to the Sargent law for β decay.
To see the evolution of the shape of the γ-ray spectra more quantitatively, we plotted the light curves of the gamma-ray flux in three bands: 70-200, 200-500, and 500-3000 keV, which cover multiple lines around 100 and 300 keV, and a prominent line around 700 keV, respectively. As indicated in the top panel of Figure 7, the flux in the higher-energy bands decreases more quickly than that in the low-energy bands. A decaying trend is also seen in the time dependency of the hardness ratio among these bands, as is indicated in the lower panel of Figure 7. The ratio drops dramatically at around 200-300 yr, indicating that the gamma-ray flux above the 500 keV band quickly decreases at this age. This phenomenon is primarily due to the decay of 125 Sb and 137m Ba listed in Table 1. Note that this result does not change even if we adopt other Y e -fraction models of DD2-125145, DD2-135135, SFHo-125145, and SFHo-135135 in Kullmann et al. (2022), as shown in Figure 7.
To compare the spectral shape of NSM remnants with other astronomical objects, we plotted the color-color diagrams in the hard X-ray band (10-500 keV) and in the hard X-ray to γray band (70-3000 keV), in the top and bottom of Figure 8, respectively. We divided the energy bandpass for these spectra into three ranges: 10-25, 25-70, and 70-500 keV for the hard X-ray band (top of Figure 8), and 70-500, 500-1,000, and 1000-3000 keV for the hard X-ray to gamma-ray band (Figure 8 bottom). Note that the divisions of the energy bands are defined so that they follow the energy bandpass of current gamma-ray instruments on board NuSTAR (Harrison et al. 2013), INTEGRAL (Winkler et al. 2003), and other observatories. For comparison, the spectral colors of other astronomical objects, calculated using the INTEGRAL catalog version 0043 8 , are also plotted in the same figures. In the hard X-ray band (the 10-500 keV band in the top of Figure 8), the spectra of NSM remnants older than t ∼ 1000 yr have spectral colors similar to those of SN remnants or active galactic nuclei, but NSM remnants younger than t ∼ 1000 yr can be distinguished from other known objects by their hard X-ray colors. In other words, the spectral color in the hard X-ray band below 500 keV is a good indicator of young NSM remnants. Furthermore, this differentiation from known objects becomes more prominent when we include the higher-energy band covering the megaelectronvolt portion of the spectrum, as is clearly indicated in the bottom of Figure 8. Note that this result does not change even if we adopt other Y e -fraction models of DD2-125145, DD2-135135, SFHo-125145, and SFHo-135135 in Kullmann et al. (2022), as shown in Figure 8. Therefore, NSM remnants have unique spectral colors in the hard X-ray to gamma-ray bands. This observation is one of the important conclusions drawn from our calculation. Note that the spectral models in the INTEGRAL catalog are simple enough that the colors of known objects in the gamma-ray band (bottom of Figure 8) are less scattered than those in the hard X-ray band (top of Figure 8). The spectral separation between NSM remnants and other objects in the bottom of Figure 8 does not change dramatically, even if we lower the low-energy threshold (70 keV in the bottom of Figure 8) to cover 20 keV, for example. However, it becomes worse if we set it higher so that everything up to a certain point; 200 keV, for example, is ignored. This implies that hard X-rays around 100 keV provide key information for distinguishing NSM remnants from other objects. Note that these results are based on the pure-nuclear gamma-rays from r-process nuclei in NSMs, and thus the synchrotron radiation from electrons that are accelerated by the shocks may contaminate the hard X-ray band for young remnants. Additionally, when taking actual observations, we must be careful to isolate the contamination of the hard X-ray spectrum that arises from other objects located behind the NSM, such as active galactic nuclei within the FOV.

Nuclear Line Emissions from Older NSM Remnants
For ages older than t > 3000 yr, nuclear lines are clearly seen in the gamma-ray spectra of NSM remnants due to the minimal Doppler-broadening effect, as is shown in Figures 4 and 5. Using the gamma-ray spectra of NSM remnants that were shown in Section 2 (i.e., the Y e distribution for the NSM case with M ej = 0.01 M e at d = 10 kpc), we selected the brightest nuclear lines in each energy band, 3-75, 75-500, and 500-4000 keV, for at least one epoch in the age range spanning t = 10-4 × 10 6 yr. Note that these energy bands are defined such that they simulate the energy bands that are observable by current and near-future  Figure 9 presents the time evolution of the brightest nuclear gamma-ray lines in these energy bands. To account for the reduction in the line sensitivities as a result of the Dopplerbroadening effect, we accumulated the photons that were within the energy resolution of ΔE = 3 ± 1 keV from the center of energy of their associated lines. This chosen value for the energy resolution is typical for semiconductor gamma-ray detectors. For reference, the evolution of lines without Doppler broadening is also shown in the figure as dashed lines. As indicated in Figure 9, the Doppler-broadening effect becomes less dominant in the hard X-ray band after a few hundred years, but it is still present until about t = 10 3 and 10 4 yr in the soft gamma-ray and the hard gamma-ray bands, respectively. Note that the reason why several lines, such as those of 126m Sb and 239 Np, increase as t approaches 10 3 -10 5 yr is that the number of parent nuclei increases in these phases.
From Figure 9, we can identify the nuclear lines that are useful as indicators for the ages of NSMs. The ages can be categorized into three epochs: t < 100, t ∼ 10 3 -10 4 , and t > 10 4 yr. In summary, if we detect the lines from 125 Sb,194 Os,227 Th,or 194 Ir, then we can determine the age of the NSM to be very young at t < 100 yr. Similarly, lines from 137m Ba in the gamma-ray band indicate that the age is around t ∼ 10 2 yr. In the age range spanning t ∼ 10 3 -10 4 yr, nuclear lines will be detected from 241 Am, 243 Am, 214 Pb, 239 Np, and/or 214 Bi. A nuclear line from 126m Sb indicates that the NSM is very old at t > 10 4 yr. In the wide age range from t = 400-10 5 yr, the line from 126 Sn stays almost constant at 10 −9 photons s −1 cm −2 for a distance of d = 10 kpc, and thus it can be used as a standard candle for measuring d.

Line Diagnostics for the Electron Fraction
In addition to the spectral colors (Section 3.1), nuclear lines can be used to identify NSM remnants among astronomical objects, especially when the remnants are of an older age. Since the NSMs are thought to have both a more neutron-rich environment and a lower Y e condition than SNe (Lattimer & Schramm 1974;Metzger et al. 2010;Wanajo et al. 2011), a new line-diagnostic method utilizing Y e values will be useful for distinguishing NSMs from SNe. In this subsection, we search for gamma-ray line diagnostics for Y e . We use the gamma-ray spectra calculated under the pure Y e conditions in the Y e = 0.10-0.45 range, whereas in the previous sections we used the mixed Y e condition for NSMs.
To identify the best candidates among the nuclear gammaray lines for the identification of Y e , we first selected the five brightest lines for each age, t = 100, 1000, 10 4 , 10 5 , and 10 6 yr, and for each Y e (= 0.10, 0.15, 0.20, 0.25, 0.30, 0.35, 0.40, and 0.45). Then among these 5 (ranks) × 5 (t) × 7 (Y e ) lines, we Note. The † and ‡ marks indicate the lines used in the Y e diagnostics in Section 3.3 for the young (t < 100 yr) and old (t > 100 yr) cases, respectively. Note that the "m" in 137m Ba indicates that it is meta-stable. Figure 6. The energy spectra with Doppler broadening, the same as the red plots in Figures 4 and 5 but normalized to the 11.5 keV flux at 30, 100, 300, 1000, 3000, 10,000, 30,000, 100,000, 300,000, and 1,000,000 yr, which are shown in black, red, green, blue, light blue, magenta, yellow, orange, yellow green, and olive green, respectively. selected the nuclear lines that appeared in two or more of the conditions for t and Y e . In total, 10 gamma-ray lines are selected and are marked as † and ‡ in Table 1 for t < 100 and t > 100 yr, respectively. Therefore, the lines from 137m Ba (661.7 keV), 85 Kr (513.9 keV), and 125 Sb (427.9 keV) in ages below 100 yr are indicators of low-, middle-, and high-Y e environments, respectively. Here, low, middle, and high are numerically defined as Y e ∼ 0.10-0.20, 0.20-0.35, and 0.35-0.45, respectively. In ages older than t = 100 yr, the lines from 225 Ra (40.0 keV), 243 Am (74.7 keV), 239 Np (106.1 keV), 213 Bi (440.5 keV), and 214 Bi (609.3 keV) are emitted from a low-Y e environment, whereas the lines from 60m Co (58.6 keV) and 126 Sn (87.6 keV) become bright in the middle-and high-Y e environments, respectively. Since the absolute flux of a single line changes with respect to t and d, the ratio between two or more lines should be a good indicator of Y e . Figure 10 summarizes the line intensities and their ratios using the 10 nuclei selected above. For simplicity, the plots for young and old ages, t = 10-10 3 and 10 3 -10 6 yr, respectively, are shown separately. In the young age range (top of Figure 10), the ratios of 137m Ba/ 85 Kr, 125 Sb/ 137m Ba, and 85 Kr/ 125 Sb become larger than unity in the low-, middle-, and high-Y e environments, respectively. These low-and high-Y e indicators (i.e., 137m Ba/ 85 Kr and 85 Kr/ 125 Sb, respectively) exhibit more prominent ratios over time since 85 Kr decays slower than 125 Sb and faster than 137m Ba, whereas the middle indicator ( 125 Sb/ 137m Ba) becomes dim after 100 yr. Note that the plots use the incident line fluxes calculated in step 2 in Section 2.1, and the reduction due to the Doppler-broadening effect is not considered. The Doppler effect is particularly significant in plots for the young age range (top of Figure 10). Quantitatively, the ratios of 137m Ba/ 85 Kr, 125 Sb/ 137m Ba, and 85 Kr/ 125 Sb change by factors of 1.57, 0.82, and 0.77, respectively, for t < 1, 000 yr. In the old age range (middle of Figure 10), the ratios of 243 Am/ 60m Co, 126 Sn/ 243 Am, and 60m Co/ 126 Sn indicate the low-, middle-, and high-Y e environments, respectively. The line from 239 Np has the same flux and time evolution as that from 243 Am (red plots) because they are in the same decay chain. Similarly, the lines from 214 Bi, 225 Ra, and 213 Bi (green plots) follow almost the same trend as those from 243 Am and 239 Np (red plots). Among them, the low-Y e indicator ( 243 Am/ 60m Co) in the old age range is valid up to t = 10 6 yr, and the middle-Y e indicator ( 126 Sn/ 243 Am) shows more significant ratios with t > 10 3 yr. On the other hand, lines for the high-Y e indicator 60m Co/ 126 Sn decay quickly and become unavailable after 10 3 yr; that is, if the ratio 60m Co/ 126 Sn is larger than unity, then the object is in a high-Y e environment with an age of t ∼ 10 3 yr. In summary, using these indicators, which become larger than unity in specific Y e conditions, we can estimate the Y e environment independently from the spectral color diagnostics shown in Section 3.1.
Finally, we checked the line ratios blended by Y e distributions of the NSM cases. The time evolution is plotted in Figure 10 (bottom). The difference in the Y e -fraction models between Wanajo et al. (2014) and Kullmann et al. (2022) does not affect the trend of the NSMs very much. These gamma-ray lines are also expected to be observed from the remnants of core-collapse SNe, which are considered to be less neutron-rich environment at Y e ∼ 0.5 (Andrews et al. 2020) than NSMs. However, our numerical calculation model in this paper has limitations in estimating gamma-ray radiation from the corecollapse SNe because the mass fractions of the r-process nuclei in the ejecta are different between the maximum Y e condition of our calculation (i.e., Y e = 0.45) and the SNe case (Y e ∼ 0.50), and the neutron-rich nuclei in the nominal core-collapse SNe are predominantly generated via the s-process rather than the rprocess. For reference, we plotted the time evolution of the line ratios in the gamma-ray spectra of Y e = 0.45, which should still reproduce well in an environment with almost-equal numbers of neutrons and protons. According to Figure 10 (bottom), we expect the low-Y e indicators (i.e., 137m Ba/ 85 Kr and 243 Am/ 60m Co) in the NSM case to become many orders of magnitude larger than those in the SNe case. In the corecollapse SNe where neutron-rich nuclei are generated via the sprocess, a relatively large amount of 85 Kr and almost no 243 Am are synthesized. Therefore, the difference in these low-Y e indicators between the NSMs and SNe cases is expected to Figure 8. (Top) A color-color diagram of the flux ratio between the 10-25 and 25-70 keV bands versus that between the 25-70 and 70-500 keV bands is shown using simulated gamma-ray spectra for the case of an NSM at a distance of 10 kpc, with the Doppler-broadening effect and an initial velocity of 0.3 c. The spectral evaluations of the NSM both with and without K X-ray emission are shown by magenta and gray lines, respectively, and those by the other Y efraction models, DD2-125145, DD2-135135, SFHo-125145, and SFHo-135135 in Kullmann et al. (2022), are shown in the light-pink lines, for reference. The color-color diagram of the X-ray objects listed in the INTEGRAL catalog (version 0043) are also plotted: X-ray binaries with the purple squares, galactic compact objects, such as cataclysmic variables, with the green squares, SN remnants with the cyan circles, active galactic nuclei with the orange circles, and unidentified objects with the red triangles.
(Bottom) Same plot as the top panel but in the 70-500, 500-1000, and 1000-3000 keV bands. become larger than those shown in Figure 10 (bottom). As for the middle-Y e indicators ( 125 Sb/ 137m Ba and 126 Sn/ 243 Am), they may not be useful in distinguishing gamma-rays from NSMs and SNe according to Figure 10 (bottom). In the sprocess environment of core-collapse SNe, almost no 137m Ba and 243 Am are synthesized and thus these middle-Y e indicators can be larger than the values in the figure. Finally, the high-Y e indicators ( 85 Kr/ 125 Sb and 60m Co/ 126 Sn) can discard the NSMs from the SNe cases as indicated in Figure 10 (bottom). In summary, the new line diagnostic method for Y e provides a tool for distinguishing between NSMs and SNe.

Discussion
In Section 2, we presented a nuclear-decay simulation using a large nuclear database, the goal of which was to estimate the gamma-ray spectra of NSMs up to the age of t = 10 6 yr. We have identified many nuclear lines, listed in Table 1, that can be used for identifying the nucleosynthesis environments of NSMs, even with the Doppler-broadening effect altering the profiles of these lines in the young age range. In Section 3, we numerically analyzed the simulated gamma-ray spectra from NSMs and found that the spectral slope in the soft gamma-ray band above 500 keV changes at around t = 200-300 yr. We also found that the spectral colors of NSMs in the hard X-ray to soft gamma-ray bands differ from those of other astronomical objects up to t = 10 5 yr old. Consequently, we can identify a gamma-ray object as an NSM remnant using the gamma-ray spectral colors (Section 3.1). Among the many nuclear lines in the spectra, we identified that the nuclear lines from 241 Am, 243 Am, 214 Pb, 239 Np, and 214 Bi are prominent for t = 10 3 -10 4 yr, and that the lines from 126 Sn and 126m Sb are prominent for t > 10 4 yr (Section 3.2). In addition, we proposed a new line diagnostic method for distinguishing Y e environments that use the line ratios of 137m Ba/ 85 K and 243 Am/ 60m Co, which become larger than unity for low-Y e objects with young and old ages, respectively (Section 3.3). This diagnostic method distinguishes NSMs from SNe. In the next section, we focus on the sensitivities in the gamma-ray band that are required for current and future megaelectronvolt gamma-ray missions that aim to search for Galactic NSM remnants.

Detectable Distance to Galactic NSM Remnants
A gamma-ray flux from the brightest line in a particular age range can be used to estimate the distance that a virtual gamma-ray instrument with a certain line sensitivity will be able to detect. If we assume the age of the NSM remnants in Figure 9, then we can estimate the limit of the distance that can be detected with the line sensitivity of a specific instrument. The top of Figure 11 summarizes the achievable limit of d for NSM remnants as a function of t and is calculated for the three energy bands 3-75, 75-500, and 500-4,000 keV. For example, instruments with a line sensitivity of 10 −7 photons s −1 cm −2 in the 3-75 keV band (red line at the top of Figure 11), such as Hitomi HXI (Takahashi et al. 2014) and NuSTAR (Harrison et al. 2013), can observe the brightest lines from NSM remnants with t < 10 3 yr at d = 3 kpc. We also checked the degradation of the distance limit due to the Doppler-broadening effect, as shown in the middle of Figure 11, but the results do not dramatically change. If the line sensitivities are the same among the three energy bands, then hard X-rays (thick lines) will be a powerful tool in the search for NSMs that are younger than t < 10 3 yr, but the gamma-ray observations (dotted or dashed lines) are better for surveying NSMs that are older than t > 10 4 yr.
The G4.8+6.2 associated with AD 1163 is one example, from the middle of Figure 11, that provides the requirement for the gamma-ray sensitivity needed to observe an NSM remnant with a known distance and age. The object is reported to be a young kilonova remnant with t ∼ 860 yr (Liu et al. 2019). If it is an NSM remnant at d ∼ 10 kpc, then a sensitivity of 10 −8 and 10 −9 photons s −1 cm −2 is required to observe G4.8+6.2 in the hard X-ray and gamma-ray bands, respectively. This sensitivity is roughly one or two (or more) orders of magnitude deeper than that of INTEGRAL IBIS (Winkler et al. 2003). If the distance is closer at d ∼ 3 kpc, then hard X-ray instruments with a sensitivity of 10 −7 photons s −1 cm −2 in the 3-75 keV band, such as Hitomi HXI (Takahashi et al. 2014) and NuSTAR (Harrison et al. 2013), are expected to be able to observe the emissions from the object.

Direct Estimation of Local NSM Rates Using Gamma-Rays
To estimate the coverage of Galactic NSM remnants observable for specific gamma-ray sensitivities, we first prepare a probability map for the existence of Galactic NSMs. This is given in the same plane as the top and middle of Figure 11 (the td plane). Since NSMs are not uniformly distributed in our Galaxy, we apply the probabilities for NSMs in the d and t spaces given by Wu et al. (2019) and multiply them to get the plot shown in the bottom of Figure 11. We assume that the NSMs are primarily concentrated around the Galactic plane. Most of the NSMs are expected to exist at around d ∼ 8 kpc and t ∼ 10 4 -10 6 yr, as has already been described in Wu et al. (2019).
We then accumulate the probabilities for the existence of NSMs (bottom of Figure 11) within the distance-limit curves (top and middle of Figure 11). As a result, we obtained the coverage of Galactic NSMs as a function of the line sensitivity; this is shown at the top of Figure 12. For example, if we survey Galactic NSM remnants with an instrument having a line sensitivity of 10 −8 photons s −1 cm −2 in the 3-75 keV band, then we expect to observe about 3% of the NSMs in our Galaxy with M ej = 0.01 M e . This value corresponds to about one object if we assume an NSM rate in our Galaxy of 30 per 10 6 yr. In addition, we performed the same procedure to estimate the NSM coverage in units of erg per second per cubic centimeter, which requires a sensitivity that is g E 1 times higher than that required for units of photons per second per cubic centimeter. Here, E γ is the photon energy (the energy of the gamma-ray line). The results are shown in the bottom of Figure 12. Therefore, instruments that can achieve a sensitivity of 10 −14 erg s −1 cm −2 in the 75-500 or the 500-4000 keV bands are expected to be able to observe one NSM remnant  Table 1 are shown in the upper panel. The reduction due to the Doppler-broadening effect is not considered. The line intensities of 137m Ba, 85 Kr, and 125 Sb are shown in red, black, and blue, respectively. Thick, dotted, and dashed lines represent the intensities at t = 10, 100, and 1000 yr. The ratios among the lines are plotted in the lower panel; the ratios between 137m Ba and 85 Kr, 125 Sb and 137m Ba, and 85 Kr and 125 Sb are shown in purple, orange, and cyan, respectively. (Middle) Same plot as the top panel but for older ages (t > 100 yr, ‡ marks in Table 1). The dependencies of 243 Am (and 239 Np), 214 Bi (and 225 Ra, 213 Bi), 60m Co, and 126 Sn are shown in red, green, black, and blue, respectively, and the thick, dotted, dashed, and dotted dash lines represent the intensity at t = 10 3 , 10 4 , 10 5 , and 10 6 yr, respectively, in the upper panel. The ratios in the lower panel between 243 Am and 60m Co, between 214 Bi and 60m Co, between 126 Sn and 243 Am, between 126 Sn and 214 Bi, and between 60m Co and 126 Sn are shown in purple, brown, orange, dark yellow, and cyan, respectively. (Bottom) The time dependencies of the line ratios of 137m Ba/ 85 Kr, 125 Sb/ 137m Ba, 85 Kr/ 125 Sb, 243 Am/ 60m Co, 126 Sn/ 243 Am, and 60m Co/ 126 Sn are shown by the thin purple, thin orange, thin cyan, thick purple, thick orange, and thick cyan lines, respectively. The NSM case by the mass fraction of Wanajo et al. (2014) and the case of Y e = 0.45 are shown by the straight and dotted lines, respectively. The results of the other NSM models, DD2-125145, DD2-135135, SFHo-125145, and SFHo-135135 in Kullmann et al. (2022), are shown in the lighter colors.
with M ej = 0.01 M e in our Galaxy under the same assumption of the NSM rate mentioned above. Similarly, a sensitivity of 10 −15 erg s −1 cm −2 is required in the hard X-ray band to observe one object with the same M ej .
The NSM rates from previous studies are summarized in Figure 13. The NSM rates are estimated using several methods, and even though the values approach each other recently, they still have non-negligible uncertainties or systematic errors that are dependent on the methods used. According to Figure 12, instruments with higher sensitivities can cover more than 10% of NSMs and should be able to observe multiple Galactic NSM remnants (meaning a sensitivity of 10 −9.5 -10 −8.5 photons s −1 cm −2 or 10 −16.5 -10 −14.5 erg s −1 cm −2 in the hard X-ray to gamma-ray bands). The actual numbers observed by future NSM surveys with highly sensitive instruments will provide direct information for the NSM rate in the local universe.

Sensitivity Requirements for Future Missions
To assess the feasibility of detecting Galactic NSM remnants using past, current, and future gamma-ray missions, the gamma-ray spectra expected from NSMs (Figures 4 and 5) Figure 11. (Top, middle) The distance limit of NSMs as a function of t is shown with the line sensitivities given in units of photons per second per cubic centimeter; see keys for details. The thick, dotted, and dashed lines represent the results for the 3-75, 75-500, and 500-4000 keV bands, respectively. The top and middle panels show the results without and with the Dopplerbroadening effect, respectively. The line photons are accumulated within ΔE = 3 keV in the middle panel. Note that the small jump in the 75-500 keV data at t ∼ 10 3 yr in the middle panel is due to the interaction between the two brightest lines when the Doppler-broadening effect is applied. (Bottom) Probabilities for the existence of NSMs, which we took as an assumption in the calculation for Figure 12, in the t-d plane. Figure 12. Coverage of NSM remnants in our Galaxy as a function of the narrow-band sensitivity is shown. In the top panel, the results are given in units of photons per second per square centimeter. The magenta and red plots assume M ej = 0.01 M e with and without the Doppler effect, respectively. Similarly, the cyan and blue plots assume M ej = 0.05 M e with and without the Doppler effect. The thick, dotted, and dashed lines represent the coverage values in the 3-75 keV, 75-500 keV, and 500-4,000 keV bands, respectively. The top axis (shown in green) represents the detectable number of NSMs under the assumption of an NSM rate in our Galaxy of 30 per 10 6 yr. The bottom panel shows the sensitivity results given in units of erg per second per square centimeter, but omitting the results without the Doppler effect. Figure 13. The NSM rate in our Galaxy estimated by previous studies as a function of the published year. The red, magenta, green, and blue plots represent the estimation from the binary or double pulsar population (Belczynski et al. 2002;Voss & Tauris 2003;Kim et al. 2010Kim et al. , 2015Chruslinska et al. 2017;Pol et al. 2019;Grunthal et al. 2021), star evolution (Belczynski et al. 2010;Dominik et al. 2012;Mennekens & Vanbeveren 2014;Artale et al. 2019;Olejak et al. 2020;Chu et al. 2022), short gamma-ray bursts (Petrillo et al. 2013;Jin et al. 2015), and gravitational wave events (Abbott et al. 2019(Abbott et al. , 2020(Abbott et al. , 2021, respectively. The values from the gravitational wave events in Gpc −3 yr −1 are converted into the Myr −1 unit under the assumption of a galactic density of 0.01 galaxy Mpc −3 . Figure 14. The 3σ sensitivities of the missions in the hard X-ray to gamma-ray bands are compared with the gamma-ray expected from NSMs that are 100, 10 4 , and 10 6 yr old, with the assumption of an ejecta mass of 0.05 M · at a 10 kpc distance with ΔE = 3 keV. Past and future missions are shown with the dotted and thick lines, respectively. The sensitivities are taken from the following references: CGRO/COMPTEL (9 yr Harrison et al. (2013), FORCE (3.5σ, 10 6 s) from Nakazawa et al. (2018), Hitomi HXI and SGD (100 ks) from Takahashi et al. (2014), SMILE3 (10 6 s, best condition) from Takada et al. (2022), and GRAMS (1 yr) from Aramaki et al. (2020). For reference, 1.0, 10 −3 , 10 −6 , and 10 −9 times the hard X-ray flux from the Crab Nebula (Kouzu et al. 2013), with a simple extension into the gamma-ray band with the single power-law spectrum, are shown in dashed lines and are noted as m-crab, μ-crab, and n-crab, respectively.