Semi-analytic modelling of the europium production by neutron star mergers in the halo of the Milky Way

TL;DR

This study models europium production in the Milky Way's halo via neutron star mergers, highlighting the impact of neutron star kick velocities on r-process element distribution and comparing results with observations.

Contribution

It introduces a semi-analytic model incorporating neutron star kick velocities to study r-process enrichment in the Galactic halo.

Findings

High neutron star kick velocities reduce europium abundance in halo stars.

Low kick velocities produce europium levels consistent with observations for certain metallicities.

Instantaneous mixing assumption may be inadequate for modeling r-process enrichment.

Abstract

Neutron star mergers (NSM) are likely to be the main production sites for the rapid (r-) neutron capture process elements. We study the r-process enrichment of the stellar halo of the Milky Way through NSM, by tracing the typical r-process element Eu in the Munich-Groningen semi-analytic galaxy formation model, applied to three high resolution Aquarius dark matter simulations. In particular, we investigate the effect of the kick velocities that neutron star binaries receive upon their formation, in the building block galaxies (BBs) that partly formed the stellar halo by merging with our Galaxy. When this kick is large enough to overcome the escape velocity of the BB, the NSM takes place outside the BB with the consequence that there is no r-process enrichment. We find that a standard distribution of NS kick velocities decreases [Eu/Mg] abundances of halo stars by $\sim 0.5$~dex compared to models where NS do not receive a kick. With low NS kick velocities, our simulations match observed [Eu/Mg] abundances of halo stars reasonably well, for stars with metallicities [Mg/H]$\geq -1.5$. Only in Aquarius halo B-2 also the lower metallicity stars have [Eu/Mg] values similar to observations. We conclude that our assumption of instantaneous mixing is most likely inaccurate for modelling the r-process enrichment of the Galactic halo, or an additional production site for r-process elements is necessary to explain the presence of low-metallicity halo stars with high Eu abundances.

Equations

$$M_{i/o}={\overline{M}}_{\mathrm{b}}(M_{\mathrm{vir}},z)-\sum _{\mathrm{i}}{M}_{\mathrm{b},\mathrm{i}}$$

$${\overline{M}}_{\mathrm{b}}(M_{\mathrm{vir}},z)=\frac{f_{\mathrm{b}}{M}_{\mathrm{vir}}}{{\left[1+0.26M_F(z)/M_{\mathrm{vir}}\right]}^3}$$

$$\Delta M_{\mathrm{reheated}}=\frac{4}{3}ϵ\frac{{\eta }_{\mathrm{SN}}{E}_{\mathrm{SN}}}{V_{\mathrm{vir}}^2}\Delta M_{\ast }$$

$$M_{\mathrm{vir}}=\frac{4}{3}\pi R_{\mathrm{vir}}^{3}200\frac{3H^2}{8\pi G}=\frac{100H^2{R}_{\mathrm{vir}}^3}{G}$$

$${\dot{M}}_{\mathrm{re}-inc.}={\gamma }_{\mathrm{eg}}\cdot M_{\mathrm{eg}}\frac{V_{\mathrm{vir}}}{R_{\mathrm{vir}}}$$

$${\dot{M}}_{\mathrm{bh},\mathrm{r}}={\kappa }_{\mathrm{AGN}}\left(\frac{M_{\mathrm{bh}}}{10^8{M}_{\odot }}\right)\left(\frac{M_{\mathrm{hg}}}{0.1M_{\mathrm{vir}}}\right){\left(\frac{V_{\mathrm{vir}}}{200{\mathrm{km}\mathrm{s}}^{-1}}\right)}^3.$$

$$\Delta M_{\mathrm{bh},\mathrm{q}}=f_{\mathrm{bh}}\left(\frac{m_{\ast }+m_{\mathrm{cg}}}{M_{\ast }+M_{\mathrm{cg}}}\right)\frac{M_{\mathrm{cg}}}{1+(280{\mathrm{km}\mathrm{s}}^{-1}/V_{\mathrm{vir}}{)}^2}$$

$${\Sigma }_{\mathrm{crit}}(M_{\odot }{\mathrm{pc}}^{-2})=0.59V_{\mathrm{vir}}({\mathrm{km}\mathrm{s}}^{-1})/R_{\mathrm{disc}}(\mathrm{kpc})$$

$$e_{\mathrm{burst}}={\beta }_{\mathrm{burst}}{\left(\frac{m_{\ast }+m_{\mathrm{cg}}}{M_{\ast }+M_{\mathrm{cg}}}\right)}^{{\alpha }_{\mathrm{burst}}}$$

$$\rho (r)=\frac{{\rho }_0}{\left(r/a\right){\left(1+r/a\right)}^2}.$$

$$v_{\mathrm{esc}}(0)=\frac{\sqrt{2}}{0.465}{V}_{\max },$$

Full text

††pubyear: 2017††pagerange: Semi-analytic modelling of the europium production by neutron star mergers in the halo of the Milky Way–References

1 Introduction

The rapid (r-) neutron capture process which leads to the formation of roughly half of the elements heavier than iron has been known since the classical paper of Burbidge et al. (1957). Yet, there is an ongoing debate about its astrophysical production sites. The two main candidates are neutron star mergers (NSM) and core-collapse supernovae (SNe), where for the latter scenario a distinction is made between the prompt explosions of massive stars with masses in the range $8-10\ M_{\odot}$ (eg. Wheeler et al., 1998) and the delayed explosions of very massive stars ($\geq 20\ M_{\odot}$, eg. Woosley et al., 1994). Those with a mass between 10 and 20 $M_{\odot}$ are not likely candidates because chemical evolution models predict little to no scatter in [Eu/Fe] abundance ratios for the stars born out of their ashes, which clearly contradicts observations (Ishimaru & Wanajo, 1999). Whereas neutrino-driven winds associated with delayed explosions seem too proton rich to be the production sites of the high-mass r-process elements (Arcones et al., 2007; Arcones & Thielemann, 2013), prompt explosions might not actually exist in reality (Argast et al., 2004). Other suggested r-process production sites are long-duration gamma-ray bursts (see e.g. Metzger et al., 2008), magnetic proto-neutron star winds (Suzuki & Nagataki, 2005) and supernovae fallback (Fryer et al., 2006) amongst others.

This paper focusses on the NSM scenario for r-process enrichment, which is arguably the most likely r-process production site, because of the following reasons:

• •

