# Neutron Star Mergers as the Main Source of R-process: Natal Kicks And Inside-Out Evolution to The Rescue

Projjwal Banerjee, Meng-Ru Wu, Zhen Yuan

► To cite this version:

Projjwal Banerjee, Meng-Ru Wu, Zhen Yuan. Neutron Star Mergers as the Main Source of R-process: Natal Kicks And Inside-Out Evolution to The Rescue. The Astrophysical Journal Letters, 2020, 902 (2), pp.L34. 10.3847/2041-8213/abbc0d hal-02981272

HAL Id: hal-02981272

https://hal.science/hal-02981272v1

Submitted on 29 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Neutron Star Mergers as the Main Source of *r*-process: Natal Kicks and Inside-out Evolution to the Rescue

Projjwal Banerjee<fnref n="1" />,*, Meng-Ru Wu<fnref n="2" />,<fnref n="3" />,<fnref n="4" />,†, and Zhen Yuan<fnref n="5" />,<fnref n="6" />,||

## Abstract

Binary neutron star mergers (BNSMs) are currently the most promising source of *r*-process thanks to the detection of GW170817. The estimated occurring frequency and the amount of mass ejected per merger indicate that BNSMs by themselves can account for all the *r*-process enrichment in the Galaxy. However, the decreasing trend of [Eu/Fe] versus [Fe/H] of disk stars for [Fe/H] $\gtrsim -1$ in the solar neighborhood is inconsistent with the flat trend expected from BNSMs with a standard delay time distribution (DTD) $\propto t^{-1}$. This has led to the suggestion that either additional sources or modification to the DTD of BNSMs is required to match the observations. We investigate the effects of natal kicks received during the birth of neutron star binaries on the chemical evolution of *r*-process element Eu in the Milky Way by combining the results from the galactic dynamics code GALPY with a one-zone Galactic chemical evolution model OMEGA. We show that when key inputs from simulations of the inside-out disk evolution are combined with natal kicks, BNSMs can naturally reproduce the observed decreasing trend of [Eu/Fe] with [Fe/H] in the solar neighborhood without the need for modification to the DTD or additional *r*-process sources.

*Unified Astronomy Thesaurus concepts: R-process (1324); Galaxy chemical evolution (580); Solar neighborhood (1509); Milky Way Galaxy physics (1056); Compact objects (288); Neutron stars (1108)*

## 1. Introduction

Following the seminal discovery of GW170817 (Abbott et al. 2017a, 2017b), binary neutron star mergers (BNSMs) have become the first, and currently the only, confirmed site for the synthesis of heavy elements by the rapid neutron capture process (*r*-process; see Cowan et al. 2019; Metzger 2019 for recent reviews). The current estimated BNSM rate of 110–3840 Gpc⁻<fnref n="3" /> yr⁻<fnref n="1" /> (Abbott et al. 2019) along with the amount of *r*-process material ejected per merger of ≈0.03–0.06 M☉ (e.g., Drout et al. 2017; Cowperthwaite et al. 2017; Kasen et al. 2017; Tanaka et al. 2017; Villar et al. 2017; Kawaguchi et al. 2018; Wanajo 2018; Wu et al. 2019) is sufficient to explain all of the *r*-process enrichment in the Galaxy (Côté et al. 2018; Hotokezaka et al. 2018; Cowan et al. 2019). However, if BNSMs are assumed to be the only source of *r*-process, then the Galactic evolution of elements such as Eu, which are primarily produced by *r*-process, is very different from elements that are produced by other sources associated with massive stars. This is due to the inherent delay between the birth of a neutron star binary (NSB) and the eventual merger with a typical delay time distribution (DTD) $\propto t^{-1}$, compared to massive stars that have negligible delay. In particular, with a DTD $\propto t^{-1}$ that is expected from population synthesis calculations (Dominik et al. 2012; Chruslinska et al. 2018; Côté et al. 2019), [Eu/Fe] remains almost constant for stars with [Fe/H] $\gtrsim -1$ in galactic chemical evolution (GCE) calculations due to the fact that both BNSMs and Type Ia supernovae (SNe Ia) have the same DTD (van de Voort et al. 2015; Komiya & Shigeyama 2016).<fnref n="7" /> In sharp contrast, the observed values of [Eu/Fe] in disk stars at the solar neighborhood show a

clearly decreasing trend with increasing [Fe/H] (Battistini & Bensby 2016). This was pointed out to be an issue by Côté et al. (2017a) and Hotokezaka et al. (2018) and has been studied in detail recently by Côté et al. (2019) (see also Schönrich & Weinberg 2019). Possible solutions to the problem include a steeper DTD ($\propto t^{-1.5}$) or a burst of merger at early times followed by a typical DTD $\propto t^{-1}$ (see also Hotokezaka et al. 2018). These, however, have been pointed out to be inconsistent with observations of short gamma-ray bursts and SNe Ia in early-type galaxies (see Côté et al. 2019 for details). Alternatively, additional source(s) of *r*-process can explain the observed trend provided that their frequency decreases with metallicity. This source could be an *r*-process site associated with the death of massive stars such as magnetorotational SNe (Winteler et al. 2012; Mösta et al. 2018), accretion disk outflow from collapsars (Siegel et al. 2019), or supernova explosions associated with hadron-quark phase transition (Fischer et al. 2020). Higher frequency of mergers resulting from neutron star–black hole binaries in the early Galaxy (lower metallicities; Mennekens & Vanbeveren 2014) and *r*-process occurring in accretion disk outflows during the common envelope phase of NS–massive star system (Grichener & Soker 2019) are some of the other possibilities. Such sources may also be needed to explain the presence of *r*-process elements in the very early Galaxy and the large scatter in their abundances as observed in very metal-poor stars with [Fe/H] $\lesssim -2.5$ (Argast et al. 2004; Wehmeyer et al. 2019), but this is still under debate (Tsujimoto & Shigeyama 2014; Hirai et al. 2015; Ishimaru et al. 2015; Shen et al. 2015; van de Voort et al. 2015; Safarzadeh et al. 2019).

