Did a Kilonova Set Off in Our Galactic Backyard 3.5 Myr ago?

The recent detection of the live isotopes 60Fe and 244Pu in deep ocean sediments dating back to the past 3–4 Myr poses a serious challenge to the identification of their production site(s). While 60Fe is usually attributed to standard core-collapse supernovae, actinides are r-process nucleosynthesis yields, which are believed to be synthesized in rare events, such as special classes of supernovae or binary mergers involving at least one neutron star. Previous works concluded that a single binary neutron star merger cannot explain the observed isotopic ratio. In this work, we consider a set of numerical simulations of binary neutron star mergers producing long-lived massive remnants expelling both dynamical and spiral-wave wind ejecta. The latter, due to a stronger neutrino irradiation, also produce iron-group elements. Assuming that large-scale mixing is inefficient before the fading of the kilonova remnant and that the spiral-wave wind is sustained over a 100–200 ms timescale, the ejecta emitted at mid-high latitudes provide a 244Pu over 60Fe ratio compatible with observations. The merger could have happened 80–150 pc away from the Earth and between 3.5 and 4.5 Myr ago. We also compute expected isotopic ratios for eight other live radioactive nuclides showing that the proposed binary neutron star merger scenario is distinguishable from other scenarios proposed in the literature.


INTRODUCTION
The production of half of the elements heavier than iron in the Universe (including all the elements heavier than lead) is due to the rapid neutron capture process, (r-process Burbidge et al. 1957).Despite much progress in the past few years, our precise understanding of the r-process and of its yields is presently limited by uncertainties of both nuclear physical and astrophysical nature (Cowan et al. 2021).The former are mostly due to the paucity of experimental measurement of exotic neutron-rich nuclei, while the latter are related to limitations in the modeling of the astrophysical sites.
Even the sites where r-process nucleosynthesis happens are still uncertain.The detection of the kilonova AT2017gfo unambiguously associated to binary neutron star (BNS) merger GW170817 provided the first direct evidence of r-process nucleosynthesis (Pian et al. 2017;Smartt et al. 2017;Kasen et al. 2017).The question whether compact binary mergers involving at least one neutron star (NS) are the only relevant site is still debated.Compact binary mergers seem to have issues in accounting for all available observables, including the abundances of r-process elements in very metal-poor stars or in ultrafaint dwarf galaxies (Côté et al. 2019;Bonetti et al. 2019).Other possible r-process sites include special types of (rare) supernovas (SNs), such as magnetorotational supernovae (Winteler et al. 2012;Mösta et al. 2018) or collapsars (Siegel et al. 2019), although their viability is debated and they could be limited to low-metallicity environments (e.g.Macias & Ramirez-Ruiz 2019;Bartos & Marka 2019).
The observation of r-process abundance patterns traceable to single events has the potential to shed light on their production site.The detection of live (i.e.undecayed) radioactive isotopes in sediments is powerful in this respect, since it features a nontrivial temporal dependence from their decay profiles (Ellis et al. 1996;Wehmeyer et al. 2023).There is a common consensus that the 60 Fe (mean lifetime: τ60 Fe ≃ 3.8 Myr) observed in terrestrial and lunar samples points to one or more explosive events happening ≲ 10 Myr ago (Mya) not far from the Earth (≲ 120 pc).The production site is in general identified with a SN, while BNS mergers were firmly disregarded (Fields et al. 2005;Fry et al. 2015;Schulreich et al. 2017).A concomitant finding in deepsea sediments of 53 Mn (τ53 Mn ≃ 5.4 Myr), an isotope usually associated with Type Ia SN events, has also been reported at a similar depth of 60 Fe (Korschinek et al. 2020).
Recently, Wallner et al. (2016); Wallner et al. (2021) reported new measurements of 60 Fe in deep-ocean sediments and ferromanganese crusts.The deduced interstellar influx shows two peaks within the last 10 Myr, the most prominent one starting ∼3.5 Mya and centered at ∼2.5 Mya, the smaller and narrower second one peaking at ∼6.5 Mya.Interestingly, they also documented the unambiguous emergence of a 244 Pu (τ244 Pu ≃ 116.3 Myr) signature, especially in association with the most recent and prominent 60 Fe peak.Based on their measurements, they reported an interstellar medium (ISM) 244 Pu fluence at Earth orbit of F244 Pu = (7.7 ± 1.6) × 10 3 atoms cm −2 and an abundance ratio of Y244 Pu /Y60 Fe = (53 ± 6) × 10 −6 for the 0-4.6 Mya time window.Wallner et al. (2021) attributed the two 60 Fe peaks to multiple nearby SNs happening within the last 10 Myr at 50-100 pc from Earth.To explain the 244 Pu abundance, they considered a pure SN origin within the Local Bubble or a combination of SNs with a previous nucleosynthesis event (e.g. a BNS merger).Wang et al. (2021) compared results from Wallner et al. (2021) with yield predictions obtained from SN and BNS merger models.They claimed that the observations are compatible with a single source located at ≲ 100 pc only by considering a SN with an enhanced rprocess production, while the scenario of a single nearby BNS merger would be unfeasible.They also proposed a two-step scenario, in which 244 Pu was produced by a rare and more distant event polluting the Local Bubble, before reaching the Earth among the debris from a nearby SN.We notice that Wang et al. (2021) considered isotropic ejecta.Moreover, the abundance yields of each source were obtained as a linear combination of a few representative trajectories obtained in simulations and fitted to reproduce metal-poor star observations.
In this Letter, we revive the single BNS merger scenario using exclusively the outcome of BNS merger simulations and the data reported by Wallner et al. (2021).We consider the case in which the merger produces a long-lived remnant.By taking into account the prolonged effect of neutrino irradiation and the resulting anisotropy in the nucleosynthesis yields, we find that the coincident detection of 60 Fe and 244 Pu in the more recent portion of the crust (≲ 4 Mya) is compatible with a BNS merger event occurring at a distance of ∼ 80 − 150 pc.We consider the outcome of six BNS merger simulations originally presented in Bernuzzi et al. (2020); Nedora et al. (2021) and part of the CORE database (Gonzalez et al. 2023).They were targeted to GW170817 (i.e.their chirp mass is 1.188 M ⊙ , see Abbott et al. 2019), that we consider as representative BNS merger event and for which we explore different mass ratios, q ∈ [0.7, 1].A summary of the models is reported in Table 1.Matter was evolved employing the WhiskyTHC code (Radice & Rezzolla 2012;Radice et al. 2014a,b), complemented by a finite temperature, compositiondependent equation of state (EOS).The EOSs used in those simulations were BLh (Logoteta et al. 2021), HS(DD2) (Hempel & Schaffner-Bielich 2010;Typel et al. 2010, hereafter DD2), SFHo (Steiner et al. 2013) and SRO(SLy4) (Douchin & Haensel 2001;Schneider et al. 2017, hereafter SLy4).Neutrino radiation was taken into account by a leakage scheme to model neutrino emission and a M0 transport scheme to account for absorption in optically thin conditions (Radice et al. 2016(Radice et al. , 2018)).These schemes were shown to describe well the most relevant features of neutrino emission and reabsorption (Zappa et al. 2023).The latter effect is crucial to predict the properties of the unbound matter, ultimately influencing the ejecta composition (e.g.Sekiguchi et al. 2015;Foucart et al. 2016;Perego et al. 2017;Radice et al. 2018).In all but one simulation, turbulent viscosity of magnetic origin was included via a large eddy scheme (Radice 2017(Radice , 2020)).Each simulation covered the innermost part of the domain with a grid of resolution of ∆x = 185 m.