NSM nucleosynthesis models robustly produce both heavy and light r-process nuclei (Goriely et al., 2011; Wanajo et al., 2014).

• •

The 244Pu interstellar medium density as a function of time is incompatible with continuous r-process production, as would be the case in the core-collapse SNe production scenario, but instead points towards a low rate / high yield production site, and perfectly matches the NSM scenario (Wallner et al., 2015; Hotokezaka et al., 2015).

• •

The recent direct detection of gravitational waves from a NSM and associated signal emitted by radioactive (r-process) material (Abbott et al., 2017a, b) favor the NSM scenario for r-process enrichment.

• •

Ultra faint dwarf galaxies and the Milky Way are likely to have the same production mechanism for r-process elements, triggered by rare events (Beniamini et al., 2016b).

However, a combination of the two different production sites could also explain the origin of r-process elements in our Galaxy (eg. Ji et al., 2016b).

We study the r-process enrichment by using europium (Eu) as its tracer, because this stable r-process element is one of the r-process elements that can be observationally detected without exceptionally high data quality (eg. Frebel, 2010). In Figure 1 we show the scatter in [Eu/Mg] as a function of [Mg/H], for stars in the Galactic halo and in nearby dwarf galaxies, made after querying the VizieR database (Frebel, 2010). Mg is chosen as a metallicity tracer to compare the Eu-abundance with, instead of the more commonly used Fe, to obtain an abundance distribution that is directly comparable to our model predictions (as we will explain in the next sections). At low metallicities (eg. [Mg/H]$<-2.5$) some halo stars have high Eu abundances, a feature that should be explained by any chemical evolution model of the r-process.

For many years it has been assumed that NSM alone are unable to reproduce the [Eu/Mg] ratios for stars with [Mg/H]$<-2.5$ (Argast et al., 2004; Wehmeyer et al., 2015), i.e. that Eu enrichment through NSM cannot take place at these low metallicities, mainly because NSM would be too rare and their average coalescence times too large. Also Matteucci et al. (2014) and Cescutti et al. (2015) showed that neutron star binaries can only be the dominant production sites of r-process elements if they merge after 1 Myr or in less than 10 Myr, respectively. However, population synthesis studies of neutron star binaries suggest a wide distribution of delay times proportional to $1/t$, with the first NSM occurring after 10 Myr or 100 Myr and $t$ the time since the formation of the neutron star binary (Belczynski et al., 2006). However, a small fraction of double neutron stars can merge on even smaller timescales of 1-10 Myr (e.g. Beniamini et al., 2016a).

Recently, Tsujimoto & Shigeyama (2014), Ishimaru et al. (2015) and Komiya & Shigeyama (2016) have shown that Eu-enhanced halo stars are expected in a Galactic halo that is formed from merging subhaloes, where the Eu originates from NSM with a coalescence time of $\sim$100 Myr. There is only a small probability that the subhaloes with small galaxies ($\sim 10^{5}M_{\odot}$) host a NSM at early times, but if this occurs the next generation of stars formed in such a galaxy would be extremely enhanced in Eu (eg. [Eu/Mg]$\gtrsim+2$), because of the large Eu yields per NSM. Therefore, this can be considered an other argument in favor of the NSM production scenario for r-process elements, as it is now generally accepted that the Galactic halo has in part grown from the mergers of smaller systems, building blocks (BB) as it were, a scenario that is based on the hierarchical growth of structures as explored in e.g. White & Rees (1978). Observational evidence supporting growth via accretion and mergers has been reported in e.g. Ibata et al. (1994), Helmi et al. (1999) and Belokurov et al. (2006) amongst others.

However, with the exception of the work of Beniamini et al. (2016a); Beniamini et al. (2018), and Hotokezaka et al. (2018), little attention has been given in the literature to the fact that neutron star binaries receive a kick velocity at the formation time of the two neutron stars, which might expel the system out of its host galaxy at early times, when the escape velocity of the host galaxy is still small. Therefore, we investigate in this paper if it is still possible to enrich the early Milky Way galaxy in r-process elements through the NSM scenario if neutron star kicks are taken into account. We use a semi-analytic galaxy formation model different from those of Tsujimoto & Shigeyama (2014) and Ishimaru et al. (2015).

This paper is structured as follows: we explain our methods in section 2, which contains one subsection summarizing the semi-analytic techniques that are used in our galaxy formation model and another subsection explaining how the r-process (Eu) enrichment of the galaxy by NSM is modelled. We show the results for a standard model in section 3, expand upon this result and discuss variations of the standard model in section 4. In section 5, we compare the effect of the assumption that the NS receive a kick at birth between BB galaxies and surviving satellites. We end the paper with a summary and discussion section (section 6).

2 Methods

The build-up of the Galactic halo through merging building BBs is modelled using the Munich-Groningen semi-analytic galaxy formation model (Kauffmann et al., 1999; Springel et al., 2001; De Lucia et al., 2004; Croton et al., 2006; De Lucia & Blaizot, 2007; De Lucia & Helmi, 2008; Li et al., 2010; Starkenburg et al., 2013a). As a backbone of this model, three of the six (A, B and C) high-resolution Milky Way mass Aquarius dark matter halo simulations are used (Springel et al., 2008). Resolution level 2 is used ($\sim 10^{4}M_{\odot}$ per particle), because only Aquarius A is available at higher resolution. The semi-analytic model does not trace the evolution of individual elements explicitly, and assumes an instantaneous recycling approximation (IRA). Therefore, it does not accurately predict the abundance of elements that are produced on long timescales, such as Fe. Also, as we explain in section 2.2, the abundance of Eu as a function of cosmic time was put into the model explicitly for this study.

2.1 The flow of baryons in the semi-analytic model

As described in detail in the literature (eg. Starkenburg et al., 2013a), the semi-analytic galaxy formation model that we use distinguishes between three different types of galaxies. Main galaxies in the centre of a main dark matter halo (type 0), satellite galaxies in the centre of a dark matter subhalo (type 1) and so called orphan galaxies that have lost their dark matter subhalo through tidal disruption (type 2). For type 0 galaxies, the flow of baryons in the model is visualized in Figure 2. In the next paragraphs, we will briefly summarize the physical processes that drive the mass exchange between the different boxes baryonic mass can be in, since these also affect the route of Eu. As soon as a galaxy becomes a satellite, the dark filled arrows in Figure 2 are no longer allowed. One arrow with the dashed lines represents two possible heating processes, of which one is only allowed for satellites.