An interesting feature that distinguishes BNSMs from other rare *r*-process sites is that they receive large natal kicks during the formation of binaries (Fong & Berger 2013;

<footnote n="1">¹Discipline of Physics, Indian Institute of Technology Palakkad, Kerala, India 678557; [projjwal.banerjee@gmail.com](mailto:projjwal.banerjee@gmail.com)</footnote>

<footnote n="2">²Institute of Physics, Academia Sinica, Taipei, 11529, Taiwan; mwu@gate.sinica.edu.tw</footnote>

<footnote n="3">³Institute of Astronomy and Astrophysics, Academia Sinica, Taipei, 10617, Taiwan</footnote>

<footnote n="4">⁴Physics Division, National Center for Theoretical Sciences, Hsinchu, 30013, Taiwan</footnote>

<footnote n="5">⁵Key Laboratory for Research in Galaxies and Cosmology, Shanghai Astronomical Observatory, Chinese Academy of Sciences, 80 Nandan Road, Shanghai 200030, People's Republic of China; sala.yuan@gmail.com</footnote>

<footnote n="6">⁶Université de Strasbourg, CNRS, Observatoire astronomique de Strasbourg, UMR 7550, F-67000 Strasbourg, France
*Received 2020 July 8; revised 2020 September 11; accepted 2020 September 27; published 2020 October 19*</footnote>

<footnote n="7">⁷We note that the decreasing trend of [Eu/Fe] at high metallicity can be reproduced if a fixed delay time for BNSM is adopted; see, e.g., Argast et al. (2004), Matteucci et al. (2014), and Wehmeyer et al. (2015).</footnote>
Behroozi et al. 2014; see, however, Beniamini & Piran 2016; Tauris et al. 2017). Consequently, the final location of mergers can be very different from their birth location where a certain fraction of NSBs are effectively lost as they do not contribute to the Galactic enrichment of heavy elements due to the fact that they merge far away from the star-forming regions (Safarzadeh & Côté 2017; Safarzadeh & Scannapieco 2017). In this Letter, we show that natal kicks have a large impact even for NSBs that do contribute to the Galactic enrichment of heavy elements. Specifically, we model two new effects due to natal kicks on GCE: the effect of kick-induced migration on the effective BNSM frequency as well as the impact on the effective DTD. We show that when these effects are combined with the inside-out formation of the Milky Way (MW; see e.g., Minchev et al. 2013; Schönrich & McMillan 2017; Frankel et al. 2019), the decreasing trend of [Eu/Fe] versus [Fe/H] for stars with [Fe/H] $\gtrsim -1$ matching the observation can be naturally obtained with BNSMs as the only r-process source with a standard DTD $\propto t^{-1}$.

## 2. Effects of Natal Kicks

It is known, that due to natal kicks, only a fraction $f_{ret} \le 1$ of the total NSBs born at a certain time interval contribute to the Galactic enrichment of heavy elements, whereas the rest merge far from the star-forming regions (Safarzadeh & Côté 2017). In addition, however, there are two other important effects caused by natal kicks that have not been taken into account previously but turn out to be crucial in modeling the GCE of r-process elements in the solar neighborhood. The first one is the kick-induced migration of NSBs within the Galactic disk. Considering a simple model of the MW disk as consisting of independent concentric rings, one can study the GCE in the vicinity of a particular radius that is described by a ring centered at the given radius. For any ring, the GCE of non-r-process elements that receive negligible contribution from BNSMs depends mostly on the local star formation rate (SFR). In contrast, for r-process elements, their GCE in a given ring depends not only on the BNSMs that are both born and merge inside the ring, but also on BNSMs that are born in other rings that migrate and eventually merge within the considered ring. In this regard, we define a useful quantity relevant for GCE calculations, $\eta(R, t)$, as

$$ \eta(R, t) = \frac{N_{\text{merge}}(R, t)}{N_{\text{born}}(R, t)}, \qquad (1) $$

where $N_{\text{born}}(R, t)$ is the number of NSBs born inside the ring centered at radius $R$ between time $t$ and $t + \Delta t$, and $N_{\text{merge}}(R, t)$ is the actual number of NSBs born in the entire disk during the same time interval but eventually merge within the ring centered at $R$. If migration due to natal kicks is neglected, only the NSB that is born inside a ring can merge within that ring such that $\eta$ is 1 for all rings. It is important to note that although $\eta(R, t)$ effectively alters the birth rate of BNSMs at a time $t$, the mergers occur later according to the DTD.

The second important effect of natal kicks is that it impacts the effective DTD of BNSMs for a given ring. This is simply due to the fact that NSBs with shorter merger times tend to be retained and coalesce within the star-forming region of the Galaxy, whereas those with longer merger times have a higher chance of escaping. This leads to lower values of average merger times and

thus results in an effective value of $\beta^{\text{eff}}(R, t) \le \beta$ for an actual DTD $\propto t^{\beta}$.

As mentioned above, the values of $\eta(R, t)$ and $\beta^{\text{eff}}(R, t)$ for a certain vicinity (ring) are influenced by the birth rate of NSBs born both inside and outside the ring. Consequently, in order to calculate their values, it is important to know both the spatial and temporal evolution of the SFR of the entire disk. We use the SFR predicted by a detailed chemodynamical simulation based on the inside-out formation of the MW disk from Minchev et al. (2013). At any given time $t$, we generate the radial coordinate of the birth locations $R_b$ of the NSBs according to a distribution $\propto R \Sigma(R, t)$, where $\Sigma(R, t)$ is the surface SFR density adapted from Minchev et al. (2013), shown in Figure 1(a). Because the starting time, $t_0$, from the simulation in Minchev et al. (2013) is when the bulge is formed, we adopt two different values of $t_0 = 1$ and 2 Gyr. For $t < t_0$, we keep the radial dependence of the SFR the same as that at $t = t_0$. The maximum value of $R_b$ for the birth location of NSBs is limited to 16 kpc corresponding to the maximum value for which SFR is provided in Minchev et al. (2013). For simplicity, we assume that all NSBs are born at a vertical height of $z = 0$.