METHODS
In addition to unbind dynamical ejecta, these simulations produced a long-lived merger remnant, lasting t end ∼ 40 − 110 ms postmerger and showing no sign of gravitational collapse up to these times.Such a merger outcome produces spiral-wave winds (Nedora et al. 2019(Nedora et al. , 2021) ) that possibly unbind an amount of matter significantly larger than the dynamical one.Moreover, the longer exposition to neutrino irradiation increases the electron fraction (Y e ) in these ejecta.
The unbound matter properties are extracted from the simulations on a sphere of coordinate radius R E = 200 M ⊙ ≃ 295 km.We use the geodesic and the Bernoulli criterion to identify the dynamical and the spiral-wave ejecta, respectively.We extract mass histograms in the space characterized by the specific entropy (s), Y e and expansion timescale (τ exp ).The latter is obtained from the radial speed and density at R E according to the method outlined in Radice et al. (2016Radice et al. ( , 2018)).We keep track of the spatial composition of the ejecta along the polar angle θ (measured with respect to the binary's orbital axis), while we marginalize over the azimuthal angle.
We compute the nucleosynthesis yields produced in each simulation convolving the ejecta properties with the outcome of nuclear network calculations performed with SkyNet (Lippuner & Roberts 2017).More details on how the isotopic masses are extracted for each ejecta component are given in Appendix A. Overall, we compute the mass of each isotope i at different polar angles as the sum of the dynamical (m dyn ej,i ) and spiral-wave wind (m wind ej,i ) ejecta contribution at 30 years after merger.Since the latter did not saturate at the end of the simulations, but had an approximately constant ejection rate, we rescale the spiral-wave wind ejecta yields by a factor f wind (t wind ) ≡ (t wind − t wstart )/(t end − t wstart ), such that m ej,i (θ, t wind ) = m dyn ej,i (θ) + f wind (t wind )m wind ej,i (θ), where t wstart = 20 ms is the spiral-wave wind onset (Nedora et al. 2019(Nedora et al. , 2021) ) and t wind ∈ [50,200 ]ms.The advantage of using such a rescaling procedure is twofold: i) to align the outcome of simulations with different durations, ii) to explore the effect of longer winds, whose duration is still comparable to the simulated one.
For a BNS merger for which θ is the polar angle pointing toward Earth, and that happened at a time t ≫ 30 yr in the past, the abundance ratio for isotopes i and j of mass numbers A i and A j is The fluence of isotope i measured on the Earth, F i , is related to the mass of the ejecta and to its radioactivity distance D rad,i by where f dust,i is the fraction of i isotopes that form dust, m iso ej,i the isotropized ejecta mass of the isotope i emitted in the direction θ and m u the atomic mass unit.For consistency with Wallner et al. (2021), we set f dust, 244 Pu = f dust, 60 Fe = 0.5.For 244 Pu, we use the fluence value reported in Table 2 of Wallner et al. (2021), while the fluence of 60 Fe is calculated assuming the same fluence over layer incorporation ratio of 244 Pu (also taken from Table 2 of Wallner et al. (2021)).
To gauge the viability of our BNS scenario, the radioactivity distance of different isotopes must be mutually compatible and has to be compared with some relevant length scales.The first is the fading radius, R fade .Upon expulsion, the ejecta expand homologously.Then they enter first a self-similar Sedov-Taylor expansion phase, and then a snow-plow phase.Finally they dissolve into the ISM upon reaching R fade (Montes et al. 2016).Using the model outlined in Beniamini et al. (2018); Bonetti et al. (2019), assuming fiducial values for the ejecta mass and speed of 0.04 M ⊙ and 0.2 c, respectively, and a neutral hydrogen density of 0.05 atoms cm −3 for the Local Bubble (Zucker et al. 2022), we estimate R fade ≃ 240 pc.The second length scale is the radius of the Local Bubble, R LB ≃ 165±6 pc (Zucker et al. 2022).We also estimate the typical timescale the BNS ejecta would need to expand to such radii.Considering the kilonova remnant models presented above and the outcome of kilonova remnant simulations (Montes et al. 2016), the remnant radius could reach ∼100 pc within ∼1 Myr.Finally, we also assume that until fading into the ISM the ejecta do not undergo large-scale mixing, so the angular dependence of the ejecta is relevant and they enrich the surrounding space in an anisotropic way.Under this assumption, the isotopic ratios emerging from the BNS event can be directly reflected into those at the Earth's orbit once accounting for the decay, since 60 Fe and 244 Pu (and all the other live radioactive isotopes) expand across the Local Bubble with the same history.