In the top left corner of Figure 2, the infalling (pristine) gas or outflowing hot gas is shown. The amount of infalling/outflowing gas $M_{\mathrm{i/o}}$ is determined by the difference between the assumed (average) baryonic mass $\overline{M}_{\mathrm{b}}(M_{\mathrm{vir}},z)$ of the galaxy with virial mass $M_{\mathrm{vir}}$ at redshift $z$, taking reionization into account, and the baryonic mass in the galaxy, i.e. the sum of the mass in the other five elliptically shaped boxes of Figure 2, $M_{\mathrm{b}}=M_{\mathrm{*}}+M_{\mathrm{cg}}+M_{\mathrm{hg}}+M_{\mathrm{eg}}+M_{\mathrm{bh}}$.

[TABLE]

where the sum is taken over all (i) galaxies in the halo. $\overline{M}_{\mathrm{b}}(M_{\mathrm{vir}},z)$ is used as approximated by Gnedin (2000),

[TABLE]

with the filtering mass $M_{F}(z)$ as described in Appendix B of Kravtsov et al. (2004). A cosmic baryon fraction $f_{\mathrm{b}}$ of 0.17 is assumed, consistent with the first-year WMAP results (Spergel et al., 2003), and reionization is assumed to take place between redshift 15 and 11.5 (see also Li et al., 2010). Since only hydrogen is affected by reionization, Eu does not follow this arrow. However, the process is shown for completeness.

Hot gas cools following the prescriptions described by De Lucia et al. (2004), Croton et al. (2006) and Li et al. (2010). The cooling rate is strongly dependent on the gas temperature and on its metallicity. The collisional ionization cooling curves of Sutherland & Dopita (1993) are used to model these dependencies. Galaxies with a virial temperature $T_{\mathrm{vir}}(z)$ below the atomic hydrogen cooling limit ($T\sim 10^{4}$K) are prevented to cool gas in our model, since it is assumed to be very inefficient at early times. We assume that the amount of Eu to the total metal budget of a galaxy that accounts for cooling is negligible at all times. All metals (including Eu) that are contained in the fraction of the hot gas that is cooling are also cooled onto the cold disk of the galaxy.

From the cold gas box, baryons can flow to three different other boxes of Figure 2 (or four if it is a main galaxy). For main galaxies (type 0), the energy injected by SN, as described by Kauffmann & Haehnelt (2000) is deposited into the ejected gas box of Figure 2.

[TABLE]

with $\Delta M_{\mathrm{*}}$ is the stellar mass formed, $V_{\mathrm{vir}}$ the circular velocity of the galaxy at the viral radius ($R_{\mathrm{vir}}$), $E_{\mathrm{SN}}=10^{51}$ erg the energy released per SN, and $\epsilon=0.05$ and $\eta_{\mathrm{SN}}=8\cdot 10^{-3}M_{\odot}^{-1}$ are the feedback efficiency and the number of SNe respectively (Li et al., 2010). We define the virial radius as the radius of a sphere with mass $M_{\mathrm{vir}}$ and mean interior density is 200 times the critical density for closure of the universe (Navarro et al., 1997):

[TABLE]

with H the Hubble parameter. Smaller galaxies have shallower potential wells and thus a lower $V_{\mathrm{vir}}$, which makes them more efficient reheaters.

Two different physical processes are behind the arrow from cold gas pointing back to hot gas (highlighted with dashes), describing reheating of the gas. One is related to the radio activity of the active galactic nucleus (AGN) of the galaxy, as described by Croton et al. (2006), and should technically be considered a suppression of the cooling flow. This recipe is applied to all galaxy types. The other heating recipe again corresponds to the energy injected by SN, equation (3), but then only applied for satellite galaxies (types 1 and 2). For those, the virial mass is taken to be the number of particles in the simulated galaxy times the particle mass. For both reheating recipes, all metals (including Eu) that are contained in the fraction of the cold gas that is reheated are moved to the hot gas phase of the host galaxy.

The baryons in the ejected gas box, $M_{\mathrm{eg}}$, can be re-incorporated into the hot gas box, following the prescription of De Lucia et al. (2004):

[TABLE]

with $\dot{M}_{\mathrm{re-inc.}}$ the amount of gas that is reincorporated during one timestep, and $\gamma_{\mathrm{eg}}=0.5$ the ejected gas reincorporation efficiency (Croton et al., 2006, table 1). The ratio $R_{\mathrm{vir}}/V_{\mathrm{vir}}$ is proportional to the dynamical time of the galaxy (Springel et al., 2001). As for all the other processes, all metals (including Eu) contained in the fraction of the ejected gas that is re-incorporated are moved to the hot gas box.

Only main galaxies have ejected gas and hot gas components. Satellites transfer their ejected mass to their hosts (the mains). This is indicated through the dashed line from Stellar Mass to Hot Gas in Figure 2.

Both cold gas and hot gas can be accreted onto the central black hole (BH). With $M_{\mathrm{bh}}$, we denote the baryonic part of the BH mass. The change in this parameter during one timestep is indicated with an overdot. From the hot gas box, quiescent accretion follows the empirical recipe described by Croton et al. (2006), i.e. the ‘radio mode’:

[TABLE]

Here, $M_{\mathrm{hg}}$ is the baryonic mass in the form of hot gas and $\kappa_{\mathrm{AGN}}=7.5\cdot 10^{-6}M_{\odot}\mathrm{yr}^{-1}$ the hot gas BH accretion rate (De Lucia & Blaizot, 2007). Cold gas accretion happens when a satellite galaxy merges with the central (main) galaxy. This recipe is also described by Croton et al. (2006) and dubbed ‘quasar mode’:

[TABLE]

with $M_{\mathrm{*}}$ and $M_{\mathrm{cg}}$ the stellar mass and mass in cold gas of the main galaxy and $m_{\mathrm{*}}$ and $m_{\mathrm{cg}}$ the stellar mass and mass in cold gas of the satellite galaxy respectively. $f_{\mathrm{bh}}=0.03$ is the cold gas BH accretion fraction. Again, all metals (including Eu) contained in the fraction of the hot or cold gas that is accreted are moved to the central BH box.

Last, but not least, we describe how stars are formed from cold gas and how metals are added to the hot and cold gas components. We assume cold gas is transformed into stars in an exponential thin disc (by the formalism of Mo et al., 1998), which is assumed to extend to three scalelengths, when its density is above a critical threshold (Kennicutt, 1989; Kauffmann, 1996):