For each NSB born at a given $R$ and $t$, we assign a kick velocity $\vec{v}_{\text{kick}}$ whose magnitude is randomly generated from an exponential distribution $\propto \exp(-v/v_0)$, with $v_0 = 90 \text{ km s}^{-1}$ similar to Behroozi et al. (2014), which is consistent with the kick velocities inferred from the observed offsets of short gamma-ray bursts by Fong & Berger (2013). As the latter study inferred a nonzero value for the lower limit of $v_{\text{kick}}$, we assume a minimum value of $10 \text{ km s}^{-1}$. The direction of $\vec{v}_{\text{kick}}$ is generated from a uniform and isotropic distribution. The randomly sampled $\vec{v}_{\text{kick}}$ is then added to the velocity of the NSB (just before the birth of the second neutron star), which is assumed to be the circular velocity corresponding to its birth radius.

In order to follow the motions of NSBs under the influence of the Galactic potential until they merge, we use GALPY (Bovy 2015) to trace their orbits. Each NSB is evolved in time until $t_{\text{merge}}$, sampled from DTD $\propto t^{-1}$ with minimum and maximum values of 10 Myr and 10 Gyr, respectively. Because GALPY only allows for a static potential, we restrict values of $t_{\text{merge}} \le 2$ Gyr to minimize the effect of disk growth on the motion of NSB. We use MWPotential2014 in GALPY as the model for the MW potential with the default value of the circular velocity $v_c(R_\odot = 8 \text{ kpc}) = 220 \text{ km s}^{-1}$ at the present time. In order to account for the time evolution of the Galactic potential, we assume that the MW potential is proportional to the virial mass $M_{\text{vir}}$ of the dark matter halo. Thus, for any given $t$, we simply scale the potential by changing the value of $v_c(R_\odot, t)$ given by

$$ v_c(R_\odot, t) = v_c(R_\odot, t_{\text{gal}}) \left( \frac{M_{\text{vir}}(t)}{M_{\text{vir}}(t_{\text{gal}})} \right)^{1/2}, \qquad (2) $$

where $t_{\text{gal}} = 13.7$ Gyr is the age of the Galaxy, and $v_c(R_\odot, t_{\text{gal}}) = 220 \text{ km s}^{-1}$ is the current circular velocity at the solar radius. $M_{\text{vir}}(t)$ is adopted from the average fit reported in Griffen et al. (2016) from simulations of 24 Milky Way-sized halos. We adopt different values for $v_c(R_\odot)$ ranging from $5-220 \text{ km s}^{-1}$ that correspond to ages of $\sim 300$ Myr to the present age of 13.7 Gyr.
<figure><img src="image_3.png" /><figcaption>**Figure 1.** (a) SFR as a function of R from Minchev et al. (2013) at various stages of disk evolution. (b) Evolution of $f_{\text{ret}}$, $\eta_{\odot}$, and $\beta_{\odot}^{\text{eff}}$ for $t_0 = 1$ and 2 Gyr. (c) Distribution of the birth radius of NSBs that merge within the solar neighborhood (shown in gray) at two different stages of the disk for $t_0 = 1$ Gyr. The total number of NSBs at each time is normalized to 100. (d) Same as (b), but using a fixed radial scale length of 3 kpc (no inside-out formation).</figcaption></figure>

<figure><img src="image_4.png" /></figure>

For each $v_c(R_\odot)$, we simulate the motion of (1–4) $\times 10^4$ NSBs according to their birth location, kick velocity, and DTD distribution described above. An NSB is considered to contribute to the Galactic enrichment if it merges within coordinates $R \le R_{\text{max}}$ and $|z| \le z_{\text{max}}$. In order to calculate the fraction $f_{\text{ret}}(t)$ of NSBs retained by the Galaxy at different times, we use fixed values of $R_{\text{max}} = 20$ kpc and $z_{\text{max}} = 5$ kpc. The particular choice of $R_{\text{max}}$ roughly corresponds to the sum of the maximum value of $R_b$ and the typical radial scale length of ~3–4 kpc. The value of $z_{\text{max}}$ was taken to be the sum of the typical disk thickness ~3 kpc found in simulations by Minchev et al. (2013) and the typical remnant radius of a BNSM that explodes a few kiloparsecs above the disk plane. The latter is estimated to be ~2–3 kpc (Thornton et al. 1998) by using a BNSM kinetic energy of ~$10^{51}$ erg and an ambient density of ~$10^{-4}$ cm$^{-3}$ that is expected at a height of a few kiloparsecs from the disk plane (Miller & Bregman 2013). For $\eta(R, t)$ and $\beta^{\text{eff}}(R, t)$, we first focus on the values $\eta_\odot(t)$ and $\beta_\odot^{\text{eff}}(t)$ for the solar ring defined by $7 \le R \le 9$ kpc. The values of $\eta_\odot$ are calculated using Equation (1) where only NSBs that merge within $z_{\text{max}}$ are considered. For $\beta_\odot^{\text{eff}}$, we first compute the average merger time of BNSMs within the solar ring that have $|z| < z_{\text{max}}$. We then use the value of the average merger time to find the corresponding $\beta_\odot^{\text{eff}}$ by assuming a DTD $\propto t^{\beta_\odot^{\text{eff}}}$.

Figure 1(b) shows the evolution of $f_{\text{ret}}$ for the entire Galaxy as well as $\eta_\odot$ and $\beta_\odot^{\text{eff}}$ for the solar ring. At early times ($t \lesssim 1$ Gyr), $f_{\text{ret}}$ ranges from 40% to 50% but increases to ~90% by $t \sim 4$ Gyr. This is similar to the values obtained by Safarzadeh & Côté (2017), who considered a purely dark matter halo with NSBs traveling along the radial direction, and with different criteria for deciding whether a BNSM contributes to the Galactic enrichment. The value of $\beta_\odot^{\text{eff}}$ is always lower than the true value of $\beta = -1$, as expected. The $\beta_\odot^{\text{eff}}$ is more negative at earlier times due to the higher escape rate of BNSMs, and approaches -1 at later times.