RESULTS
In Fig.
(1), we present Y 60 Fe /Y 244 Pu for a merger happening t = 3.5 Mya.For the q = 1 and q = 0.7 BLh models there exists a relatively broad range of polar angles at mid-high latitude (20 • ≲ θ ≲ 60 • ) for which the ratio matches the observed value.The angular intervals are 40 • ≲ θ 1 ≲ 60 • and 20 • ≲ θ 2 ≲ 40 • for q = 1 and q = 0.7, respectively.For matter expelled at those latitudes, 244 Pu is mostly produced in the dynamical ejecta, while the 60 Fe comes from the spiral-wave wind (see also Fig. ( 4)).In particular, for 20 • ≲ θ ≲ 70 • the distribution of 60 Fe in the spiral-wave wind ejecta is rather flat and not very sensitive to the mass ratio.Instead, the 244 Pu distribution decreases moving from the equator to the poles but with a shallower profile for the unequal mass case.Models other than the BLh ones fail to reproduce the observed ratio.We notice however BLh (q = 0.7) DD2 (q = 1.0)DD2* (q = 0.83) SFHo (q = 0.7) SLy4 (q = 0.7) Figure 1.Isotopic ratio between 60 Fe and 244 Pu for all the models listed in Table 1 as a function of the polar angle θ.The bands represent the variability in the wind duration, t wind ∈ [50, 200] ms, with larger (smaller) ratios related to a longer (shorter) duration.The horizontal band corresponds to the measured ratio (Wallner et al. 2021).that all the simulations share a similar behavior and some of them still produce a ratio which is not too far from the observed one.In the case of the DD2 models, the spiral-wave wind is not rich enough in 60 Fe with respect to the plutonium-rich dynamical ejecta.In the case of the SFHo and SLy4 models, the amount of 60 Fe is comparable to the one in the BLh q = 0.7 case, but the amount of 244 Pu is one order of magnitude larger.In the rest of the analysis we will focus on the BLh models that match the measured isotopic ratio.
In Fig.
(2) we present the radioactivity distance for both 60 Fe and 244 Pu.For each BNS model, we consider a variable spiral-wave wind duration and we span the angular intervals θ 1 and θ 2 presented above.For both models, there exist a relatively broad t wind interval in which the radioactivity distances of the two isotopes are mutually compatible.Depending on the model, D rad is compatible with an explosion happening between ∼80 and 150 pc from Earth.These values are roughly comparable with the Local Bubble radius.Crucially, they are sufficiently distant from the Earth to avoid life extinction (D rad ≳ 10 pc, see e.g.Perkins et al. 2024) but not too distant for the kilonova remnant to dissolve before reaching the Earth (i.e.D rad ≲ R fade ).
Given the possible ∼ 1 Myr delay between the merger and the arrival of the ejecta on the Earth, we test the robustness of our results with respect to the explosion time.In BLh (q = 0.7) Figure 2. Radioactivity distances as a function of t wind for the BLh models with q = 1.0 (left) and q = 0.7 (right).The BNS merger is assumed to occur t = 3.5 Mya.The bands account for different polar angles, i.e. θ1 ∈ [40 • , 60 • ] for q = 1, θ2 ∈ [20 • , 40 • ] for q = 0.7.The horizontal lines mark R LB and R fade .