[TABLE]

The star formation rate is then proportional to the amount of cold gas available. See also De Lucia & Helmi (2008) and Starkenburg et al. (2013a). Furthermore, during a galaxy merger, a collisional starburst recipe (Somerville et al., 2001) is triggered, where a fraction

[TABLE]

of the combined cold gas of the two galaxies ($m_{\mathrm{cg}}+M_{\mathrm{cg}}$) is turned into stars. The parameters $\alpha_{\mathrm{burst}}$ and $\beta_{\mathrm{burst}}$ are set to be 0.7 and 0.56 respectively (Croton et al., 2006). Due to our adopted IRA and assumed Chabrier (2003) initial mass function (IMF), 43% of the cold gas is returned immediately to the cold gas component during star formation, hence the two-way pointing arrow in Figure 2. Also, 43% of the metals (including Eu) in the cold gas available for star formation are returned immediately to the cold gas, the other 57% ends up in stars. After a starburst and instantaneous recycling, new metals (other than Eu, which is injected through a different channel - see section 2.2) become available in the two gas phases. We follow the recipe of Li et al. (2010), who suggested a two-state value of the fraction of metals deposited directly into the hot gas of the main galaxy in the dark matter halo (which is the galaxy itself in case it is a type 0). For galaxies with a virial mass below $5\cdot 10^{10}M_{\odot}$, 95% of the metals are added to this box and 5% to the cold gas box of the galaxy that produced the metals. Higher mass galaxies deposit all their new metals into the cold gas phase of their galaxy. We assume a metal yield of 0.03 (Croton et al., 2006).

2.2 Adding the europium

First we run the semi-analytic model once without any enrichment in europium, and record the total stellar mass formed in each galaxy, as well as the maximum circular velocity of a particle in the potential well of the galaxy at each time step. Then, assuming a NSM rate, we know in post-processing for each star formation episode how many NS mergers happen after a delay time picked out of an assumed delay time distribution (DTD). In our standard model, we assume a NSM rate of $10^{-4}$ per $M_{\odot}$ stars formed. This assumption is based on the NS merger rate inferred from gravitational wave (GW) observations (Côté et al., 2018). It is about a factor $10$ higher than the rate assumed by Tsujimoto & Shigeyama (2014), who deduce a NSM rate of one per $\sim 1000$ core-collapse SNe111In our semi-analytic model we assume $8\cdot 10^{-3}$ core-collapse SNe per $M_{\odot}$ stars formed (Li et al., 2010)., and about a factor $100$ higher than that of Cescutti et al. (2015), who assume a NSM rate of a few times $10^{-6}$ per $M_{\odot}$ stars formed.

In our standard model, we assume a DTD proportional to $1/t$, with the first NSM occurring after 10 Myr, which we will abbreviate as 10 Myr $+1/t$ in the remainder of this paper. As already mentioned in the introduction, this DTD is compatible with binary population synthesis studies (Belczynski et al., 2006).

Assuming a NFW profile for the density of dark matter haloes at radius $r$ (Navarro et al., 1996),

[TABLE]

where $\rho_{0}$ and $a$ are fit parameters, it can be shown that the escape velocity that a particle needs to have in order to escape from the centre of the potential well is

[TABLE]

where $V_{\mathrm{max}}$ is the maximum of the circular velocity of the galaxy. At the scale radius $a$, the escape velocity is a factor $\sqrt{\ln(2)}=0.833$ lower.

We assign a kick velocity to each NS binary at the formation time of the two NS, following the kick velocity distribution of Arzoumanian et al. (2002): a two-component velocity distribution (Maxwellians) with characteristic velocities of 90 and 500 km/s, with 40% of the population being drawn from the distribution with the lowest characteristic velocity. Since we assume that the system remains bound after receiving the kick, we assign half the velocity picked from this distribution as the kick velocity of the binary. We compare this velocity with the escape velocity from the centre of the galaxy hosting the NS binary (equation 11). If the kick velocity is smaller than the escape velocity, we assume that the neutron star binary will be inside the galaxy at the time of the merger. Otherwise, the binary is discarded as a potential source of Eu enrichment. We also investigate the effect of applying the low velocity Maxwellian distribution of Arzoumanian et al. (2002) only, i.e. that with a characteristic velocity of 90 km/s. This assumption is hereafter referred to as the “low kick” model, which is compared to the “no kick” (standard kick model) and “with kick” (assuming the full kick velocity distribution) models. The choice for the low kick model is motivated by the fact that many double neutron stars are found to have small NS kick velocities (Tauris et al., 2017). Furthermore, Verbunt & Cator (2017) have recently shown that a description with two Maxwellians, with distribution parameters $\sigma_{1}=77$ and $\sigma_{2}=320$ km/s, is significantly better than that of Arzoumanian et al. (2002).

For those binaries that do not leave the host galaxy after receiving their kick, we assume a europium yield of $1.5\cdot 10^{-5}M_{\odot}$/NSM in our standard model. This is the upper limit of the Eu yield inferred from GW observations (Côté et al., 2018) and slightly larger than what was assumed by Matteucci et al. (2014, 2015), who assume Eu yields in the range $10^{-5}-10^{-7}M_{\odot}$, but a factor $10$ smaller than the estimate of Komiya & Shigeyama (2016). We assume that halo stars are only formed in BBs. If a NSM happens when its hosting BB has already merged with the central galaxy, the gas enriched by this NSM might form new stars in the disk or inner spheroid, but since this paper focusses on the halo component where no new star formation takes place this is not relevant for our research. Since the amount of Eu is negligible compared to the total NSM metal yield, which is already added to its metal budget through the “metals from stars to gas” channel (see section 2.1), we also do not add the Eu to the gas metal budget of the main galaxy’s spheroid in this scenario.

Knowing exactly when a NSM happens in each galaxy, we now re-run the semi-analytic model and incorporate the chemical evolution of Eu in the ensemble of galaxies. The europium will be instantaneously mixed into the gas of each galaxy, which is used for the formation of new stars. Similar to the enrichment of the galaxy in other metals (see section 2.1), we assume that 95% of the newly produced Eu is mixed with the hot gas in the main galaxy in the dark matter halo if the galaxy that produced the Eu has a virial mass below $5\cdot 10^{10}M_{\odot}$, and 5% is added to the cold gas budget of their own galaxy222In case the galaxy is not a satellite, the main galaxy in the dark matter halo is also the galaxy itself.. All Eu is added to the cold gas budget of the galaxy itself in case it has a virial mass above $5\cdot 10^{10}M_{\odot}$.