The evolution of $\eta_\odot$, on the other hand, is noticeably different from $f_{\text{ret}}$ and $\beta_\odot^{\text{eff}}$. It peaks at early times at $t \simeq t_0$ with a value reaching $\simeq 1.8$ and decreases with time for $t > t_0$. This is a direct consequence of the SFR from Minchev et al. (2013) based on the inside-out formation of the Galaxy. As can be seen from Figure 1(a), the surface SFR is higher at the center but drops sharply with $R$ for the first ~0.5 Gyr after $t_0$, with a typical scale length of $\lesssim 1.5$ kpc. Subsequently, the scale length increases to $\gtrsim 3$ kpc within ~2 Gyr of disk evolution. As a result of the steep drop of SFR with $R$ during the first ~0.5 Gyr of disk evolution, a substantial number of BNSMs that were originally born closer to the center of the Galaxy merge within the solar ring. This can be clearly seen in Figure 1(c), which shows the distribution of the birth radii $R_b$ of BNSMs that
merge within the solar ring at two different stages of the disk evolution. At early stages of the disk evolution ($t_d \lesssim 0.5$ Gyr), $\gtrsim 50\%$ of the BNSMs that merge within the solar ring originated from $R \le 4$ kpc, with the peak of the distribution of $R_b$ at $R = 1–3$ kpc. As the scale length increases with time, at $t_d \sim 4.5$ Gyr, only $\sim 5\%$ of BNSMs that originated from $R \le 4$ kpc contribute to the solar ring, and the peak of the distribution lies within the ring at $R = 7–9$ kpc. Because the radial dependence of the SFR stays constant for $t \le t_0$, $\eta_\odot$ decreases slightly with decreasing values of $t$ as an increasing fraction of NSBs are able to escape due to the shallower Galactic potential.

An important point to note here is that for a given distribution of natal kick velocities, only $\eta_\odot$ is uniquely sensitive to the radial distribution of SFR, whereas $f_{ret}$ and $\beta_{\odot}^{\text{eff}}$ are mostly sensitive to the Galactic potential. To illustrate this, we calculate $f_{ret}$, $\eta_\odot$, and $\beta_{\odot}^{\text{eff}}$ assuming a surface SFR with a fixed radial scale length $R_d = 3$ kpc, i.e., $\Sigma \propto e^{-R/R_d}$, throughout the Galactic evolution and show the results in Figure 1(d). When compared to Figure 1(b), it can be seen clearly that whereas the evolution of $f_{ret}$ and $\beta_{\odot}^{\text{eff}}$ remain qualitatively unchanged, the evolution of $\eta_\odot$ changes dramatically. In this case, the evolution of $\eta_\odot$ becomes very similar to $f_{ret}$, which increases gradually instead of decreasing with time.

### 3. Impact of Natal Kicks on GCE Calculations

In order to study the impact of natal kicks on GCE, we use the one-zone chemical evolution code OMEGA (Côté et al. 2016, 2017b) to model the solar ring centered at $R_\odot = 8$ kpc with a thickness of 2 kpc, and take the closed box approximation for simplicity. The code tracks the contributions from low- and intermediate-mass stars taken from Karakas (2010) as well as massive stars taken from Kobayashi et al. (2006), where 50% of stars from 20–40 $M_\odot$ are assumed to explode as hypernovae. Stars with initial masses between 3–8 $M_\odot$ that form white dwarfs are assumed to contribute toward SNe Ia with a DTD $\propto t^{-1}$ and a minimum delay time of $\sim 40$ Myr (corresponding to the lifetime of an 8 $M_\odot$ star). The number of SNe Ia per unit of stellar mass formed, $N_{\text{Ia}}$, is fixed at $2 \times 10^{-3}$ with the SNe Ia yields adopted from Iwamoto et al. (1999) (for more details see Côté et al. 2016 and Ritter et al. 2018). A fraction $f_{\text{BNSM}}$ of massive stars are assumed to lead to BNSMs with a DTD $\propto t^\beta$ and ejecta mass $m_{\text{ej}}^{\text{BNSM}}$. We fix the value of $f_{\text{BNSM}} = 0.01$ and $\beta = -1$ when the effects due to natal kicks discussed in Section 2 are neglected. The effects of natal kicks are included by replacing $f_{\text{BNSM}}$ with $f_{\text{BNSM}} \times \eta_\odot(t)$, and $\beta$ with $\beta_{\odot}^{\text{eff}}(t)$, respectively. The yields of r-process elements in the BNSM ejecta are assumed to follow the solar r-process pattern for mass number $A \gtrsim 80$ from Arnould et al. (2007), which amounts to a mass fraction of $1.04 \times 10^{-3}$ for Eu in the ejecta. Because varying the ejecta mass $m_{\text{ej}}^{\text{BNSM}}$ amounts to an overall scaling of the Eu yield, its value is chosen between 0.01 and 0.02 $M_\odot$ that best fits to the data. We note that recent nucleosynthesis studies of BNSM outflows at different phases that produce a wide range of r-process nuclides (Wanajo et al. 2014; Just et al. 2015; Wu et al. 2016; Siegel & Metzger 2017) reported similar Eu yields as the number adopted here.

The SFR for $t \ge t_0$ in our calculation is adopted directly from Minchev et al. (2013) using the surface SFR at 8 kpc that is multiplied by the area of the ring. It is known from existing calculations (e.g., Côté et al. 2019) that with a constant value of $f_{\text{BNSM}}$, i.e., neglecting the effect of natal kicks, the curve of [Eu/Fe] has a flat trend with [Fe/H]. Thus, when $f_{\text{BNSM}}$ is