Pu)
BLh (q = 1.0) -twind = 200 ms BLh (q = 0.7) -twind = 100 ms Figure 3. Isotopic ratio between 60 Fe and 244 Pu as a function of time for the BLh, q = 1.0 model with t wind = 200 ms (blue) and for the BLh, q = 0.7 model with t wind = 100 ms (orange).The bands represent the variability in the polar angle considering the same intervals as in Fig. (2).The horizontal and vertical bands show instead the measured ratio (Wallner et al. 2021) and a ±1 Myr uncertainty about the fiducial value t = 3.5 Mya, respectively.
patible radioactivity distances can also accommodate a relatively large explosion time uncertainty (±1 Myr).

DISCUSSION
Note-The intervals span uncertainties in the polar angle and t wind .The BNS merger is assumed to occur t = 3.5 Mya.
Our results show that the coincident excess of 60 Fe and 244 Pu observed in deep ocean sediments, dating back to ∼ 3 − 4 Mya, can be explained as the result of a single BNS merger that happened 80-150 pc away from our Solar System.
The difference between our results and the ones presented in several previous papers (e.g.Fry et al. 2015;Wang et al. 2021Wang et al. , 2023) ) can be understood in terms of some specific features that characterize our BNS models.First of all, the merger remnant must consist in a massive NS, not collapsing to a black hole over a time scale of 100-200 ms, to produce a significant spiral-wave wind ejecta in addition to the dynamical ejecta.The presence of both these two components is essential, since 60 Fe is synthesized in the former, while a significant amount of 244 Pu in the latter.Depending on the EOS stiffness and on the colliding NS masses, such an outcome is expected to be relatively frequent and possibly larger than 50% of the cases (Margalit & Metzger 2019), especially if the recent detections of massive NSs will be confirmed (Fonseca et al. 2021;Riley et al. 2021;Romani et al. 2022).Moreover, in our analysis we retain information about the angular distribution of the ejecta and we find that the precise conditions to match the observations are realized only in matter expelled at mid-high latitudes, i.e. for a viewing angle 30 • ≲ θ ≲ 50 • , with an angular width of ∆θ ≈ 20 • .The corresponding solid angle fraction is ∆Ω/4π = 2 sin θ sin (∆θ/2) ≈ 0.35 sin θ, which ranges between 0.18 and 0.27.Despite not being realized in the majority of cases, the probability of observing a BNS merger in those conditions is not negligible and not even small.
Since our models disfavor viewing angles very close to the poles, the relativistic jet that could have originated from such an event would not have hit the Earth due to its small opening angle (θ jet ≲ 6 • , see e.g.Fong et al. 2015;Perkins et al. 2024).Moreover, the presence of a relatively broad range of θ still allows the possibility for the ejecta to mix, at least over an angular scale ≲ ∆θ, due to the lateral expansion that becomes relevant once the strong blast wave has converted much of its kinetic energy into thermal energy (Montes et al. 2016).Large scale, turbulent mixing is expected to occur only once the expansion speed has reached the ISM sound speed (∼ 10 km s −1 ) and the kilonova remnant starts to fade away, on timescales of a few hundreds Myr (Hotokezaka et al. 2015;Beniamini & Hotokezaka 2020;Kolborg et al. 2023).
It must be noticed how BNS mergers (and, more in general, events in which r -process nucleosynthesis occurs) are expected to be rare (Hotokezaka et al. 2015;Abbott et al. 2023), making their nearby occurrence in the recent past an exceptional event.Additionally, our analysis does not rule out the single supernova origin or the two-step model discussed in previous works, e.g.Wallner et al. (2021); Wang et al. (2021Wang et al. ( , 2023)); Wehmeyer et al. (2023).The identification of other relevant isotopic ratios could be the key to discriminate between the different scenarios, as suggested in Wang et al. (2021Wang et al. ( , 2023)).To this end, in Table 2 we provide the isotopic ratios with respect to 244 Pu of eight other live radioactive isotopes for the two representative BNS models discussed in Fig.
( 2) and ( 3).Within each model, for 93 Zr, 107 Pd and 129 I, the values of the isotopic ratios are proportional to t wind , while for 135 Cs, 182 Hf and all the actinides the dependence on t wind is weak or even negligible.Among the different ratios, the largest values are observed for 107 Pd, followed by 93 Zr and 129 I, lower by one order of magnitude.These trends are different compared with the values presented in Wang et al. (2021Wang et al. ( , 2023)), for which the largest ratio is always realized for 129 I.For these isotopes, as well as for 135 Cs and 182 Hf, the values extracted from our models are intermediate between the larger values obtained by magnetically driven SN and the smaller ones obtained by considering enhanced r -process SN wind models.The ratios extracted for the actinides are similar to the ones reported by Wang et al. (2021Wang et al. ( , 2023)), confirming that they have a low discriminating power.We also look at the production of 53 Mn, that we find to occur only for very specific thermodynamics conditions (Y e ≳ 0.45).Our BNS models are not able to reproduce the 53 Mn over 60 Fe ratio of 2:1 in the interstellar dust predicted by Korschinek et al. (2020), due to the very small amount of ejecta with such high electron fraction (≲ 10 −6 M ⊙ ).However, more recent BNS merger simulations employing more detailed neutrino transport (Zappa et al. 2023;Espino et al. 2023) suggest the presence of a significant amount of ejecta in the high Y e tail, especially in the case of long-lived remnants, matching the conditions required to produce 53 Mn.
In our analysis we focused on the coincident 60 Fe and 244 Pu peaks observed in the youngest deep oceanic crust (3-4 Mya), for which the amount and quality of data are more significant.Due to the paucity of expected nearby BNS mergers, it is implausible that a similar event can also explain the previous, smaller 60 Fe peak, especially if associated with a nonnegligible amount of 244 Pu.However, several alternative solutions were previously discussed for the interpretation of the older peak, which apply also to our scenario, including the fact that the older peak could originate from outside the Local Bubble (e.g.Schulreich et al. 2023).
We stress that the single BNS event scenario can explain the observational data without any need for tuning.Indeed the only free parameter in our model is the spiral-wave wind duration, and we conservatively vary it over a range comparable to the duration of our simulations.To better address the viability of the kilonova scenario, more realistic models would be necessary.But, if confirmed, our analysis shows that the merger of a BNS system could have happened in the Solar neighborhood no earlier than ∼ 4 Mya.

