Irina Dvorkin
email: dvorkin@iap.fr

Frédéric Daigne

Stephane Goriely

Elisabeth Vangioni

Joseph Silk

The impact of turbulent mixing on the galactic r-process enrichment by binary neutron star mergers

Keywords:

We study the enrichment of the interstellar medium (ISM) with rapid neutron capture (r-process) elements produced in binary neutron star (BNS) mergers. We use a semi-analytic model to describe galactic evolution, with merger rates and time delay distributions of BNS mergers consistent with the latest population synthesis models. In order to study the dispersion of the relative abundances of r-process elements and iron, we applied a turbulent mixing scheme, where the freshly synthesized elements are gradually dispersed in the ISM. We show that within our model the abundances observed in Milky Way stars, in particular the scatter at low metallicities, can be entirely explained by BNS mergers. Our results suggest that BNS mergers could be the dominant source of r-process elements in the Galaxy.

 Pian et al. 2017;Shappee et al. 2017). The high opacities and the total ejected mass of ∼ 0.05 M inferred from the observed spectral evolution of this event are consistent with model predictions (Li & Paczyński 1998;Metzger et al. 2010;Kasen et al. 2017). Moreover, r-process nucleosynthesis in BNS mergers was shown to produce a solar abundance pattern (Goriely et al. 2011;Bauswein, Goriely & Janka 2013;Just et al. 2015), at least for nuclei with mass number A > 90. On the other hand, under the hypothesis that all or most of the r-process material in the Universe is produced in BNS mergers, the observed r-process abundances can be used to infer the merger rate (Vangioni et al. 2016). The result of this estimate is indeed consistent with the merger rate measured directly from the two BNS mergers discovered in GW, which further supports the idea that most of the r-process elements were formed in these astrophysical sites.

However, a detailed comparison of the observed abundances with model predictions reveals a more complex picture. Population synthesis models generally predict long time delays between star formation episodes and BNS mergers, on the order of hundreds of Myr and up to several Gyr (Chruslinska et al. 2017;Belczynski et al. 2018;Giacobbo & Mapelli 2019). These long time delays come about because the initial separations of the newly formed NSs can be quite large, and energy loss to GW emission (and therefore orbital decay) is slow. This characteristic property of BNS mergers is difficult to reconcile with the obsevations of ultrametalpoor but r-process-rich stars in the MW halo and the Reticulum II and Tucana III ultrafaint dwarf galaxies (Ji et al. 2016;Hansen et al. 2017). This early enrichment seems to suggest a common origin of iron-group and r-process elements in these stars, as pointed out in several studies (Argast et al. 2004;van de Voort et al. 2015;Safarzadeh & Scannapieco 2017;Côté et al. 2018;Naiman et al. 2018;Safarzadeh, Sarmento & Scannapieco 2019;van de Voort et al. 2020), although note that some earlier models managed to reproduce observations with BNS mergers as the only source of rprocess elements (Shen et al. 2015;Komiya & Shigeyama 2016;Hirai et al. 2017). Moreover, Bonetti et al. (2019) show that the mass fraction of r-process elements that remains in dwarf galaxies can be as low as 10 -3 , further reducing their abundance in low-metallicity systems (although see Beniamini, Dvorkin & Silk 2018). The role of BNS mergers in high-metallicity environments is also a matter of debate. Côté et al. (2019) argue that the decreasing trend of r-process elements like Europium (Eu) relative to iron abundance observed in MW disc stars is inconsistent with the large time delays expected in the BNS scenario. On the other hand, Banerjee, Wu & Yuan (2020) argue that this problem can be resolved by including the effect of natal kicks of the NSs.

It is important to stress that in order to explain both the mean observed r-process abundance and its scatter relative to the iron-group element abundances, chemical evolution models need to include also a description of the dispersion of r-process elements in the interstellar medium (ISM), which can itself depend on the properties (rate, spatial distribution, explosion energy etc.) of r-process producing sites (Beniamini et al. 2018;Naiman et al. 2018;Brauer et al. 2019;Safarzadeh et al. 2019;Tsujimoto & Baba 2019).

In this paper, we study the contribution of BNS mergers to the total r-process element abundances in a MW-like galaxy using a semianalytic galaxy evolution model that includes a turbulent mixing scheme. We focus here on Eu since it is produced almost exclusively through the r-process, and is readily observable in the spectra of stellar atmospheres. Our model is based on the work described in Vangioni et al. (2016), where many of the uncertainties are studied in detail. In the present paper, we extend the model to a full semianalytic description of the growth of a MW-like galaxy, starting from small building blocks at early epochs. We then focus on the comparison between instantaneous and turbulent mixing, as well as the impact of the time-delay distribution of BNS mergers.

The structure of the paper is as follows. In Section 2, we discuss the main components of our model and the resulting BNS merger rates. We present the case of instantaneous mixing in Section 3, where we focus on the differences between various time delay distributions. Results for turbulent mixing, where both the mean abundances and the observed scatter are reproduced, are discussed in Section 4. The sensitivity to the model parameters is studied in Section 5. We conclude in Section 6.

M O D E L D E S C R I P T I O N

Merger trees of dark matter haloes

We use the modified GALFORM code (Cole et al. 2000;Parkinson, Cole & Helly 2008) to construct merger trees of dark matter (DM) haloes that host the galaxies. This algorithm implements the excursion set formalism, a key aspect of which is the conditional mass function: the fraction of mass f(M 1 |M 2 ) from haloes of mass M 2 at redshift z 2 that consisted of smaller haloes of mass M 1 at an earlier redshift z 1 . The conditional mass function is given by Lacey & Cole (1993):

f (M 1 |M 2 )d ln M 1 = 2 π σ 2 1 (δ 1 -δ 2 ) (σ 2 1 -σ 2 2 ) 3/2 exp - 1 2 (δ 1 -δ 2 ) 2 (σ 2 1 -σ 2 2 ) d ln σ 1 d ln M 1 d ln M 1 , (1)
where σ i = σ 2 (M i ) is the variance of the linear perturbation field smoothed on scale M i , and δ i is the critical density for spherical collapse at redshift z i . This conditional mass function is used to calculate the mean number of progenitors of mass M 1 at redshift z 1 + dz 1 of a halo of mass M 2 at z 2 = z 1 :

dN dM 1 = 1 M 1 df dz 1 z 1 =z 2 M 2 M 1 dz 1 .
(2)

In the modified GALFORM algorithm the following substitution is made:

dN dM 1 → dN dM 1 G(σ 1 /σ 2 , δ 2 /σ 2 ), ( 3
)
where G is referred to as a perturbing function, to be calibrated by comparison with N-body simulations. Parkinson et al. (2008) showed that by fitting the outcome of the algorithm to the results of the Millennium Simulation they obtained halo abundances that are consistent with the Sheth-Tormen mass function (Sheth & Tormen 2002). Thus, the perturbing function G expresses the uncertainty in the choice of the correct mass function. In this work, we use the following parametrization of the perturbing function: Parkinson et al. (2008). By specifying the required mass resolution M res , we can now calculate the mean number of progenitors with masses M 1 in the range M res < M 1 < M 2 /2:

G = G 0 (σ 1 /σ 2 ) γ 1 (δ 2 /σ 2 ) γ 2 , with the parameters G 0 = 0.57, γ 1 = 0.38, γ 2 = -0.01 taken from
P = M 2 /2 Mres dN dM 1 dM 1 (4)
as well as the fraction of mass of the final halo in unresolved objects:

F = Mres 0 dN dM 1 M 1 M 2 dM 1 . ( 5
)
Starting with the specified mass and redshift, the algorithm proceeds back in time. The redshift step dz 1 is chosen so that P 1 to ensure that the halo is unlikely to have more than two progenitors at the earlier redshift. The next step is to generate a uniform random number R⊂(0, 1). If R > P, no merger has occurred in this time-step and the halo mass is reduced to M 2 (1 -F; that is, the halo has grown by accretion of this unresolved mass rather than by mergers). If R ≤ P, a merger has occurred, and a random value of the progenitor mass M 1 is generated using the distribution in equation ( 2). The merging masses in this case are M 1 and M 2 (1 -F) -M 1 (there is still accretion). This process is repeated until the progenitor mass at some node reaches the resolution mass (and so cannot be further subdivided) or until the limiting redshift (typically z = 15) is reached.