replaced by $f_{\text{BNSM}} \times \eta_\odot$, the curve of [Eu/Fe] is expected to follow the trend of $\eta_\odot$. Because the downward trend of [Eu/Fe] starts at [Fe/H] $\sim -0.8$, we assume that this coincides with $t = t_0$. For $t < t_0$, a constant value of SFR is chosen such that [Fe/H] reaches $\sim -0.8$ at $t = t_0$. Figure 2(a) shows the resulting SFR for $t_0 = 1$ Gyr (SFR1) and $t_0 = 2$ Gyr (SFR2). The initial gas mass is calculated by requiring [Fe/H] to reach 0.2 by the end of the evolution at $t = 13.7$ Gyr, which gives values of $10.5 \times 10^9 M_\odot$ and $9.4 \times 10^9 M_\odot$ for $t_0 = 1$ and 2 Gyr, respectively. We note here that the resulting evolution of $\alpha$ elements, such as Mg, with the above choice of values agrees well with the observed trend.

In Figure 2(b), we show the evolution of [Eu/Fe] as a function of [Fe/H], using both SFRs, for three different cases that illustrate effects due to natal kicks: (i) completely neglecting the effect of natal kicks, i.e., including neither $\eta_\odot(t)$ nor $\beta_{\odot}^{\text{eff}}(t)$ and using $\beta = -1$; (ii) including only $\eta_\odot(t)$ with $\beta = -1$; and (iii) including both $\eta_\odot(t)$ and $\beta_{\odot}^{\text{eff}}(t)$. When the effects of natal kicks are completely ignored, the flat trend of [Eu/Fe] for [Fe/H] $\gtrsim -1.5$ (see the dashed-dotted curves) consistent with the findings of Côté et al. (2019) is recovered.

In contrast, when only $\eta_\odot(t)$ is included with $\beta = -1$ in case (ii), it is clear that the trend of [Eu/Fe] follows that of $\eta_\odot(t)$, as expected. Specifically, [Eu/Fe] first increases with [Fe/H] for $t \le t_0$ (corresponding to [Fe/H] $\lesssim -0.8$), and then decreases for $t > t_0$ ([Fe/H] $\gtrsim -0.8$). The slope for the decreasing curve for [Fe/H] $\gtrsim -0.8$, however, is slightly flatter than the observed data. Finally, when both $\eta_\odot(t)$ and $\beta_{\odot}^{\text{eff}}(t)$ are included (solid curves) in case (iii), the smaller values of $\beta_{\odot}^{\text{eff}}(t) < -1$, i.e., steeper DTD, help to counter the increasing values of $\eta_\odot$ for [Fe/H] $\lesssim -0.8$ yielding a flat curve. On the other hand, for $t > t_0$, as the values of $\beta_{\odot}^{\text{eff}}$ continue to be lower than $-1$, it helps to steepen the slope of the [Eu/Fe] curve further that is primarily caused by the decreasing values of $\eta_\odot$. Overall, this leads to a very good agreement with the observed trend. We note that $m_{\text{ej}}^{\text{BNSM}} = 1.3 \times 10^{-2} M_\odot$ is chosen here to match the observed data for the case when both $\eta_\odot$ and $\beta_{\odot}^{\text{eff}}$ are included.

From the above discussion, it is evident that both $\eta_\odot$ and $\beta_{\odot}^{\text{eff}}$ are important for the evolution of [Eu/Fe]. In particular, the decreasing trend of $\eta_\odot$ for $t > t_0$ due to the inside-out formation of the MW disk is crucial in producing the decreasing trend in [Eu/Fe] for [Fe/H] $\gtrsim -0.8$. To reinforce this, we perform additional GCE calculations using the values of $\eta_\odot$ and $\beta_{\odot}^{\text{eff}}$ with a fixed scale length of 3 kpc shown in Figure 1(d), and show the resulting [Eu/Fe] evolution in Figure 2(c). In this case, the monotonically increasing $\eta_\odot$ results in an increasing trend of [Eu/Fe] when $\beta = -1$ is used. When $\beta_{\odot}^{\text{eff}}$ is included, it can at best counter the negative impact of $\eta_\odot$ to yield a flat curve for [Fe/H] $\gtrsim -0.8$ similar to calculations that neglect the effect of kick altogether (see Figure 2(c)).

### 4. Discussion and Conclusions

In this Letter we studied the effects of natal kicks of NSBs on the GCE of elements like Eu that are almost entirely produced by the r-process, assuming that BNSMs are the sole source of r-process with a standard DTD $\propto t^{-1}$. We find that natal kicks alter both the effective DTD and the occurring frequency of BNSMs in the solar neighborhood. In particular, the effect is amplified when the birth locations of NSBs are sampled according to an SFR that is consistent with the inside-out formation of the MW disk.
<figure><img src="image_5.png" /><figcaption>**Figure 2.** (a) Star formation rates SFR1 and SFR2. (b) Evolution of [Eu/Fe] with [Fe/H] for SFR1 and SFR2 for cases including $\eta_{\odot}$ and $\beta_{\odot}^{\text{eff}}$, including only $\eta_{\odot}$ ($\beta = -1$), and without including any kick effects. (c) Same as (b), but with $\eta_{\odot}$ and $\beta_{\odot}^{\text{eff}}$ from Figure 1(d), calculated using a fixed scale length of 3 kpc (see the text). In all cases $m_{\text{ej}}^{\text{BNSM}} = 1.3 \times 10^{-2} M_{\odot}$ is used. Observational data are adapted from Battistini & Bensby (2016).</figcaption></figure>

During the first ~1 Gyr of the disk formation, the solar vicinity gets a large enhancement in the effective BNSM frequency from NSBs that originate from the inner parts of the disk but merge within the solar neighborhood due to kick-induced migration. This enhancement decreases with time and