A. NUCLEOSYNTHESIS CALCULATION IN BNS EJECTA
In order to compute the total isotopic masses expelled in each BNS simulation, we exploit the outcome of an extensive set of nuclear network calculations performed on a grid of 11700 Lagrangian tracer particles.The grid spans broad ranges in the space of initial thermodynamic conditions parametrized by (s, Y e , τ exp ), and is identical to the one used in Perego et al. (2022).We employ the SkyNet nuclear network (Lippuner & Roberts 2017) to evolve in time the number abundances of a wide set of nuclear species, including 60 Fe and 244 Pu, depending on the specific Lagrangian particle's history.The initial nuclear distribution follows from nuclear statistical equilibrium (NSE) conditions that are determined by the (s,Y e ) values of the trajectory, at a temperature fixed to 8 GK (high enough to guarantee the validity of the NSE assumption).The composition is then evolved along the analytic density profile used in Lippuner & Roberts (2015), parametrized by τ exp .The network is run until the final time of ∼ 31.7 years, using the same input nuclear physics as in Perego et al. (2022).To obtain the polar angle dependent yields produced in the merger event we compute the convolution of the mass histograms of the ejecta (as described in the main text and directly extracted from the simulations) with the abundances obtained on the (s, Y e , τ exp ) space by the nuclear network evolution.We proceed separately for the dynamical and spiral-wave wind ejecta components.The total isotopic masses in the ejecta are then obtained by summing the two contributions, after having rescaled the spiral-wave wind one for the corresponding t wind as discussed in Sec.(2).3).In the proposed BNS scenario, 60 Fe is synthesized over a wide angular range and mainly in the spiral-wave wind ejecta.These ejecta are characterized by a relatively high Y e , an effect of the prolonged neutrino irradiation from the central long-lived remnant.A similar behav-ior is observed for 93 Zr and 107 Pd, listed in Table 2.The production of 244 Pu instead, representative of the heavier isotopes in the table, peaks in the equatorial region (θ ≳ 60 • ) and originates from the fraction of cold, neutron-rich matter present in the two ejecta components.
Fig. (3) we present the temporal evolution of Y 60 Fe /Y 244 Pu .The results show that our BNS models that can reproduce the observed isotopic ratio and com-