2.3 Metallicities

The total metallicity of a stellar population at each time step is traced by our model as the parameter $Z_{\mathrm{stars}}$, the ratio of mass in metals over the total mass in stars. Because of the IRA, this value can thus best be compared with the abundance of elements which are mainly originating in (short-lived) Supernovae type II, such as the $\alpha$-elements. In this paper we show predicted $\log[Z_{\mathrm{stars}}/Z_{\odot}]$ values of our model, with $Z_{\odot}=0.0142$ the reference metallicity (Asplund et al., 2009). These can thus be thought of as [Mg/H] values. Similarly, $\log[(X_{\mathrm{Eu}}/Z)_{\mathrm{stars}}/(X_{\mathrm{Eu}}/Z)_{\odot}]$ values, with $X_{\mathrm{Eu},\odot}=3.6\cdot 10^{-10}$ (Asplund et al., 2009), can be thought of as [Eu/Mg]. That is why we chose to plot these abundance ratios against each other in Figure 1, instead of the more commonly used [Fe/H] and [Eu/Fe].

3 Results for a standard model

3.1 The accreted spheroid

In the upper left-hand panel of Figure 3, we show the [Eu/Mg] versus [Mg/H] map predicted by our model for the accreted spheroid of Aquarius halo A-2, assuming that the NS do not receive a kick velocity at their formation time (see section 2.2). Colors indicate the stellar mass per bin, on a logarithmic scale, with darker colors representing more stars.

The accreted spheroid contains bulge stars and halo stars. Their combined metallicity distribution function (MDF), is shown with a green line in the lower left-hand panel of Figure 3. Since the stars selected from the VizieR database and presented in Figure 1 that we want to compare our model [Eu/Mg] vs. [Mg/H] maps with are all halo stars, we want to divide this accreted spheroid into a bulge part and a halo part. To do this, we use the same technique as van Oirschot et al. (2017b). The semi-analytic model underestimates the number of accreted spheroid stars with metallicities $-3\leq\log[Z_{\mathrm{stars}}/Z_{\odot}]<-1$. Therefore, we assume all stars in these and lower metallicity bins to be halo stars, and consequently all accreted bulge stars to have metallicities $\log[Z_{\mathrm{stars}}/Z_{\odot}]>-1$. The fraction of halo stars in these $\log[Z_{\mathrm{stars}}/Z_{\odot}]>-1$ bins of the accreted spheroid MDF is set to follow a single Gaussian fit to the observationally determined MDF by An et al. (2013), which is normalized to the $-1.5\leq\log[Z_{\mathrm{stars}}/Z_{\odot}]<-1.0$ bin of the accreted spheroid MDF and shown with a blue line in the lower left-hand panel of Figure 3. The red line is the resulting stellar halo MDF of the semi-analytic model (SAM).

The downscaling of the stellar mass in the $\log[Z_{\mathrm{stars}}/Z_{\odot}]>-1$ bins of the [Eu/Mg] vs. [Mg/H] map of accreted spheroid Aq-A-2 is done linearly for all age bins, and the resulting map is shown in the right-hand panel of Figure 3. The same scale for the colormap is used as for the [Eu/Mg] vs. [Mg/H] map in the left panel. The observed halo stars shown in Figure 1 are overplotted with blue circles here (i.e. stars known to be in satellite galaxies are left out). Open circles represent stars for which only upper limits are available, filled circles represent stars which have standard errors in the determined [Eu/Mg] values.

It is a known shortcoming of the Munich-Groningen semi-analytic model that it underestimates the number of low metallicity halo stars. Therefore, it is not surprising that our model cannot explain the large number of observed stars with [Mg/H] < -1.5. More importantly, hardly any of the few low-metallicity stars that are predicted by our model, have [Eu/Mg] > 0 abundances, whereas these are observed in the Galactic halo.

In the next section we give an example of how the Galactic halo is built up from BBs, in order to investigate if we can confirm the findings of Tsujimoto & Shigeyama (2014), Ishimaru et al. (2015) and Komiya & Shigeyama (2016) that Eu-enhanced halo stars naturally occur in a Galactic halo that is formed from merging subhaloes. In the remainder of this paper we will not select out the accreted bulge stars from modelled spheroids [Eu/Mg] vs. [Mg/H] maps, because we will focus our comparisons with data solely on the lower metallicities where we expect to have accreted halo stars only. However, the same method could be applied as in Figure 3 to end up with [Eu/Mg] vs. [Mg/H] maps of accreted halo stars only.

3.2 Building blocks of the Galactic spheroid

In this paper, we consider the accreted component of the stellar spheroid, which is assumed to be built from a few main progenitor galaxies, in agreement with the findings of many others (eg. Helmi et al., 2002, 2003; Font et al., 2006a; Cooper et al., 2010; Gómez et al., 2013). In order to study the contribution of one of these large BBs of the stellar spheroid to the [Eu/Mg] vs. [Mg/H] map in depth, we show a part of the galaxy merger tree of Aquarius halo A-2 in Figure 4 (bottom right). Only the main (simulated Milky Way) galaxy and one major BB are shown, respectively with numbers 0 and 1. BB 1 itself is built from four other BBs (labelled 2, 3, 5 and 6), of which some in turn are also built from building blocks.

In Figure 4 we also show the [Eu/Mg] vs. [Mg/H] maps of these BB galaxies as a function of cosmic (lookback) time. Each panel contains a [Eu/Mg] vs. [Mg/H] map, with successive time snapshots of the same BB vertically ordered. The labels of the axes are left out for clarity, but they are identical to those of Figure 3. However, the chosen colormap which again indicates the stellar mass in each bin of a [Eu/Mg] vs. [Mg/H] map, is slightly different from the one used in Figure 3. This is because Figure 4 also contains the stars in the disc of the main (Milky Way) galaxy of Aquarius A-2 in its leftmost panels (hence this colormap extends to $10^{11}M_{\odot}$), whereas Figure 3 only contains accreted spheroid stars (and its colormap extends only to $10^{9.5}M_{\odot}$). The structure of the galaxy merger tree (bottom right) is matched in the main part of the figure. The [Eu/Mg] vs. [Mg/H] maps of BBs that have not yet formed any stars, or those of the BBs that merged with other BBs in a particular time step, are not shown.