gives rise to the decreasing [Eu/Fe] trend for [Fe/H] $\gtrsim -0.8$. Additionally, natal kicks also lead to an effective DTD steeper than $\propto t^{-1}$, which further helps steepen the curve for [Eu/Fe] versus [Fe/H] for [Fe/H] $\gtrsim -0.8$ as well as flatten it for [Fe/H] $\lesssim -1$. When these two effects of natal kicks, namely, $\eta_{\odot}$ and $\beta_{\odot}^{\text{eff}}$, are combined together with the inside-out formation of the MW disk, the decreasing [Eu/Fe] at [Fe/H] $\gtrsim -0.8$ can be naturally reproduced without the need of additional r-process sources or modifications to the DTD $\propto t^{-1}$ for BNSM.

Because the exact evolution of [Eu/Fe] could depend on the assumed kick velocity distribution or the minimum merger time $t_{\text{merge}}^{\text{min}}$, we also performed additional calculations taking different values of $v_0$ and the minimal kick velocity, as well as the $t_{\text{merge}}^{\text{min}}$ to test the robustness of our conclusion. The corresponding results in Appendices A and B clearly show that the decreasing trend of [Eu/Fe] is not very sensitive to these parameters and the overall trend that matches the data well can be similarly reproduced. The effect of changing the minimum delay time for SN 1a is explored in Appendix C, which shows that the results are essentially unchanged.

Although we adopted a simple one-zone closed box model for the GCE using OMEGA for the solar vicinity, similar calculations using the same code have shown the mean trends are captured well when compared to more sophisticated GCE calculations (Côté et al. 2019). We also explored the effects of including outflow and inflow and found that they have a negligible impact on the results (see Appendix D). Nevertheless, detailed calculations for the chemodynamical evolution of the MW that take into account the migration of NSBs due to natal kicks are highly desirable. Such calculations, however, are likely to be computationally demanding and expensive.

An interesting consequence of natal kicks is that the evolution of $\eta(R, t)$ in particular has a strong radial dependence and is thus different for each ring. Figure 3 shows the evolution of $\eta$ and $\beta^{\text{eff}}$ for $t_0 = 2$ Gyr for different rings where the strong radial dependence of $\eta$ is evident in contrast to $\beta^{\text{eff}}$, which is roughly similar for all rings. As can be seen from the figure, the maximum value of $\eta_{\text{max}}$ that occurs at $t \approx t_0$ as well as the value of $\eta$ at any given time is lower for rings closer to the Galactic center. Because $\eta$ directly impacts the value of [Eu/Fe], its value is also expected to be lower for rings closer to the center and vice versa. Additionally, for $t > t_0$, the slope of $\eta$ becomes flatter (steeper) for rings closer (farther) to the Galactic center. Although the exact trend would depend on the details of the SFR and gas mass relevant for the ring, this is expected to result in a flatter (steeper) downward slope of [Eu/Fe] versus [Fe/H] for rings closer (farther) than the solar neighborhood. Such a prediction, in principle, can be verified if Eu is measured in a sufficient number of stars over a range of [Fe/H] at other locations of the MW.

Other effects that have already been shown to be important in explaining the metallicity distribution are radial gas flows and migration of stars in the disk (see, e.g., Schönrich & Binney 2009a, 2009b; Minchev et al. 2013; Schönrich & McMillan 2017). Although the details are complicated and beyond the scope of this Letter, the net result of such mixing is that it tends to weaken the radial metallicity gradient slightly. Clearly, this will also impact the [Eu/Fe] trend for the solar neighborhood due to some amount of mixing of stars and gas from the inner (outer) regions with a slightly flatter (steeper) [Eu/Fe] versus [Fe/H] curve from the inner (outer) region. On the other hand, because of the radial metallicity gradient, at any
<figure><img src="image_6.png" /><figcaption>**Figure 3.** (a) Evolution of $\eta$ with $t_0 = 2$ Gyr for rings at different distances from the center. (b) Corresponding evolution of $\beta^{\text{eff}}$.</figcaption></figure>

given time, the value of [Fe/H] will be higher in the inner region and lower in the outer region. Thus, for the solar ring, radial mixing will bring in stars with lower [Eu/Fe] from the inner regions, but with higher [Fe/H]. Exactly the opposite is applicable for stars coming from the outer region. Whereas the impact of stellar migration can only be calculated with detailed chemodynamical calculations, effectively, this would lead to some scatter around the mean trend (Tsujimoto & Baba 2019).

Finally, our results show that BNSMs with a DTD $\propto t^{-1}$ alone is sufficient to explain both the origin as well as the evolution of r-process elements in the Galaxy for [Fe/H] $\gtrsim -2$. This, however,

does not rule out other sources for r-process, but rather indicates that their contribution is likely subdominant compared to BNSM. We note, however, that this conclusion does not apply to the origin of r-process observed in very metal-poor stars formed in the early Galaxy where additional sources could still be required (Wehmeyer et al. 2015).

The authors thank Benoit Côte for help with using OMEGA, Ivan Minchev for his valuable inputs and comments, and the anonymous referee for a detailed and constructive report. M.-R. W. acknowledges support from the Academia Sinica by grant No. AS-CDA-109-M11, Ministry of Science and Technology, Taiwan under grant No. 108-2112-M-001-010, and the Physics Division, National Center of Theoretical Science of Taiwan. Z.Y. acknowledges the support from Special Funding for Advanced Users through LAMOST FELLOWSHIP and the Shanghai Sailing Program (Y955051001).

## Appendix A
Dependence on Natal Kick Velocity Distribution

The results presented in the main text are computed using the same distribution for the kick velocity $\propto \exp(-v/v_0)$ with $v_0 = 90 \text{ km s}^{-1}$ with a minimum value of $10 \text{ km s}^{-1}$. Here, we explore the effect of changing the distribution of kick velocity on the results. First, we repeat our calculations using $v_0 = 60 \text{ km s}^{-1}$ and $120 \text{ km s}^{-1}$. The results are qualitatively similar, as can be seen from Figure 4. The slope of [Eu/Fe] for [Fe/H] $\gtrsim -0.8$ is slightly flatter (steeper) for lower (higher) average kick velocity. The slight change in the slope is caused by slightly higher values of $\beta_{\odot}^{\text{eff}}$ for lower average kick velocity and vice versa. On the other hand, for lower average kick velocity, the NSBs migrate slightly less, but fewer of them escape from the potential. The two effects counterbalance each other such that $\eta_{\odot}$ remains roughly unchanged. Overall, the results are only weakly sensitive to the average kick velocity, with higher values resulting in even better fits to the observed data.