Halo masses are stored only at specified redshifts. The final information is therefore contained in a list of haloes and the redshifts of their mergers. Specifically, at each node we obtain the masses of the merging haloes M 1 and M 1 as well as the mass accreted on to them and the final big halo M 2 :

M unres = M 2 -M 1 -M 1 . (6)
In the remainder of the paper, we build merger trees for a DM halo with mass 10 12 M at z = 0 in order to match a MW-like structure [START_REF] Binney | Galactic Astronomy[END_REF][START_REF] Bland-Hawthorn | [END_REF] and a resolution mass of 10 9 M . We checked that the results are converged at this resolution. The effect of the total DM halo mass is further discussed in Section 5.1.

Gaseous phases and star formation

As in other semi-analytic galaxy evolution models (e.g. All other processes are calculated as described in the text. The CGM acts as a reservoir for gas outflows from the galaxy.

obtained by solving a set of evolution equations on the branch (i.e. between mergers). The galaxy and its immediate surroundings are modelled as shown schematically in Fig. 1. The galaxy contains four reservoirs that correspond to the different phases: the circumgalactic medium (CGM), hot gas that occupies the halo, cold gas, and stars whose (time-dependent) masses are M CGM , M hot , M cold , and M * , respectively. The baryonic phases occupy a DM halo with mass M DM . Furthermore, each baryonic phase a is characterized by a metallicity Z a . The CGM acts as a reservoir that contains gas outflows. Its mass at the beginning of each time-step is defined as the mass contained in an annulus between the virial radius R vir of the halo and f CGM × R vir , where f CGM can be varied. The density profile outside the halo is assumed to be an extrapolation of the one inside (i.e. a Navarro-Frenk-White (NFW; Navarro, Frenk & White 1996) for DM, see below).

The gas flowing into the galaxy is equal to the sum of all the unresolved haloes that 'merged' with the galaxy during the time-step. For a time-step t step (determined by the merger tree resolution) the inflow rate is given by

Ṁinflow = f b M unres t step , ( 7
)
where f b is the baryon fraction and M unres is given by equation ( 6).

We assume that a fraction α in flows directly from the intergalactic medium (IGM) into the galaxy (for example, via cold flows), as shown by the solid red arrow in Fig. 1, and a fraction 1α in flows into the CGM, as shown by the red-dashed arrow.

The DM has an NFW profile with concentration c (Dutton & Macciò 2014):

ρ DM (r) = ρ s r r s -1 1 + r r s , ( 8
)
where r s = R vir /c, R vir is the virial radius, ρ s = M vir /4πr 3 s f (c) and f(c) = ln (1 + c)c/(1 + c). We take c = 10 for all the haloes in the tree.

The gas that flows into the galaxy is assumed to be immediately shock heated to the virial temperature T vir = GMm p /k B R vir and to have an isothermal profile:

ρ hot (r) = ρ 0 exp - 27 2 β 1 - ln (1 + r/r s ) r/r s (9) with β = 8πμm p Gρ s r 2 s 27k B T vir ( 10
)
where μ is the mean molecular weight of the gas, k B is the Boltzmann constant, m p the proton mass, and ρ 0 is a normalization given by the total gas mass in the halo.

Radiative cooling of the hot gas occurs on a time-scale:

t cool = 3ρ hot k B T vir 2μm p n 2 e (T vir , Z hot ) , ( 11
)
where is the cooling function (Sutherland & Dopita 1993) and n e is the electron density. Then, the rate of 'flow' of gas from the hot phase to the cold one is

Ṁcooling = 4π R 0 r 2 ρ hot t cool dr . ( 12
)
To model the rate at which stars form from cold gas, we use a power-law scaling with a cut-off at high gas masses (Kennicutt 1998):

ψ = A SFR M cold M 0,SFR γ e -M cold /Mq ( 13
)
where A SFR , α, and M q are to be specified. The cut-off at M q is meant to mimic the decline in the star formation rate (SFR) of the MW towards the present time (Snaith et al. 2015). This value for the damping mass is motivated by the observational and theoretical results (e.g. Bauer et al. 2013;Schaye et al. 2015); we further discuss this parametrization in Section 5.

The initial mass function of newly formed stars φ(m) is given by a Salpeter power law φ(m) ∝ m -p with p = 2.35 (Salpeter 1955) and normalized as follows:

msup m inf mφ(m) dm = 1 . ( 14
)
Stars return some of their mass into the cold phase in the form of winds, SN ejecta, and planetary nebulae. This 'recycled' mass flow is given by

e rec = msup m inf dm φ(m)ψ(t -τ (m, Z))(m -m rem ) ( 15
)
where Z is the metallicity, m is the stellar mass, m rem (m, Z) is the remnant mass, which we take from Fryer et al. (2012) and τ (m, Z) is the stellar lifetime, taken from Maeder & Meynet (1989) and Schaerer (2002).

In addition, SN-powered winds heat the cold gas at a rate

Ṁheating = 2 heating v 2 rot m f 8 M dm φ(m)ψ(t -τ (m, Z))E kin (m), ( 16
)
where E kin (m) is the energy released by the explosion of a star of mass m taken from Woosley & Weaver (1995) and v rot is the rotation velocity of the DM halo. A fraction f ex of this reheated gas is expelled out of the halo as galactic wind.

To summarize, the flow of gas between the different phases is given by the following equations: Note that

ṀCGM = (1 -α in ) Ṁinflow + f ex Ṁheating (17) Ṁhot = α in Ṁinflow -Ṁcooling + Ṁheating (1 -f ex ) (18) Ṁcold = Ṁcooling -Ṁheating -ψ + e rec ( 19
) Ṁ * = ψ -e rec ( 20
)
ṀCGM + Ṁhot + Ṁcold + Ṁ * = Ṁinflow . ( 21
)
Model parameters used in our calculation are listed in Table 1. We explore variations in DM halo mass (and stellar mass) as well as mass and time resolution in Section 5.

With the parameters summarized in Table 1, our reference model for a 10 12 M DM halo results in M * = 10 10 M in stars at z = 0, close to the bulge mass in the MW estimated from observations (Licquia & Newman 2015). In addition, the SFR at z = 0 is 1.22 M yr -1 in our model, compared with the observed value of 1.65 M yr -1 in the MW. The stellar mass in our model is thus lower than the total stellar mass in the MW (around 6 × 10 10 M ). We further discuss this difference in Section 5.1, where we show that our main conclusions are insensitive to the total stellar mass at z = 0. We note here that the metallicity distribution function (MDF) in our model provides a reasonably good fit to the global MDF of the MW (see Section 5.1).

Formation and merger rates of binary neutron stars

The evolution of individual binary systems is best described with population synthesis models that follow the evolution of binary massive stars (e.g. Andrews et al. 2015;de Mink & Belczynski 2015;Belczynski et al. 2018;Chruslinska et al. 2018;Artale et al. 2019;Chruslinska, Nelemans & Belczynski 2019;Giacobbo & Mapelli 2019). The results of these models are then implemented in galaxy evolution models to obtain the merger rate. Since this calculation involves many different spatial and temporal scales, each model necessarily involves various simplifications in the description of BNS formation or galaxy evolution. The outcome of the most recent models is generally in agreement with the local rate at z = 0 of 1090 +1720 -800 Gpc -3 yr -1 inferred from the GW observations of two BNS mergers GW170817 and GW190425 (Abbott et al. 2017a(Abbott et al. , 2020)). Note that BNS can also form in globular clusters (Ye et al. 2020) or nuclear clusters (Fragione et al. 2019), however, the formation rates in these channels are expected to be several orders of magnitude below the observed value.

Another observational constraint can be inferred from the observations of the Galactic double binary pulsars (Kim, Perera & McLaughlin 2015;Pol, McLaughlin & Lorimer 2020). The latest estimate of R g = 37 +24

-11 Myr -1 from Pol et al. (2020) is consistent with the local rate measured from GW observations, although note that the translation of this Galactic rate into the local rate requires an understanding of the dependence of the BNS merger rate on galactic properties.

In this work, we use a phenomenological approach to describe BNS formation rate. At each time-step of the calculation and for a given stellar mass at birth m and metallicity Z (chemical evolution is discussed below), we calculate the remnant mass m rem (m, Z) using the results of Fryer et al. (2012). Namely, the remnant mass is calculated separately for stars with initial masses below 11 M and 11 M < m < 30 M (we assume that all stars with initial masses above 30 M collapse to a BH). The lower mass limit for NS formation is taken to be m lower NS = 9 M + 0.9 log 10 (Z met /Z ) M ( 22)

for metallicities of Z met > 10 -3 Z and m lower NS = 6.3 M otherwise. Thus, all stars in the initial mass range of m lower NS < m < 11 M form an NS. The mass remnant of stars in the initial mass range of 11 M < m < 30 M is given (in the delayed explosion scenario) by:

m rem = 1.1 + 0.2e (m-11)/4 -(2 + Z met )e 0.4(m-26) ( 23
)
where both m and m rem are in units of the solar mass. We then take any remnant with m rem < 3 M to be an NS. Below m lower NS all stars become white dwarfs (WDs).

Then, the formation rate of NSs is given by

R NS (t) = msup m inf NS (m, Z)ψ(t -τ (m, Z))φ(m) dm , ( 24
)
where NS (m, Z) is the probability to form an NS, taken here to be 1 if the identity of the remnant is an NS according to the fitting function in Fryer et al. (2012) where t min and t max are specified below. Note that the factor α BNS accounts for the formation of binaries.

In the following, we consider three models for the delay time distribution P(t delay ; see Fig. 2). Our Model A corresponds to the standard 1/t normalized distribution, i.e.

P (t delay ) = At -1 delay , t delay ∈ (t min , t max ), (26)
where A = 1/log (t max /t min ). This distribution is characterized by a mean delay time given by t delay = (t maxt min )/log (t max /t min ). The choice t max = t H 13.8 Gyr is a natural cut-off since these models are not evolved past the current age of the Universe. Note that an extrapolation of this distribution to even longer delay times can be counterbalanced by simply decreasing the total merging fraction α BNS . This functional form is a generic choice in population studies, as it is often obtained in population synthesis models (e.g. Dominik et al. 2012).

Due to the crucial character of the delay time distribution at early times, i.e. in the vicinity of t min , we also consider here another type of distributions that allows for a more detailed account of BNS mergers with small delay times without the artificial cut-off at t = t min . This Model B is described by the functional

P (t delay ) = B ln t delay /t min t delay , t delay ∈ (t min , t max ) ( 27
)
with a normalization constant B = 2/(ln (t max /t min )) 2 leading to a mean time delay of t delay = Bt max (ln (t max /t min ) + t min /t max -1). This functional form provides a good fit to the recent population synthesis model of Chruslinska et al. (2018), at least for t delay 100 Myr, and has the particularity to include a soft cut-off at low t delay . This specific ln (t)/t behaviour of the distribution allows us to consider a non-negligible fraction of BNS mergers at early times, as discussed below without affecting the 1/t law found by population synthesis models (see Fig. 2). Finally, we also implemented a time delay distribution in tabulated form, taken directly from fig. 5 in the recent detailed population synthesis model of Chruslinska et al. (2018), which we refer to as Model C. Note that the actual minimal delay time remains uncertain in such an approach and may be an artefact of the binning (linear in time) in Chruslinska et al. (2018). We adopt for this model an arbitrary value of t min = 50 Myr.

All three models are compared in Fig. 2 assuming the same value of t max = 13.8 Gyr and a value of t min tuned, so that the three distributions have the same mean delay time of t delay = 2.37 Gyr. As seen in Fig. 2, despite their similar mean time delay, these models can have a very different behaviour at early times, hence a different fraction of mergers with very short time delays. Such a fraction of systems with a time delay t delay ≤ t * is given by

F delay (< t * ) = t * tmin P (t delay ) dt delay . ( 28
)
For example, for the parameters shown in Fig. 2, the fraction of systems with delay times below t * = 100 Myr amounts to 0.14 for Model A, 0.27 for Model B, and ∼0.11 for Model C. As shown in Vangioni et al. (2016; see in particular their fig. 3), population synthesis models (Belczynski, Kalogera & Bulik 2002) predict that a significant fraction of BNS mergers takes place with short time delays, and more specifically that about 20-25 per cent of systems may merge with time-scales lower than 100 Myr. Such a significant fraction of fast merging systems is properly described by our Model B shown in Fig. 2 but is underestimated by Models A and C. At each time-step t, we calculate the total mean number of BNS mergers N m = R merger t using equation 25. In the instantaneous mixing case (described in Section 3), we use N m in our chemical evolution equations (see Section 2.4). However, as BNS mergers are rare events, N m can be small and even much smaller than 1 in many time-steps, especially at high redshifts for low-mass haloes. This led us to develop a more realistic approach, namely the turbulent mixing model, where we draw the number of actual mergers from a Poisson distribution with mean N m As a consequence, for certain time-steps there are no mergers at all. Then, for each of these mergers, we draw its position in the galaxy, assuming a 3D Gaussian distribution with variance equal to the half-light radius R e . The latter is computed as R e = 0.02R vir (Zanisi et al. 2020). As an output of this calculation, we obtain a list of BNS mergers, their host haloes, and the time and spatial position of each merger. The results for our turbulent mixing case are described in Section 4.

Finally, in both cases we assume that each BNS merger produces a (constant) mass m r of r-process elements, and more specifically a Eu mass m Eu = 10 -4 M , as predicted by realistic BNS merger simulations of both the dynamical and wind ejecta (Just et al. 2015). Note that the only real free parameter governing the evolution of Eu in our model is the product α BNS m Eu We checked that the mean evolution is insensitive to changes in α BNS , m Eu that keep their product constant and discuss in Section 4 the impact of α BNS on the dispersion in the turbulent mixing scenario.

The merger rates in our MW-like galaxy model, for the three time delay distributions discussed above, are shown in Fig. 3. We note that in our model, all the small haloes that form at high redshifts merge with each other or fall on to the larger halo, and thus act as building blocks for the final, MW-like structure. Therefore, each point in Fig. 3 corresponds to a specific (sub-)halo and shows the instantaneous BNS merger rate at the cosmic age that corresponds to its formation (colour coded).

We also compare our results with the new constraints from Galactic double NSs (Pol et al. 2020; the magenta line) and two population synthesis models: Chruslinska et al. (2018; the black vertical line) and Artale et al. (2019; the colour-coded lines). For the comparison with the model in Artale et al. (2019), which gives the merger rate as a function of stellar mass and redshift, we compute the merger rate for each subhalo in our model with stellar mass M * > 10 7 M (the limit given by Artale et al. 2019) and for two redshifts z = 0.1 and 2 (colour coded). The result in Chruslinska et al. (2018) is given only for z = 0, and we plot it for M * = 10 10 M . We note that our results are slightly below the value observed from Galactic double NSs at z = 0, and the corresponding value in the Chruslinska et al. (2018) model. On the other hand, our results are in general agreement with Artale et al. (2019), given as a function of stellar mass. In particular, we observe that the merger rate for a given stellar mass slightly increases with redshift. We note that we made no particular attempt to fit the exact physical parameters of the MW as discussed in Sections 2.2 and 5.

Chemical evolution

First, we outline the equations of chemical evolution for elements produced in stars, which follow those in Daigne et al. (2004Daigne et al. ( , 2006)), although note the differences due to different definitions and the introduction of gas heating and cooling.

For each element i, the rate at which it is ejected from stars to the cold gas phase is

e i (t) = msup m inf dmφ(m)ψ(t -τ (m, Z))m ej i (m, Z), ( 29
)
where m ej i (m, Z) is computed from stellar yieldstaken from Woosley & Weaver (1995). In addition, BNS mergers produce Eu with the rate: e Eu (t) = m Eu R merger (t) where R merger (t) is given in equation ( 25) and m Eu is the mass in Eu produced in each merger (see Table 1). In the instantaneous mixing model, we treat this contribution in the same way as the other species, whereas in the turbulent mixing case we follow the scheme described in detail below. We assume that Eu is produced only in BNS mergers.

We define by X i, hot , X i, cold , X i, CGM , and X i, IGM the mass fraction of element i in the cold, hot, CGM, and IGM phases, respectively. We assume that matter in-falling into the galaxy, Ṁinflow , does not contain metals. The evolution of the metal mass fractions is consequently given by Ẋi,cold (t

) = 1 M cold Ṁcooling X i,hot -X i,cold -X i,cold e rec + e i , ( 30
) Ẋi,hot (t) = 1 M hot Ṁheating (1 -f ex ) X i,cold -X i,hot + α in (X i,IGM -X i,hot ) Ṁinflow , and (31)
Ẋi,CGM (t) = 1 M CGM Ṁheating f ex X i,hot -X i,CGM + (X i,IGM -X i,CGM )(1 -α in ) Ṁinflow . ( 32
)
We calculate the evolution of iron abundance as in Vangioni et al. (2016); more specifically, we consider two sources of iron, i.e. corecollapse SNe and Type Ia SNe. For the latter, we assume that a fraction α WD = 0.09 of the WDs formed in our model ends up as Type Ia SNe, where the delay time between WD formation and SNIa explosion, t d, wd , is distributed as P (t d,wd ) ∝ t -1 d,wd with t min, wd = 100 Myr. We assume that each SNIa produces 0.5 M in iron. A discussion on the impact of these assumptions can be found in Vangioni et al. (2016). With the model parameters summarized in Table 1, we obtain a rate of SNIa of 0.4 per century at z = 0, consistent with observations in the MW [START_REF] Matteucci | The Origin of the Galaxy and Local Group[END_REF]. The evolution of metallicity is further discussed below.

Our model allows us to reconstruct the chemical enrichment history of each galaxy. Indeed, the identity of each halo, its stellar mass and metallicity are saved at various time-steps along the merger tree. Then, by going backwards in the tree and tracing the times of past mergers we can reconstruct the time of formation of a certain stellar mass and subsequently its metallicity.

Turbulent mixing of r-process elements

In order to model the production of r-process elements in BNS mergers and their subsequent diffusion in the surrounding ISM, we use the turbulent mixing scheme proposed by Hotokezaka, Piran & Paul (2015) and recently, Beniamini & Hotokezaka (2020). These authors propose an analytic description of the turbulent diffusion process that leads to the mixing of metals in the ISM. This mixing proceeds in two stages: first, the newly synthesized metals expand with the blast wave until the end of the Sedov-Taylor phase. We neglect this short-lived phase. Then, metals continue to diffuse due to turbulence in the ISM. The turbulent diffusion coefficient is given by (Beniamini & Hotokezaka 2020)

D = αc s H = α c s 10km s -1 H 100 pc kpc 2 Gyr -1 , ( 33
)
where c s is the sound speed in the ISM, H is the scale height, and α is the efficiency parameter. Typical diffusion length at time t is then given by √ Dt. Since we do not describe the 3D shape of the galaxy, in the following we treat the diffusion coefficient D directly as the free parameter of the model that we assume to be constant in all structures. Below, we explore value of D in the range of 0.001-0.01 kpc 2 Gyr -1 .

By including the turbulent diffusion effect, Beniamini & Hotokezaka (2020) obtained promising results for the scatter in the r-process element abundances observed in the MW, although for their main results they considered the case of negligible delay times between formation and merger of the BNS, constant SFR and focused on metallicity values of [Fe/H] > -1. In this work, we implement turbulent mixing in our full galaxy evolution model and show that it can explain the observed scatter across the entire metallicity range.

Our calculation proceeds as follows. As described in Section 2.3, at each time-step t of the calculation we obtain a list of BNS mergers (drawn from a Poisson distribution with mean N m = R merger t, see equation 25), their host galaxies, and the (3D) positions r i and times t i of all the events. We assume that each of these mergers produces a fixed amount of Eu m Eu , which then disperses in the ISM. If at a later time a star forms in the vicinity of the merger event, it will be enriched by the Eu found in the ISM. If the star forms farther away and/or much later than the BNS merger, the abundance of Eu it will inherit will correspond to the mean value at its location. Taking into account all the contributions from BNS mergers in a given galaxy, we can express the Eu density at a given location, hereafter referred to as 'observer point', r obs at time t obs by (see equation 2 in Beniamini & Hotokezaka 2020)

ρ Eu ( r obs , t obs ) = i m r exp -| r obs -r i | 2 4D•(t obs -t i ) (4πD • (t obs -t i )) 3/2 , ( 34
)
where D is the diffusion coefficient, r i and t i are the positions and times of the BNS merger, respectively, and the sum runs over all the mergers in the given galaxy such that t obs > t i . Note that we use a simplified version of the expression given in Beniamini & Hotokezaka (2020), in particular we assume spherical diffusion. We take a number of observer points proportional to the stellar mass formed along the branch, between mergers. We take a maximum of 500 observer stars per halo (we verified that the results are converged in this limit). Since the calculation involves random draws for the number of mergers, their positions and times, each realization of the turbulent mixing scheme for a given merger tree provides a different evolution of Eu abundance, contrary to the homogeneous case, as we further discuss below. We note that we do not employ the turbulent mixing scheme for SNII or SNIa, but rather assume that the metals produced by all SNe are immediately mixed in the entire galaxy. This assumption is justified by the fact that these events are very frequent relative to BNS mergers.

Galaxy mergers

During a merger event n smaller galaxies fall on to a larger galaxy.

We assume that the stars and DM are simply added up, while the cold gas of the smaller galaxies is heated and added to the hot gas reservoir:

M DM = M DM,host + n j =1 M DM,j , ( 35
)
M CGM = M CGM,host + n j =1 M CGM,j , ( 36
)
M hot = M hot,host + n j =1 M hot,j + M cold,j , and (37)
M * = M * ,host + n j =1 M * ,j . ( 38
)
In the turbulent mixing scheme for r-process elements abundance, we assume that, after each merger event, all the gas of the in-falling galaxies is completely mixed, and in particular the r-process elements are distributed evenly in the whole ISM. This mean abundance (immediately after the merger) is calculated under the instantaneous mixing approximation discussed above and therefore takes into account mass-loss due to galactic outflows.

E U RO P I U M A B U N DA N C E : I N S TA N TA N E O U S M I X I N G

In this section, we explore the Eu abundance as a representative rprocess element. As discussed above, we build merger trees for a halo of 10 12 M that results in a stellar mass of M * = 10 10 M at z = 0, and compare our results with the compilation of Eu and Fe abundances from the SAGA data base 1 (Suda et al. 2008) where we took all the stars in the MW and dwarf galaxies present in the data base and removed duplicates. We assume throughout a constant BNS merger fraction of α BNS = 4 × 10 -3 , m Eu = 10 -4 M . Our chemical evolution model results in a metallicity of [Fe/H] -0.1 at z = 0, slightly higher than Vangioni et al. (2016) but still lower than the value reached by disc stars, which we further discuss in the following Section. We also verified that the evolution of α-elements in our model is consistent with observations in the entire metallicity range.

In this section, we assume an instantaneous mixing of Eu in the ISM and focus on the predictions obtained with different time delay distributions. Even in this instantaneous mixing model, some scatter in the [Eu/Fe] versus [Fe/H] plane is expected, especially at early times, due to the fact that small haloes can have very different chemical enrichment histories.

Fig. 4 shows the results for Model A (i.e. the t -1 time delay distribution) characterized by different values of t min . Each point corresponds to a single halo. The increase of the cut-off time t min clearly sets the galactic Eu enrichment to higher metallicity, as already widely studied in many galactic chemical evolution models 1 http://sagadata base.jp/ Figure 4. Evolution of the Eu abundance in all the sub-haloes in the tree with time delay distribution given by Model A: P(t) ∝ 1/t. The blue circles are obtained with t min = 1 Myr, the orange triangles with t min = 10 Myr, and the green diamonds with t min = 100 Myr. In all cases, t max = 13.8 Gyr is adopted. Observations from the SAGA data base are shown by the grey crosses. 4 but for Model B time delay distribution (P(t) ∝ log (t/t min )/t). The results are qualitatively similar to those in Fig. 4. [START_REF] Matteucci | [END_REF]Côté et al. 2019;Simonetti et al. 2019). In addition, we observe that longer time delays result in slightly larger dispersion in the [Eu/Fe] versus [Fe/H] plane (the mean delay times for the three cases considered here are 1.4, 1.9, and 2.8 Gyr for t min = 1, 10, and 100 Myr, respectively) and that the scatter increases towards low metallicities. The general shape of the evolution observed in Fig. 4 is qualitatively consistent with the conclusions of van de Voort et al. ( 2020), who used an hydrodynamical simulation to model a MW-like galaxy and a phenomenological description of BNS formation similar to the one considered here, in particular a t -1 time delay distribution in their fiducial model. However, we note that van de Voort et al. ( 2020) obtained a significantly larger scatter than the one shown in Fig. 4, as discussed below.

Next, we explore other time delay distributions, not treated in previous works on r-process abundances. The results for Model B (log (t)/t time delay distribution) are shown in Fig. 5 and are quite similar to those obtained with Model A. In this case, the dispersion at low metallicities and for large t min is slightly larger than in the case of Model A, but the mean value is significantly lower.

Finally, in Fig. 6 we compare models A, B, and C, keeping t max = 13.8 Gyr and varying t min so as to obtain the same mean value of t delay . It is clear from Fig. 6 that Eu abundance as a function of metallicity depends primarily on this mean value t delay , especially towards higher metallicites. There are, however, significant differences at the low-metallicity end: Model B appears to be much more efficient in enriching small haloes at an early stage. This feature can be traced back to our comparison of the time delay models in Fig. 2. As discussed above, for the same mean value t delay , in Model B there are about 27 per cent of mergers with a delay time below 100 Myr, compared with only ∼ 11 per cent in Model C. These very rapid mergers efficiently enrich the small haloes at early times. Indeed, if we choose values of t min that give the same fraction of rapid mergers in all three models, as shown in Fig. 7 (where rapid is defined as a delay below 100 Myr, i.e. F delay (< 100Myr) = 0.11, see Eq.28), we obtain very similar evolution at low metallicities, as expected. However in this case the models slightly diverge towards higher metallicities, where mergers with long time-scales become important.

We also note the evolution of Eu abundance is determined by the product α BNS m Eu . We verified that the results remain unchanged if both α BNS and m Eu are varied so as to keep their product constant.

To summarize, we draw two conclusions from Figs 4-6. First of all, we observe that our model, which includes BNS mergers as the only sources of r-process elements, reproduces some of the Eu abundances at low metallicities, even though the mean time delays between stellar formation and BNS mergers are of the order of a few Gyrs. This result depends on the BNS merger rate in high-redshift low-mass galaxies, which is at present not constrained observationally (see Fig. 3).

Our second major conclusion is that we are not able to reproduce the large observed scatter in the Eu abundance under the assumption of an instantaneous mixing, even by taking into account the formation history of the galaxy from small building blocks. Below, we will study this dispersion with an analytic diffusion model. The results for the three models are very similar for the lowest metallicity points, but slightly diverge at the higher metallicity end.

E U RO P I U M A B U N DA N C E : T U R B U L E N T M I X I N G

In the following, we explore the consequences of the turbulent mixing scheme, described in Section 2.4.1. As discussed above, in this case we obtain multiple values of [Eu/Fe] abundance for each sub-halo, corresponding to the different 'observer' points (which mimic the stars that form at different locations in the galaxy). For each subhalo, we therefore obtain a distribution of [Eu/Fe] and [Fe/H] values. We bin these values in 20 metallicity ([Fe/H]) bins of width 0.3 dex and show the results in Figs 8-12. Note that, in this approach, we have an additional parameter, which is the diffusion coefficient D. In the following, we explore how this diffusion parameter affects our results.

Since the total number, locations and times of the mergers in each halo are drawn from probability distributions, the evolution of Eu abundance varies for different realizations for a given merger tree. Below, we show one realization for each choice of the model parameters, but we checked that these results are representative in that both the median values in each metallicity bin and the percentiles do not change significantly among the different realizations.

As can be seen in Figs 4 and8, the median values of our turbulent mixing model do not necessarily coincide with those shown in the previous section within the instantaneous mixing regime. The reason for this is that in the instantaneous mixing case, the average was performed over the whole ISM in the galaxy. In the turbulent mixing scheme, we first calculate the Eu and Fe abundances of our 'observer' points, and then average over those.

In the upper panel of Fig. 8, we show our results for Model A with a diffusion coefficient of D = 0.005 kpc 2 Gyr -1 : the 2D-binned data in the [Eu/Fe]-[Fe/H] plane (the orange heatmap) and the 5-95 percentiles in each bin of [Fe/H] (the red lines). We observe that the distribution in our model is in reasonably good agreement with the data up until [Fe/H]∼-0.5. The most interesting feature of this plot are the high [Eu/Fe] values reached even at very low metallicities: these are contributed by observer stars that happen to be born in a recently enriched region. Thus, our semi-analytic model coupled to a turbulent diffusion scheme confirms the conclusions of Shen et al. (2015) obtained with a numerical model of a MW-like galaxy. Beniamini & Hotokezaka (2020) who introduced this turbulent diffusion scheme, used it in the context of a simplified galaxy model and focused on high metallicities, where they obtained a much higher dispersion than observed. Here, we show that this model, when coupled to a detailed galaxy evolution calculation, reproduces the observed abundances in the low to intermediate metallicity ranges probed by observations. The model predicts extremely Eurich (although rare) stars even at [Fe/H] ∼ -4, however, there are very few observations at these metallcities, and in the following we will concentrate on the range [Fe/H] > -3.5.

However, the model shown in the upper panel of Fig. 8 does not reproduce the decline in [Eu/Fe] observed at high metallicities for stars that reside in the disc. Indeed, there is currently a debate in the literature, since Côté et al. (2019) argued that BNS mergers cannot reproduce this observational result and could therefore not be the only source of r-process elements, while Banerjee et al. ( 2020) recently showed that this problem can be alleviated by including the effect of natal kicks of the NSs. Our model lacks spatial resolution and is therefore not well adapted to explore this population of stars, but we can roughly identify them with the stars born along the last branch of our merger tree, i.e. the stars born after the last major merger event, which in our model takes place at the cosmic age of 5.5 Gyr. To investigate the impact of the disc geometry, we slightly change our prescriptions after this last major merger event: we assume that all the stars in the last branch of the merger tree are born in the disc with scale height h = 400 pc (Bland-Hawthorn & Gerhard 2016). On the other hand, the BNS mergers, even those occurring at later times, can take place anywhere, including the halo, since their progenitor stars formed at much earlier epochs. Therefore, we assume that only BNS that merge in the disc contribute to its chemical enrichment (these BNS could have formed before or after the last major merger). This last assumption limits the number of mergers that contribute to r-process elements in the disc, which leads to the decline in [Eu/Fe] as seen in the lower panel of Fig. 8, in much better agreement with the data at high metallicities. Precisely, the transition (last major merger event at 5.5 Gyr) corresponds to a metallicity of [Fe/H] = -0.6 (indicated by the black vertical lines in Fig. 8). This value slightly depends on our assumptions for the Type Ia supernova rate, but a lower limit is [Fe/H] = -1, i.e. the contribution of SNII alone to the iron budget at this age.

This geometrical effect is therefore a promising avenue to understand the observed evolution of [Eu/Fe] at high metallicities, although other mechanisms discussed in the literature, such as a metal-dependent efficiency of forming BNS (Simonetti et al. 2019), could lead to a similar effect. However, we stress that we do not attempt to describe all the observed properties of the disc, such as the differences in the ages and metallicities of MW disc stars (Hayden et al. 2017) and their variation as a function of radius (Toyouchi & Chiba 2018). Clearly, a full spatial model of the different components of the Galaxy is needed in order to describe the chemical evolution of the disc. We also stress that, as can be seen in Fig. 8, the results for metallicities below [Fe/H] = -0.6, which are our main subject in this work, are not affected by this 'disc' prescription, as expected. In all the plots below the results are obtained with this new prescription.

Next, we explore the effect of the diffusion coefficient D. In Fig. 9, we show the predictions for our Model A with t min = 10 Myr and varying D in the range 0.001 -0.01 kpc 2 Gyr -1 . Both values of D = 0.005 and 0.01 kpc 2 Gyr -1 are in equally good agreement with the data, and we conclude that these values are consistent with the observed scatter.

Finally, we study the effect of the BNS merger time delay distribution. In Fig. 10, we compare the results for Models A and B, adopting t min = 10 Myr and D = 0.005 kpc 2 Gyr -1 in both cases. We find that the dispersion does not change significantly, strengthening our conclusion that the scatter stems primarily from the diffusion process and therefore essentially depends on the diffusion coefficient D. On the other hand, the median value of the relative [Eu/Fe] abundance is lower for Model B, as already discussed above (see Figs 4 and5). With a minimal delay time t min = 10 Myr, Model B is characterized by a mean delay time t delay = 3.3 Gyr, significantly larger than the one characterizing Model A of 1.9 Gyr. We conclude that, keeping all the other model parameters fixed, the time delay distribution mainly affects the mean abundances and the threshold metallicity at which the r-enrichment takes place, while the diffusion coefficient D affects the scatter.

In Fig. 11, we compare Models A, B, and C where the respective t min parameter is adjusted so as to give the same value of the mean For both models, the diffusion coefficient is D = 0.005 kpc 2 Gyr -1 and the minimal delay time is t min = 10 Myr. Model A is characterized by a mean delay time t delay = 1.9 Gyr and Model B by 3.3 Gyr.

t delay (see Fig. 6). We observe that there is very little difference between these models except for the lowest metallicity bins, similar to the instantaneous mixing case shown in Fig. 6. Similarly, in Fig. 12 we compare Models A, B, and C where the respective t min parameter is adjusted so as to give the same fraction of rapid mergers (F delay (< 100Myr) = 0.11; see Fig. 7). We find that models A and C produce very similar results, but the Eu abundances in Model B are lower because of the larger (on average) merger times.

As mentioned above, the evolution of Eu abundance in the homogeneous mixing model is determined by the product α BNS m Eu . This is also the case for the median value of [Eu/Fe] turbulent mixing model. We find nevertheless that in the lowest metallicity bins there are less Eu-enriched stars if α BNS is lower.

S E N S I T I V I T Y TO M O D E L PA R A M E T E R S

In this section, we perform a sensitivity study in order to check that our main conclusions are sufficiently robust. We focus on two aspects of the model: DM (and stellar) mass and the resolution of the merger tree.

Galaxy mass: DM halo and stars

As we discussed above, our reference model reproduces several important observable features of the MW (Section 2.2), but it has slightly smaller stellar mass and it does not reach supersolar metallicities. One may wonder therefore whether our main results would hold in a model with, for example, a higher stellar mass. Our fiducial DM merger tree is calculated for a z = 0 halo of M halo = 10 12 M , which corresponds to current estimates (Smith et al. 2007;[START_REF] Bland-Hawthorn | [END_REF]. However, there is a significant uncertainty in this value, the different estimates varying by a factor of a few, depending on the technique used. Below, we explore two additional models with M halo = 5 × 10 11 M and M halo = 2.6 × 10 12 M , corresponding to the minimal and maximal halo masses from table 8 of [START_REF] Bland-Hawthorn | [END_REF]. Naturally, the choice of DM halo mass at z = 0 affects all the other masses in the model, including the mass of stars. Specifically, while our fiducial model reaches a stellar mass of M * (z = 0) = 10 10 M , the model with the smallest DM halo reaches M * (z = 0) = 6.3 × 10 9 M and the one with the largest DM halo 6.3 × 10 10 M . Note that the latter case is fully consistent with the stellar mass of the MW. In this latter case we also modified the SFR cut-off mas M q (see Table 1) to M q = 2 × 10 10 M to allow for star formation at the later epochs of galaxy evolution. The parameter M q used in our SFR law in equation ( 13) is a free parameter of the model. In this study, we did not attempt to fit a specific star formation history; we therefore leave the determination of this parameter and the resulting uncertainty to future work.

First, we explore the MDF in our three test models. In Fig. 13, we compare our results with observational data from the Pristine Survey, which explored the MDF of the MW halo (Youakim et al. 2020). Specifically, the data for all stars are taken from their Fig. 1. As can be seen in Fig. 13, all three models are in reasonable agreement with the total observed distribution. Therefore, we conclude that our model is adequate to describe halo stars. However, we caution that none of the models seems to faithfully reproduce the disc stars, as is evident from Fig. 14, where we plot the oxygen abundance as a function of metallicity: there is no visible change of slope at high metallicities. This difficulty to describe the high-metallicity population is potentially related to the lack of spatial resolution in our model, as already mentioned in Section 4.

In Figs 15 and 16, we compare the abundances in Eu and Fe for the three DM halo masses, either assuming our homogeneous or our turbulent mixing scheme. In both cases, the three models produce almost identical results, despite the differences in total stellar masses, In all cases we use the time delay distribution P(t) ∝ 1/t with t min = 10 Myr (Model A above).

as expected for relative abundances. The main difference we observe is in the model with the largest halo mass, which reaches higher metallicities than the other models. From Fig. 15, we see that the models with halo masses of 5 × 10 11 , 10 12 , and 2.6 × 10 12 M reach metallicities of [Fe/H] = -0.12, -0.09, and 0.29, respectively.

Based on this analysis we conclude that our main results, in particular the amount of scatter in [Eu/Fe], are not sensitive to the total stellar mass at z = 0 in our model.

Merger tree resolution

As mentioned above, the fiducial mass resolution for our merger tree was taken to be M res = 10 9 M . While we cannot go to much lower masses for computational reasons, we performed a test where the resolution was reduced to M res = 7 × 10 8 M . Another resolutionrelated parameter is the number of tree levels that are stored for post-processing. As a reminder, we first build the DM merger tree backwards in time with a very small time-step that ensures that we capture all the haloes with masses above the resolution mass. However, we only store haloes at specific output times, specifically at N level = 19 for our fiducial model. We then evolve the model forward in time, adding gas to the haloes and solving the equations for gas inflow/outflow, SFR, etc., as described in Section 2.2. These evolution equations are solved 'along the branch' between the times of saved output. When we reach the time of saved output, that means that a merger occurred, the smaller haloes fall into the bigger one and we assume the ISM is immediately mixed. Therefore, increasing the number of saved levels N level would, on the one hand, allow a more detailed description of all the subhaloes, but on the other hand introduce a more frequent mixing on the scale of each galaxy. To test the influence of our assumptions, we ran a model with N level = 29 levels.

The results for these resolution tests are shown in Figs 17, 18, and 19. While the stellar mass differs slightly among these models (within a factor ∼3), it converges to the same value by z = 0, and the abundances of Eu and Fe are not at all sensitive to these changes. In particular, there is no difference between the models with N level = 19 versus N level = 29 merger tree levels.

D I S C U S S I O N

The results presented above suggest that BNS mergers can be the dominant r-process producing site, and in particular our model can reproduce both the mean value and the scatter in [Eu/Fe] observed in low to intermediate metallicity stars in the MW, while a very crude 'disc' prescription allows also to reproduce the high-metallicity part. We also showed that the two key model parameters that govern this evolution are the time delay distribution of BNS mergers and the diffusion coefficient. These conclusions are in contrast with some of the previous studies, according to which BNS mergers are not able to produce r-process elements in sufficient amounts at early times (Argast et al. 2004;van de Voort et al. 2015;Côté et al. 2018;Naiman et al. 2018;van de Voort et al. 2020), but in qualitative agreement with Shen et al. (2015), Komiya & Shigeyama (2016), and Hirai et al. (2017).

What could be the source of this discrepancy between different models? As a representative case, we compare our model with van de Voort et al. (2020). In this study, the authors use an hydrodynamical model of a MW-like galaxy and a phenomenological model for BNS mergers. To calculate the BNS merger rate they define a constant BNS merger efficiency per unit of stellar mass formed of 3 × 10 -6 M -1 . In our model the equivalent quantity R merger / is not constant (see equations 13 and 25), but varies with redshift and among sub-haloes, ranging from 10 -7 to 10 -4 M -1 . We also observe that in our model this quantity increases with metallicity, but at a given metallicity the scatter reaches about two orders of magnitudes. van de Voort et al. ( 2020) define a time delay distribution for the BNS mergers with P (t delay ) ∝ t -1 delay , as in our Model A, with t min = 30 Myr. The merger rate at z = 0 for their fiducial model is 1.7 10 2 Myr -1 at the high end of the value inferred from Galactic double NSs, whereas we obtain 10 Myr -1 , slightly lower than the low end observational value. We also note that we use m Eu = 10 -4 M for the mass of Eu produced in each BNS merger, while van de Voort et al. ( 2020) use m Eu = 4.7 10 -5 M (but note that the mean evolution in our model is only determined by the product α BNS m Eu ). Another key difference between our models is the SFR, which in our model is always slightly below the one in the MW, while van de Voort et al. ( 2020) use values a factor of a few above the measured values in the MW.

To conclude, we listed above the main potential sources of differences, but a comprehensive comparison between the different models present in the literature is beyond the scope of this study. We suggest that the differences are manifested in the merger rate of BNS in small, low-metallicity galaxies, and that r-process abundances can indeed be used as an indirect probe of these environments, but we leave a more complete analysis to future work.

An additional observational feature that could help distinguish between the different scenarios is the decrease in the scatter in [Eu/Fe] towards intermediate metallicities (Cescutti et al. 2015). This effect could be caused by the fact that sources with different lifetimes (fast/slow BNS mergers and rare SNe) contribute to the r-process enrichment.

C O N C L U S I O N S

The question of the origin of r-process elements, while witnessing major advances in the past few years, still remains open. On the one hand, the spectacular multimessenger observations of the coalescing BNS 170817 (Abbott et al. 2017a,b;Coulter et al. 2017;Chornock et al. 2017) has confirmed the existence of kilonovae and the associated production of r-process elements in such systems. On the other hand, chemical evolution models struggle to reproduce the r-process abundance observed in MW stars, both ultrametal-poor stars as well as young disc stars.

In this paper, we studied the abundance of r-process elements with a semi-analytic galactic chemical evolution model, focusing on BNS mergers. We explored several assumptions regarding the time delay distribution of BNS mergers, and implemented a turbulent mixing scheme to study the dispersion of r-process elements in the ISM. Our main conclusions can be summarized as follows:

(i) BNS mergers can be a dominant r-process producing site, and in particular both the merger rates and time delays predicted by population synthesis models are compatible with the Eu abundances observed in MW stars.

(ii) The assembly history of our model galaxy from small building blocks has only a weak impact on the scatter in the Eu abundances.

(iii) Turbulent mixing of the freshly synthesized elements in the ISM leads to a significant dispersion in the [Eu/Fe]-[Fe/H] plane at low metallicities, in good agreement with the observed scatter for a moderate diffusion coefficient.

We investigated several functional forms of the BNS merger time delay distribution, and found that the shape of the distribution at the lower end plays a critical role in the resulting r-process abundances in low-mass, low-metallicity haloes. This result underlines the importance of studying fast-merging BNS scenarios with population synthesis models. More multimessenger BNS detections may provide new constraints in the future on such a fast-merging population as the afterglow should be easier to detect with the denser close environment expected in this case (Duque et al. 2020).

We also wish to stress several caveats to our model, that may have an important effect on the final results. First of all, we do not take into account the birth kicks of the NSs. Indeed, newly born NSs can acquire velocities of up to hundreds of km/sec. As a result, not only the binary may be potentially disrupted, which is already taken into account in the low value of α BNS , but even if it survives it may find itself far away from its birth site. The r-process elements produced during the merger will therefore be released in a metallicity environment very different from the one the binary itself was born in (Banerjee et al. 2020). However, we expect that this effect will not have a strong influence on the fast-merging BNS that enrich the low-metallicity environments (Andrews & Zezas 2019).

The second important caveat concerns the mass released in rprocess elements in each BNS merger event. In this work, we assumed that the NS masses were constant, and each BNS merger produced the same mass in Eu. By comparing the kilonova associated with GW170817 and a sample of short GRB afterglow light curves, Gompertz et al. (2018) find some evidence for diversity in the emission properties of kilonovae that could not be explained by the viewing angle alone. This could point to differences in the ejected mass, which depends on the masses of the binary components (Shibata & Hotokezaka 2019). Diversity is especially expected in the case of BH-NS mergers, the contribution of which is neglected in this study (see the discussion in Vangioni et al. 2016), as their rates are expected to be significantly smaller than those of BNS (Abbott et al. 2019).

We do not include the interactions of a central massive BH with the surrounding gas that can cause galactic outflows and affect the abundances in the ISM. This effect can be important, especially in view of the growing evidence that feedback from central BHs can play a role in dwarf galaxies (Koudmani, Henden & Sijacki 2020).

Finally, we emphasize that we did not attempt to produce an accurate model of the MW, and indeed our model galaxy has lower stellar mass and fewer high-metallicity stars than the MW. However, as discussed in Section 5, our main results, and in particular the amount of scatter in [Eu/Fe], are not sensitive to the stellar mass at z = 0. Better understanding the precise MDF of the MW and specifically the enrichment of disc stars would require a more detailed treatment of the various Galactic components where our model lacks spatial resolution. This is also the case for understanding the observed decreasing trend of [Eu/Fe] with [Fe/H] at high metallicity, which are mainly located in the disc. To illustrate this aspect, we included in our model a simple prescription to restrict the turbulent mixing in the disc at late epochs and managed to reproduce the observed evolution of [Eu/Fe] at [Fe/H] > -0.6.

In view of these caveats, we conclude that the role of BNS mergers as sources of r-process elements may be more prominent than according to other works, but is still not entirely clear. In future work, we plan to include other possible sources associated with rare types of core-collapses. Moreover, further progress in the modeling of fast-merging BNS is needed in order to fully understand r-process abundances in low-metallicity environments.

Figure 1 .

 1 Figure 1. Evolution of the galaxy along a merger tree branch. The infall from IGM into structures is given by the merger tree (shown by the red arrows).All other processes are calculated as described in the text. The CGM acts as a reservoir for gas outflows from the galaxy.

Figure 2 .

 2 Figure 2. Time delay distributions between formation and merger of the BNS. Models A, B are analytical, while Model C uses the data from fig. 5 in Chruslinska et al. (2018). The value of t min for Models A and B were chosen to give the same mean time delay as Model C with t min = 50 Myr, i.e. t delay = 2.37 Gyr.

Figure 3 .

 3 Figure 3. BNS merger rates versus stellar mass in the subhalos in the tree (the coloured dots) with time delay distribution according to Model A with t min = 10 Myr. The colour code corresponds to the age of the subhalo (see colourbar). The vertical lines show the range of predictions in the models of Chruslinska et al. (2018) for an MW-like galaxy (plotted here for M * = 10 10 M ) and the galactic rates (Pol et al. 2020) based on the observation of the double pulsar. Also shown are the predictions computed from the model in Artale et al. (2019) for z = 0.1 and z = 2.

Figure 5 .

 5 Figure 5. Same as Fig. 4 but for Model B time delay distribution (P(t) ∝ log (t/t min )/t). The results are qualitatively similar to those in Fig. 4.

Figure 6 .

 6 Figure 6. Comparison of [Eu/Fe] versus [Fe/H] predictions obtained for the three time delay distribution models given in Fig. 2, namely P(t) ∝ 1/t with t min = 44 Myr (Model A; the blue circles), P(t) ∝ log (t/t min )/t with t min = 0.4 Myr (Model B; the orange triangles), and Chruslinska et al. (2018;Model C; the green diamonds). The mean delay time is the same in all three models, namely t delay = 2.37 Gyr. Observed values are shown by the grey crosses. The results for the three models are very similar, except for the lowest metallicity points.

Figure 7 .

 7 Figure 7.Comparison of the three time delay models with t min values chosen to produce the same fraction of rapid mergers F delay (< 100Myr) = 0.11: P(t) ∝ 1/t with t min = 55 Myr (Model A; the blue circles), P(t) ∝ log (t/t min )/t with t min = 10 Myr (Model B; the orange triangles), andChruslinska et al. (2018; Model C; the green diamonds). Observed values are shown by the grey crosses. The results for the three models are very similar for the lowest metallicity points, but slightly diverge at the higher metallicity end.

Figure 8 .

 8 Figure 8. Eu abundance in the turbulent mixing model with the P(t) ∝ 1/t time delay distribution (Model A) with a minimal time delay of t min = 10 Myr (see Fig. 4) and diffusion coefficient D = 0.005 kpc 2 Gyr -1 . The 5 per cent and 95 per cent percentiles are shown by the red lines. Panel (a): results with the turbulent mixing scheme as described in Section 2.4.1. The model reproduces early enrichment in Eu, but overproduces Eu at high metallicities [Fe/H] -0.6. Panel (b): same model with the inclusion of a special treatment of the last branch of the merger tree mimics the evolution of disc stars and results in a better agreement with observations (see text for details). The black vertical lines at [Fe/H] = -0.6 mark the beginning of the last branch. The small differences between the panels at low metallicity are only due to different realizations of the turbulent mixing model.

Figure 9 .

 9 Figure 9. Same as Fig. 8 with t min = 10 Myr and varying diffusion coefficients (indicated in the legend). The dotted lines show the 5 per cent and 95 per cent percentiles of each case, while the solid lines show the medians. The scatter at low metallicities decreases with D while the median is not affected.

Figure 10 .

 10 Figure 10. Same as Fig. 8, comparing Model A (blue) and Model B (orange) for the delay time distribution. The dashed lines show the 5 per cent and 95 per cent percentiles of each model, while the solid lines show the medians.For both models, the diffusion coefficient is D = 0.005 kpc 2 Gyr -1 and the minimal delay time is t min = 10 Myr. Model A is characterized by a mean delay time t delay = 1.9 Gyr and Model B by 3.3 Gyr.

Figure 11 .

 11 Figure11. Same as Fig.8, comparing predictions based on Model A (blue), B (orange), and C (green) for the time delay distribution, where the respective values for t min are chosen so as to give the same mean delay time t delay = 2.37 Gyr (see Fig.6). The dashed lines show the 5 per cent and 95 per cent percentiles of each model. For all models, the diffusion coefficient is D = 0.005 kpc 2 Gyr -1 .

Figure 12 .

 12 Figure 12. Same as Fig. 8, comparing predictions based on Model A (blue), B (orange), and C (green) for the time delay distribution, where the respective values for t min are chosen so as to give the same fraction of rapid mergers F delay (< 100Myr) = 0.11. The dashed lines show the 5 per cent and 95 per cent percentiles of each model. For all models, the diffusion coefficient is D = 0.005 kpc 2 Gyr -1 .

Figure 13 .

 13 Figure 13. MDF with varying DM halo masses at z = 0 (and different stellar masses, see text). Blue: reference model (see Table 1) with M DM (z = 0) = 10 12 M , M * (z = 0) = 10 10 M . Orange: M DM (z = 0) = 5 × 10 11 M , M * (z = 0) = 6.3 × 10 9 M . Green: M DM (z = 0) = 2.6 × 10 12 M , M * (z = 0) = 6.3 × 10 10 M . Data are taken from the Pristine Survey (Youakim et al. 2020).

  Figure 13. MDF with varying DM halo masses at z = 0 (and different stellar masses, see text). Blue: reference model (see Table 1) with M DM (z = 0) = 10 12 M , M * (z = 0) = 10 10 M . Orange: M DM (z = 0) = 5 × 10 11 M , M * (z = 0) = 6.3 × 10 9 M . Green: M DM (z = 0) = 2.6 × 10 12 M , M * (z = 0) = 6.3 × 10 10 M . Data are taken from the Pristine Survey (Youakim et al. 2020).

Figure 14 .

 14 Figure 14. Comparison in the [O/Fe]-[Fe/H] plane between models with different DM haloes at z = 0. The blue circles: reference model (see Table 1) with M DM (z = 0) = 10 12 M , M * (z = 0) = 10 10 M . The orange triangles: M DM (z = 0) = 5 × 10 11 M , M * (z = 0) = 6.3 × 10 9 M . The green crosses: M DM (z = 0) = 2.6 × 10 12 M , M * (z = 0) = 6.3 × 10 10 M .

  Figure 14. Comparison in the [O/Fe]-[Fe/H] plane between models with different DM haloes at z = 0. The blue circles: reference model (see Table 1) with M DM (z = 0) = 10 12 M , M * (z = 0) = 10 10 M . The orange triangles: M DM (z = 0) = 5 × 10 11 M , M * (z = 0) = 6.3 × 10 9 M . The green crosses: M DM (z = 0) = 2.6 × 10 12 M , M * (z = 0) = 6.3 × 10 10 M .

Figure 15 .

 15 Figure 15. Comparison between models with different DM halo mass at z = 0 (and different stellar masses, see text) for our homogeneous mixing scheme.In all cases we use the time delay distribution P(t) ∝ 1/t with t min = 10 Myr (Model A above).

Figure 16 .

 16 Figure 16. Same as Fig. 15 within the turbulent mixing model. The solid and dashed lines show the median and 5-95 percentiles, respectively, for each model.

Figure 17 .

 17 Figure 17. Convergence of the stellar mass as a function of the merger tree parameters.

Figure 18 .

 18 Figure 18. Convergence in the [Eu/Fe]-[Fe/H] plane as a function of the merger tree parameters: homogeneous mixing. In all cases, we use the time delay distribution P(t) ∝ 1/t with t min = 10 Myr (Model A above). The blue circles: reference model with 19 tree levels, M res = 10 9 M . The orange triangles: 29 tree levels. The green crosses: M res = 7 × 10 8 M .

Figure 19 .

 19 Figure 19. Same as Fig. 18 within the turbulent mixing model. The solid and dashed lines show the median and 5-95 percentiles, respectively, for each model.

Table 1 .

 1 Parameters used in the galaxy evolution model (see text for details).

	Name	Description	Value
	α in	IGM fraction falling directly into galaxy	0.8
	A SFR	SFR efficiency (M yr -1 )	0 . 0 3
	M 0, SFR	Pivot mass for SFR (M )	1 0 9
	M q	Cut-off mass for SFR (M )	8 × 10 9
	γ	SFR slope: dln SFR/dln M cold	2.0
	heating	Heating efficiency	0.01
	m inf	Lower limit for stellar masses (M )	0 . 1
	m sup	Upper limit for stellar masses (M )	100
	f ex	Fraction of reheated gas as galactic winds	0.05
	α WD	Fraction of WD that merge as SNIa	0.09
	t min, WD	Minimal time delay of SNIa (Myr)	100
	m Fe	Fe mass produced in each SNIa (M )	0 . 5
	α BNS	Fraction of NS in merging binaries	4 × 10 -3
	m Eu	Eu mass produced in each merger (M )	1 0 -4

MNRAS 506, 4374-4388 (2021) Downloaded from https://academic.oup.com/mnras/article/506/3/4374/6325186 by CNRS user on 05 May 2023

This paper has been typeset from a T E X/L A T E X file prepared by the author.MNRAS 506, 4374-4388 (2021) Downloaded from https://academic.oup.com/mnras/article/506/3/4374/6325186 by CNRS user on 05 May 2023

AC K N OW L E D G E M E N T S

We thank the referee Gabriele Cescutti for constructive comments that helped to improve the manuscript. The authors acknowledge support from the Centre National de la Recherche Scientifique -International Research Project (CNRS-IRP) project 'Origin of heavy elements in the Universe: compact stars and nucleosynthesis'. ID acknowledges support from the 'Tremplins nouveaux entrants' program of Sorbonne Université. FD, ID, and EV acknowledge financial support from the Centre National d' Études Spatiales (CNES). SG acknowledges financial support from FNRS (Belgium). This work was partially supported by the Fonds de la Recherche Scientifique -FNRS and the Fonds Wetenschappelijk Onderzoek -Vlaanderen (FWO) under the EOS Project No O022818F.

DATA AVA I L A B I L I T Y S TAT E M E N T S

The data underlying this article will be shared on reasonable request to the corresponding author.