From this figure, it is clear that all except the lowest mass BB of this major BB (i.e. number 2, which has a maximum stellar mass of $10^{3.5}M_{\odot}$) form stars with Eu. This is a result of the assumed NSM rate of $10^{-4}$ per $M_{\odot}$ stars formed, since that makes the threshold BB mass for forming stars which are enriched in Eu $\sim 10^{4}M_{\odot}$. Furthermore, we see from this figure that most of the few accreted spheroid stars with enhanced Eu abundances in this model (shown above the dashed line in each panel), as well as those with low-metallicities ($\log[Z_{\mathrm{stars}}/Z_{\odot}]<1.5$) were not born in (BBs of) this major BB.

Then of course the question arises, how can we explain the large number of halo stars with [Eu/Mg] > 0 abundances? We will investigate the effect of different assumptions on the NS kick, the NSM rate or the DTD in the next section, but it might very well be our assumption of instantaneous mixing that results in these low [Eu/Mg] abundances. Eu-enhanced stars naturally arise due to the inhomogeneity of gas enriched by NSM, as was also found on the subgrid level of the cosmological magnetohydrodynamical IllustrisTNG simulations recently (Naiman et al., 2018).

4 The effect of NS kicks, the NSM rate and the DTD

Although our standard model does predict some halo stars with high Eu-abundances for stars with low [Mg/H] values, shown above the dashed line in the [Eu/Mg] vs. [Mg/H] maps presented in the previous section, compared with observations (right-hand panel of Figure 3) the number density of Eu-enhanced stars predicted by the SAM applied to Aquarius halo A-2 is too low. Therefore, we will investigate variations of our standard model in this section. Hereafter, we define enrichment in Eu to take place if all of the following three conditions are met:

• •

A NSM takes place in the BB.

• •

The coalescence time of the NS binary is shorter than the timescale on which the BB merges with the central galaxy.

• •

The NS kick velocity is smaller than the escape velocity of the BB.

The assumed NSM rate affects the first condition, the DTD affects the second condition and the assumed binary NS kick velocity distribution affects the third condition. We will show the effect of variations in all of these three variables.

4.1 NS kick velocities

In our standard model, we assume that NS binaries do not receive a kick, i.e. the third of the above conditions is always met. Therefore, introducing a kick velocity naturally decreases the fraction of Eu-enriched halo stars. To investigate by how much, we show in Figure 5 the fraction of BBs with Eu-enrichment by NSM as a function of BB stellar mass, for three different assumptions on the NS kick velocity. We also show the effect of two different assumptions on the DTD for each of these three models.

With our standard model assumptions (no kick) $\sim 90$% of the highest mass BBs ($m_{*}>10^{8.5}M_{\odot}$) and $\sim 80$% of the lowest mass BBs ($m_{*}<10^{5.5}M_{\odot}$) is enriched in Eu. Assuming longer delay times through a 100 Myr $+1/t$ DTD in this case give $5-10$% lower enrichment fractions, because fewer NS binaries merge before the BB merges with the central galaxy in that model. Introducing the 2-peak Arzoumanian et al. (2002) kick velocity distribution, drastically decreases this fraction, to less than 50% of the highest mass BBs and to only 10% of the lowest mass BBs. This result is almost identical if we take a 100 Myr $+1/t$ DTD instead of our standard 10 Myr $+1/t$. The enrichment fractions corresponding to the low kick models (only the low velocity peak of the Arzoumanian et al. (2002) distribution) fall in between the other lines for low mass BBs, and overlap with the lines for no NS kick for the highest mass BBs because these highest mass BB have escape velocities which are always larger than these low NS kick velocities.

In Figure 6, the [Eu/Mg] vs. [Mg/H] maps are shown for the full accreted spheroids of Aq-A-2, Aq-B-2 and Aq-C-2 from left to right respectively. The top, middle and bottom panels show the difference between assuming that the NS binary receives no kick, a low kick, or a higher kick. A DTD of 10 Myr $+1/t$ is assumed. Figure 6 shows that there is indeed a larger Eu-enrichment in case the NS do not receive a kick. Roughly, the distribution of stars in these [Eu/Mg] vs. [Mg/H] maps shifts 0.5 dex downwards when we assign a NS kick velocity. Furthermore, this figure gives an impression of the stochastic scatter between the Aquarius haloes. One cannot expect to find a galaxy in one of the Aquarius haloes that is an exact replica of our Milky Way, because only six simulations were done (of which we study three in this paper). In particular, our semi-analytic model “no kick” and “with low kick” applied to Aquarius halo B results in an [Eu/Mg] vs. [Mg/H] map that is in very good agreement with the observations. The majority of r-process enriched stars in this halo originate in a single relatively massive ($10^{7.45}M_{\odot}$) BB, in which a relatively large starburst occurred at late times ($10$% of its stars were formed in a burst about $9.5$ Gyr ago, whereas most halo stars were formed about $11$ Gyr ago). Our model thus predicts quite some variation between galaxies in terms of r-process enrichment.

4.2 Delay time distribution

The effect of assuming a different DTD on the [Eu/Mg] vs. [Mg/H] map of stellar spheroid Aq-A-2 is shown in Figure 7. A constant delay time of 1 Myr may be unrealistically short to be compatible with population synthesis studies (Belczynski et al., 2006), but this is what has been assumed before to explain Eu abundances in the Milky Way galaxy with a NSM enrichment scenario (Matteucci et al., 2014). The models in which a fixed delay time of 1 Myr is assumed yield [Eu/Mg] vs. [Mg/H] maps where most stars have [Eu/Mg]$\sim 0$ with little scatter, since NSM are treated almost like SNe Type II with this assumed DTD and fairly high NSM rate. The bulk of the stars is shifted to higher [Eu/Mg] abundances by about $\sim$0.5 dex compared to our standard model (10 Myr $+1/t$). Thus the effect is similar to that of assuming a kick velocity for the binary NS.

Assuming a longer DTD on the other hand (100 Myr $+1/t$), results in the majority of stars being spread out over the low [Eu/Mg], high [Mg/H] abundance bins (eg. $0\gtrsim$[Eu/Mg]$\gtrsim-2$ for $0\gtrsim$[Mg/H]$\gtrsim-1$), which is clearly not where the majority of the observed halo stars are (see Figures 1 and 3).