We further explore the effect due to the contribution of NSBs with lower kick velocities below $v \lesssim 10 \text{ km s}^{-1}$ that could arise from binaries that involve low-mass core-collapse SNe. We, again, repeat our calculations using $v_{\min} = 0 \text{ km s}^{-1}$ with $v_0 = 90 \text{ km s}^{-1}$. Compared to the default case with $v_{\min} = 10 \text{ km s}^{-1}$, the fraction of NSBs with $v \lesssim 20 \text{ km s}^{-1}$ is increased from ~10% to ~20%. As shown in Figure 5, this only leads to a marginally lower values of $\eta_{\odot}$ and $\beta_{\odot}^{\text{eff}}$ relative to the default case. Overall, the results qualitatively remain unchanged.
<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 4.** Effect of changing the distribution of $v_{\text{kick}}$. (a) Same as Figure 1(b), but with $v_0 = 60$ km s$^{-1}$. (b) Same as Figure 2(c), but with $v_0 = 60$ km s$^{-1}$. (c) Same as (a), but with $v_0 = 120$ km s$^{-1}$. (d) Same as (b), but with $v_0 = 120$ km s$^{-1}$ with $m_{\text{ej}}^{\text{BNSM}} = 1.35 \times 10^{-2} M_{\odot}$.</figcaption></figure>

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 5.** (a) Same as Figure 1(b), but with $v_{\text{min}} = 0$ km s$^{-1}$. (b) Same as Figure 2(c), but with $v_{\text{min}} = 0$ km s$^{-1}$.</figcaption></figure>

## Appendix B
### Dependence on Minimum BNSM Merger Time

We also explore the dependence of our results on the choice of $t_{\text{min}}^{\text{merge}}$, by repeating our calculations with $t_{\text{min}}^{\text{merge}} = 30$ Myr. Figure 6 shows that the corresponding results are very similar

to the calculations with $t_{\text{min}}^{\text{merge}} = 10$ Myr shown in Figure 2. As before, $\beta_{\odot}^{\text{eff}}$ helps to flatten the curve for [Fe/H] $\lesssim -0.8$ by countering the increasing values of $\eta_{\odot}$. Above [Fe/H] $\sim -0.8$, $\beta_{\odot}^{\text{eff}}$ acts in tandem with the $\eta_{\odot}$ to produce a negative slope for [Eu/Fe] that matches the observed data very well. Thus, the results are not sensitive to the choice of $t_{\text{min}}^{\text{merge}}$.
<figure><img src="image_13.png" /><figcaption>**Figure 6.** (a) Same as Figure 1(b), but with $t_{\text{merge}}^{\text{min}} = 30$ Myr. (b) Same as Figure 2(c), but with $t_{\text{merge}}^{\text{min}} = 30$ Myr.</figcaption></figure>

<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 7.** Same as Figure 2(c), but with minimum SN Ia delay time of 400 Myr with $N_{\text{Ia}} = 10^{-3}$.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Figure 8.** Same as Figure 2(c), but with inflow and outflow as described in Appendix D.</figcaption></figure>

# Appendix C

## Dependence on Minimum Delay Time of SN Ia

As shown in Hotokezaka et al. (2018), the minimum delay time of SN Ia can also affect the evolution of [Eu/Fe] versus [Fe/H]. In this appendix, we explore this effect by changing the minimum delay time from the default 40 Myr to 400 Myr, similar to the values considered in Hotokezaka et al. (2018). Figure 7 shows the results where the parameter $N_{\text{Ia}}$ is reduced to $10^{-3}$ in order to ensure the final value of [Fe/H] ~ 0.2 remains unchanged. The resulting curve for [Eu/Fe] versus [Fe/H] is almost unchanged compared to our default model that clearly shows that the results are not sensitive to the variation of the minimum delay time of SNe Ia.

# Appendix D

## Dependence on Gas Inflow and Outflow

As mentioned in the main text, our results are based on a closed box model for the GCE. In this appendix, we explore the effect of including the outflow and inflow of gas on the [Eu/Fe] versus [Fe/H] trend. The outflow is assumed to be proportional to the SFR, whereas the inflow rate is assumed to be

proportional to the outflow rate (see Côté et al. 2017b for details). The proportionality constant for the outflow rate and the ratio between the inflow and outflow rates are taken to 0.1 and 1, respectively. These values are consistent with numerical simulation of gas outflow and inflow of the Galactic disk (Kim & Ostriker 2018). We use the same initial gas masses of 10.5 $M_{\odot}$ and 9.4 $M_{\odot}$ for SFR1 and SFR2, respectively, as in the main text. The resulting trend of [Eu/Fe] versus [Fe/H] in this case, shown in Figure 8, is similar to the closed box calculations, indicating that inflow and outflow do not affect the results.

## ORCID iDs

Projjwal Banerjee https://orcid.org/0000-0002-6389-2697

Meng-Ru Wu https://orcid.org/0000-0003-4960-8706

Zhen Yuan https://orcid.org/0000-0002-8129-5415

## References

Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017a, PhRvL, 119, 161101

Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017b, ApJL, 848, L12

Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2019, PhRvX, 9, 031040
Argast, D., Samland, M., Thielemann, F. K., & Qian, Y. Z. 2004, A&A, 416, 997

Arnould, M., Goriely, S., & Takahashi, K. 2007, PhR, 450, 97

Battistini, C., & Bensby, T. 2016, A&A, 586, A49

Behroozi, P. S., Ramirez-Ruiz, E., & Fryer, C. L. 2014, ApJ, 792, 123

Beniamini, P., & Piran, T. 2016, MNRAS, 456, 4089

Bovy, J. 2015, ApJS, 216, 29