Figure 4 .
Figure 4. 60 Fe and 244 Pu masses ejected in BLh models with q = 1.0 (left) and q = 0.7 (right) as a function of the polar angle θ.The contribution from the dynamical and spiral-wave wind ejecta are shown separately.The band represents the variability in the wind duration, t wind ∈ [50, 200] ms, with larger (smaller) isotopic masses related to a longer (shorter) duration.
Fig. (4) shows the amount of 60 Fe and 244 Pu produced at 30 years after merger in the two representative BNS models of Fig. (2) and Fig. (

Table 1 .
Summary of the Properties of the BNS Merger Models

Table 2 .
Selected Radioisotope Ratios with Respect to 244 Pu for the BLh q = 1.0 and q = 0.7 BNS Models The authors thank B. Wehmeyer for helpful discussions, the Computational Relativity (CORE) collaboration for providing simulation results, and the European Centre for Theoretical Studies in Nuclear Physics and Related Areas (ECT * ) for hosting the MICRA2023 workshop, where this work was conceived.They also thank the INFN for the usage of computing and storage resources through the tullio cluster in Turin.This work has been supported by STRONG-2020 "The strong interaction at the frontier of knowledge: fundamental research and applications" which received funding from the European Union's Horizon 2020 research and innovation programme under grant agreement No 824093.The work of AP is partially funded by the European Union under NextGenerationEU.PRIN 2022 Prot.n. 2022KX2Z3B.