4.3 NSM rate

Three different NSM rates are assumed in Figure 8. These are combined with Eu yield assumptions such that the Galactic disc stars have solar Eu abundances at [Mg/H]=0 in our models, i.e. the production of the NSM rate and the Eu yield always equals $1.5\cdot 10^{-9}$/NSM. With high NSM rates ($10^{-3}M_{\odot}^{-1}$) and low Eu yields ($1.5\cdot 10^{-6}M_{\odot}$/NSM), the effect of continuous Eu production is studied, similar to how it would be in a supernova enrichment scenario. We also investigate a scenario in which NSM are a factor 10 rarer than in our standard model ($10^{-5}M_{\odot}^{-1}$), but yield a large amount of Eu ($1.5\cdot 10^{-4}M_{\odot}$/NSM).

In the left-hand panels of Figure 8, we show the models with an increased Eu yield and correspondingly lower NSM rate, compared to our standard model. Here, only the most massive BBs are expected to form Eu-enriched stars. This can also be seen from the number of stars without any Eu, that end up in the lowest [Eu/Mg] bin ([Eu/Mg]$<-3.5$, below the solid line) of each [Eu/Mg] vs. [Mg/H] map. There are much more of these in the left-hand panels of Figure 8 (low NSM rate) than in the right-hand panels (high NSM rate). The difference between the middle and the right-hand panels seems not very large.

5 Comparing BBs with surviving satellites

Defining galaxies that contain stars with [Eu/Mg]>0 to be Eu-rich, and those with some Eu-enriched stars but no stars above the solar [Eu/Mg] abundance as Eu-poor, we show the average number of galaxies in these two categories as a function of galaxy stellar mass in Figure 9, for the three different assumptions on the NS kick. The average is taken over the three Aquarius haloes A, B and C. We distinguish BBs (filled gold histograms) and surviving satellites (transparent histograms with thick blue edges) to investigate the effect of the timescale on which the galaxies merge with the central galaxy: surviving satellites per definition have not (yet) merged with the central galaxy in the simulation. If a galaxy is not fully disrupted, it is counted as partially stripped and divided beween the BB and surviving satellite categories according to the fraction of the total stellar mass that is stripped. For example, if 70% of its initial mass is stripped, the galaxy is counted as 0.7 BB stripped material and 0.3 surviving satellite galaxy.

We do not trust our model’s prediction for galaxies with a stellar mass below a few times $10^{4}M_{\odot}$ due to the limiting resolution of our simulation. Largely overlapping with these, the galaxies with less than five star formation snapshots are left out of Figure 9, because the neglect of pre-enrichment from the very first generation of stars in our model is most clear in these galaxies (Starkenburg et al., 2013a). This can be thought of as a mass cut since more massive galaxies sustain star formation for more extended periods and show star formation in more than five snaphots. The number five is a consequence of the time resolution of our code, which is redshift dependent but in the order of tens to hundreds of Myr. The ultra-faint dwarf galaxy Reticulum II (Ret II), which was recently discovered to be extremely enriched in r-process elements (Ji et al., 2016a, b), falls below this mass resolution limit (Simon et al., 2015). Although we assume a NSM rate of only one per $10^{4}M_{\odot}$ stars formed, due to Poisson scatter galaxies with stellar masses as low as that of Ret II can have Eu-rich stars in the NSM r-process enrichment scenario. In the left-hand panels of Figure 9, it can be seen that this is most likely to be the case if the NS binary receives no kick, or a kick with a low velocity.

The fact that there is a large difference between the number of stars in the three categories between the three assumptions on the NS kick, is related to our definition of Eu-rich stars, i.e. those with [Eu/Mg]>0. By introducing a kick velocity, the distribution shifts to lower [Eu/Mg] abundances, by about $\sim$0.5 dex. Those galaxies with [Eu/Mg] abundances just above zero move to values just below zero when a kick is introduced, which makes them fall into a different category.

Figure 9 also predicts a difference between the Eu enrichment of the stellar halo of the Milky Way and that of its surviving satellite galaxies. Surviving satellites more often are expected to have Eu-rich stars than those that are fully disrupted. An explanation for this can be found in both the time the galaxy takes to merge with the central galaxy, and in the average timespan of star formation, that is longer for surviving satellite galaxies than for disrupted BBs (van Oirschot et al., 2017a). Galaxies with a longer star formation timespan for a similar final stellar mass have a lower star formation rate, a slower build-up of the metallicity (Mg) budget due to Type-II supernovae, and correspondingly more stars with high [Eu/Mg] abundance. A similar argumentation is used to explain the high alpha-over-iron element abundances in the stellar halo compared to that of surviving satellite galaxies (eg. Font et al., 2006a, b; Geisler et al., 2007; Robertson et al., 2005).

In the top left panel of Figure 9 we see that in our standard model, 27% of the modelled Milky Way spheroid stars with [Eu/Mg]>0 comes from BBs that still have a surviving counterpart. For the models “with kick” this ratio is as high as 54%, as can be seen in the bottom left panel of Figure 9. However, it is unlikely that many stars with high Eu abundances in the Milky Way halo are stripped from surviving satellites, since the Eu-enhanced are local field halo stars for which no known stream or parent is known. Furthermore, these stars are [$\alpha$/Fe] enhanced, which the satellite galaxies that have survived are not.

6 Summary and discussion

In this paper, we have modelled the europium production by neutron star mergers in the halo of the Milky Way with the Munich-Groningen semi-analytic galaxy formation model. In particular, we have investigated the effect of the kick velocity the NS binary receives upon its formation, which may lead to a NS merger outside the BB galaxy hosting the double NS. Although this may lead to Eu enrichment in the inner spheroid or disk of the main galaxy in the simulation, in our models this situation will not lead to enrichment of the stellar halo, which is assumed to be formed from merging BBs only. The effect is dependent on the NS kick velocity and the escape velocity of the hosting BB. The NS kick velocity influences the amount of enrichment of the Galactic spheroid in Eu: in [Eu/Mg] vs. [Mg/H] stellar density maps, stars have $\sim 0.5$ dex lower [Eu/Mg] abundances if a NS kick drawn from the Arzoumanian et al. (2002) velocity distribution is assumed (Figure 6). Of the BBs with a stellar mass below $10^{5.5}M_{\odot}$, only 10% will be enriched in Eu when a NS kick is assumed, compared to $\sim 80$% when no NS kick is assumed (Figure 5). For more massive BBs the enrichment is always larger, as these have deeper potential wells from which NS binaries have to escape.