Chruslinska, M., Belczynski, K., Klencki, J., & Benacquista, M. 2018, MNRAS, 474, 2937

Côté, B., Belczynski, K., Fryer, C. L., et al. 2017a, ApJ, 836, 230

Côté, B., Eichler, M., Arcones, A., et al. 2019, ApJ, 875, 106

Côté, B., Fryer, C. L., Belczynski, K., et al. 2018, ApJ, 855, 99

Côté, B., O'Shea, B. W., Ritter, C., Herwig, F., & Venn, K. A. 2017b, ApJ, 835, 128

Côté, B., Ritter, C., O'Shea, B. W., et al. 2016, ApJ, 824, 82

Cowan, J. J., Sneden, C., Lawler, J. E., et al. 2019, arXiv:1901.01410

Cowperthwaite, P. S., Berger, E., Villar, V. A., et al. 2017, ApJL, 848, L17

Dominik, M., Belczynski, K., Fryer, C., et al. 2012, ApJ, 759, 52

Drout, M. R., Piro, A. L., Shappee, B. J., et al. 2017, Sci, 358, 1570

Fischer, T., Wu, M.-R., Wehmeyer, B., et al. 2020, ApJ, 894, 9

Fong, W., & Berger, E. 2013, ApJ, 776, 18

Frankel, N., Sanders, J., Rix, H.-W., Ting, Y.-S., & Ness, M. 2019, ApJ, 884, 99

Grichener, A., & Soker, N. 2019, ApJ, 878, 24

Griffen, B. F., Ji, A. P., Dooley, G. A., et al. 2016, ApJ, 818, 10

Hirai, Y., Ishimaru, Y., Saitoh, T. R., et al. 2015, ApJ, 814, 41

Hotokezaka, K., Beniamini, P., & Piran, T. 2018, IJMPD, 27, 1842005

Ishimaru, Y., Wanajo, S., & Prantzos, N. 2015, ApJL, 804, L35

Iwamoto, K., Brachwitz, F., Nomoto, K., et al. 1999, ApJS, 125, 439

Just, O., Bauswein, A., Pulpillo, R. A., Goriely, S., & Janka, H. T. 2015, MNRAS, 448, 541

Karakas, A. I. 2010, MNRAS, 403, 1413

Kasen, D., Metzger, B., Barnes, J., Quataert, E., & Ramirez-Ruiz, E. 2017, Natur, 551, 80

Kawaguchi, K., Shibata, M., & Tanaka, M. 2018, ApJL, 865, L21

Kim, C.-G., & Ostriker, E. C. 2018, ApJ, 853, 173

Kobayashi, C., Umeda, H., Nomoto, K., Tominaga, N., & Ohkubo, T. 2006, ApJ, 653, 1145

Komiya, Y., & Shigeyama, T. 2016, ApJ, 830, 76

Matteucci, F., Romano, D., Arcones, A., Korobkin, O., & Rosswog, S. 2014, MNRAS, 438, 2177

Mennekens, N., & Vanbeveren, D. 2014, A&A, 564, A134

Metzger, B. D. 2019, LRR, 23, 1

Miller, M. J., & Bregman, J. N. 2013, ApJ, 770, 118

Minchev, I., Chiappini, C., & Martig, M. 2013, A&A, 558, A9

Mösta, P., Roberts, L. F., Halevi, G., et al. 2018, ApJ, 864, 171

Ritter, C., Côté, B., Herwig, F., Navarro, J. F., & Fryer, C. L. 2018, ApJS, 237, 42

Safarzadeh, M., & Côté, B. 2017, MNRAS, 471, 4488

Safarzadeh, M., Ramirez-Ruiz, E., Andrews, J. J., et al. 2019, ApJ, 872, 105

Safarzadeh, M., & Scannapieco, E. 2017, MNRAS, 471, 2088

Schönrich, R., & Binney, J. 2009a, MNRAS, 396, 203

Schönrich, R., & Binney, J. 2009b, MNRAS, 399, 1145

Schönrich, R., & McMillan, P. J. 2017, MNRAS, 467, 1154

Schönrich, R. A., & Weinberg, D. H. 2019, MNRAS, 487, 580

Shen, S., Cooke, R. J., Ramirez-Ruiz, E., et al. 2015, ApJ, 807, 115

Siegel, D. M., Barnes, J., & Metzger, B. D. 2019, Natur, 569, 241

Siegel, D. M., & Metzger, B. D. 2017, PhRvL, 119, 231102

Tanaka, M., Utsumi, Y., Mazzali, P. A., et al. 2017, PASJ, 69, 102

Tauris, T. M., Kramer, M., Freire, P. C. C., et al. 2017, ApJ, 846, 170

Thornton, K., Gaudlitz, M., Janka, H. T., & Steinmetz, M. 1998, ApJ, 500, 95

Tsujimoto, T., & Baba, J. 2019, ApJ, 878, 125

Tsujimoto, T., & Shigeyama, T. 2014, A&A, 565, L5

van de Voort, F., Quataert, E., Hopkins, P. F., Kereš, D., & Faucher-Giguère, C.-A. 2015, MNRAS, 447, 140

Villar, V. A., Guilochon, J., Berger, E., et al. 2017, ApJL, 851, L21

Wanajo, S. 2018, ApJ, 868, 65

Wanajo, S., Sekiguchi, Y., Nishimura, N., et al. 2014, ApJL, 789, L39

Wehmeyer, B., Fröhlich, C., Côté, B., Pignatari, M., & Thielemann, F. K. 2019, MNRAS, 487, 1745

Wehmeyer, B., Pignatari, M., & Thielemann, F. K. 2015, MNRAS, 452, 1970

Winteler, C., Käppeli, R., Perego, A., et al. 2012, ApJL, 750, L22

Wu, M.-R., Barnes, J., Martínez-Pinedo, G., & Metzger, B. D. 2019, PhRvL, 122, 062701

Wu, M.-R., Fernández, R., Martínez-Pinedo, G., & Metzger, B. D. 2016, MNRAS, 463, 2323