We assumed the velocity of the center of mass of the binary to be half the kick velocity, which is slightly oversimplifying the effect of the binary evolution. Neglecting the effect of the kick after the supernova of the first star, and assuming the star that undergoes the second supernova drags along the other star, which is assumed to be of equal mass, we arrive at the factor 2. Once a complete binary stellar evolution model is combined with the semi-analytical galaxy formation model, a more accurate value for the system velocity after the supernova of the second star can be assumed.

Our simulations indicate that our standard assumption on the DTD (proportional to $1/t$, with the first merger occurring after 10 Myr) results in [Eu/Mg] vs. [Mg/H] maps which match that of the observed stellar halo population reasonably well, for stars with metallicities [Mg/H]$\geq-1.5$. Although we confirm the result of Tsujimoto & Shigeyama (2014), Ishimaru et al. (2015) and Komiya & Shigeyama (2016) that a stellar halo built from merging building block galaxies naturally hosts stars with high Eu abundances at low metallicities, our model predicts too few low-metallicity spheroid stars. The MDF bins with metallicities between $-3\leq\log[Z_{\mathrm{stars}}/Z_{\odot}]<-1$ contain too few stars compared to the single-Gaussian fit to the observed MDF of the Galactic halo by An et al. (2013). It is exactly in this metallicity range, where surprisingly many stars are found to be enriched in Eu (Frebel, 2010). Contrary to the findings of others, eg. Komiya & Shigeyama (2016) or Naiman et al. (2018), we do not find many Eu-rich stars (i.e. those with [Eu/Mg]>0) at low metallicity in our model of the Galactic stellar halo, although for Aquarius halo B there are stars with [Eu/Mg]>0.5 at $-3\leq\log[Z_{\mathrm{stars}}/Z_{\odot}]<-1.5$ in our simulations when we assume that NS binaries do not receive a kick velocity or receive a low kick velocity.

A delay time proportional to $1/t$, but with the first merger occurring after 100 Myr (instead of 10 Myr as assumed in our standard model) results in too many stars with high metallicity and low Eu abundance (Figure 7). A shorter, fixed delay time of 1 Myr, which Matteucci et al. (2014) concluded to be necessary for a pure NSM-enrichment scenario for the Galactic stellar halo, results in [Eu/Mg] vs. [Mg/H] maps with little scatter around [Eu/Mg]$\sim 0$. Furthermore, population synthesis studies favor longer delay times (Belczynski et al., 2006). Reducing or neglecting the NS kick velocity has a similar effect as assuming a shorter delay time, with respect to increasing the Eu-abundances of the modelled halo stars. To enhance the r-process enrichment of a simulated galaxy, one could either increase the Eu yield per NSM or assume a higher NSM rate. In Figure 8 we resolve this degeneracy in high NSM rate and low r-process yield versus low rate and high yield. We find that assuming a lower NSM rate and higher Eu yield compared to our standard model results in a [Eu/Mg] vs. [Mg/H] map that is less in agreement with observations of Eu-enhanced halo stars. A higher NSM rate and lower Eu yield on the other hand, does not solve the issue that there are too few Eu-enhanced halo stars with $\log[Z_{\mathrm{stars}}/Z_{\odot}]([Mg/H])<-1.5$ in our standard model.

Although we do not assume instantaneous recycling of Eu, our neglect of the short ($\sim 10$ Myr) but nonzero recycling time of alpha-elements such as Mg could partly explain the discrepancy with the observed stellar halo stars. It is more likely however, that the absence of these halo stars with high Eu abundances at low metallicities indicates that our assumption of instantaneous (i.e. homogeneous) mixing is incorrect. From the chemical signatures of the metal-poor stars in dwarf galaxies there have been hints of inhomogeneous mixing processes, for instance of SNe Type Ia products (see, for example Venn et al., 2012; Starkenburg et al., 2013b; Jablonka et al., 2015), although the energies and outflow of new material of SN Ia are quite different from the energies and outflow of NSM products. This paper provides another line of evidence that the gas in dwarf galaxies and BBs is probably inhomogeneously mixed at early times. Furthermore, there are several works that suggest that the IMF in ultra-faint dwarf galaxies can be very different from the standard one (e.g. Geha et al., 2013). Since a very different IMF can change the neutron star density in such a stellar population, this leads to even more inhomogeneity in r-process enrichment.

An alternative conclusion one could draw, is that an additional production mechanism for r-process elements is necessary to explain the discrepancy between our models and the observations. This mechanism could enhance halo stars in r-process elements already at early times, when the metallicity of the halo stars is low and many building blocks of the stellar halo still have shallow potential wells from which NS binaries easily escape after receiving a standard kick velocity. This scenario is supported by the recent observation of a only mildly r-process enhanced star in the ultra-faint dwarf galaxy Tucana III (Hansen et al., 2017), contrary to the highly r-process enhanced stars observed in Ret II.

Finally, the star formation history of (the BBs of) the stellar halo could play a role, since we found that Aquarius halo B is more in line with observations than haloes A and C (see Figure 6, which gives an impression of the stochastic scatter between the Aquarius haloes). In a comparison between the Eu-enrichment in our modelled BB galaxies of the Galactic spheroid and the surviving satellite population, we also find that the surviving satellites on average more often have stars with [Eu/Mg]>0 abundances than the fully disrupted BBs and that low metallicity spheroid stars with high Eu abundances are often stripped from a satellite galaxy with a still surviving counterpart (Figure 9). However, this is unlikely to be the case in the Milky Way halo, since there is no observational evidence supporting it.

van de Voort et al. (2015) found that stars with low [Fe/H] and high [Eu/Fe] are at large galactocentric radius in a cosmological zoom-in simulations of a Milky Way-mass galaxy, but they are formed in-situ at high redshift and moved there through dynamical effects. Furthermore, Naiman et al. (2018) find no correlation between europium abundances and assembly history, in Milky Way-like galaxies in the IllustrisTNG simulations.

Acknowledgements

The authors are indebted to the Virgo Consortium, which was responsible for designing and running the halo simulations of the Aquarius Project and the L-Galaxies team for the development and maintenance of the semi-analytical code. In particular, we are grateful to Gabriella De Lucia and Yang-Shyang Li for the numerous contributions in the development of the code. We also thank Amina Helmi for valuable comments.