Neutron star binary orbits in their host potential: effect on early r-process enrichment

TL;DR

This study investigates how neutron star binary orbits within their host galaxy's gravitational potential influence the distribution and enrichment of r-process elements, revealing challenges to the merger scenario based on observed elemental abundances.

Contribution

It provides a detailed analysis of neutron star binary dynamics in different galaxy types and their impact on r-process enrichment, highlighting the limitations of the merger scenario in explaining observed abundances.

Findings

Neutron star binaries tend to merge outside shallow dwarf galaxies, causing dilution of r-process elements.

There is a large scatter in r-process enrichment predictions due to geometrical dilution and small number statistics.

Observed europium abundances challenge the merger scenario, suggesting additional r-process sources or faster mergers.

Abstract

Coalescing neutron star binary (NSB) systems are primary candidates for $r$-process enrichment of galaxies. The recent detection of $r$-process elements in ultra-faint dwarf (UFD) galaxies and the abundances measured in classical dwarfs challenges the NSB merger scenario both in terms of coalescence time scales and merger locations. In this paper, we focus on the dynamics of NSBs in the gravitational potentials of different types of host galaxies and on its impact on the subsequent galactic enrichment. We find that, for a $\sim t^{-1}$ delay time distribution, even when receiving a low kick ($\sim 10~{\rm km~s^{-1}}$) from the second supernova explosion, in shallow dwarf galaxy potentials NSBs tend to merge with a large off-set from the host galaxy. This results in a significant geometrical dilution of the amount of produced $r-$process elements that fall back and pollute the host galaxy gas reservoir. The combination of dilution and small number statistics produces a large scatter in the expected $r$-process enrichment within a single UFD or classical dwarf galaxy. Comparison between our results and observed europium abundances reveals a tension that even a systematic choice of optimistic parameters in our models cannot release. Such a discrepancy could point to the need of additional $r$-process production sites that suffer less severe dilution or to a population of extremely fast merging binaries.

Tables

Table 1. Table 1: Summary of the most relevant properties of the different dwarf galaxy models employed in our study (full details about MW-like parameters can be found in Barros et al. ( 2016 ) , tables 3-4). M b subscript 𝑀 b M_{\rm b} and M DM subscript 𝑀 DM M_{\rm DM} are the baryonic and dark matter masses, T SF subscript 𝑇 SF T_{\rm SF} is the duration of the star formation, M ∗ ​ ( T SF ) subscript 𝑀 subscript 𝑇 SF M_{*}(T_{\rm SF}) is the mass converted in star at T SF subscript 𝑇 SF T_{\rm SF} , R d subscript 𝑅 d R_{\rm d} and z d subscript 𝑧 d z_{\rm d} the scale radius and height of the disk, N CCSN ​ ( T SF ) subscript 𝑁 CCSN subscript 𝑇 SF N_{\rm CCSN}(T_{\rm SF}) the average number of CCSN exploded by T SF ) T_{\rm SF}) , while the two values of N NSB subscript 𝑁 NSB N_{\rm NSB} represent the minimum and the maximum average number of NSB formed. The least (most) optimistic case is obtained considering one NSB every 50 (1350) CCSNe.

Model	$M_b$	$M_{\mathrm{DM}}$	$T_{\mathrm{SF}}$	$M_{\ast }\left(T_{\mathrm{SF}}\right)$	$R_d$	$z_d$	$N_{\mathrm{CCSN}}\left(T_{\mathrm{SF}}\right)$	$N_{\mathrm{NSB}}\left(T_{\mathrm{SF}}\right)$
	$\left[M_{\odot }\right]$	$\left[M_{\odot }\right]$	[Gyr]	$\left[M_{\odot }\right]$	$\left[\mathrm{kpc}\right]$	$\left[\mathrm{kpc}\right]$		min	max
Dwarf 1	${10}^5$	${10}^7$	1.0	$6.0\times {10}^3$	0.07	0.014	56	0.04	1.12
Dwarf 2	${10}^6$	${10}^8$	1.6	$9.5\times {10}^4$	0.15	0.030	886	0.65	17.7
Dwarf 3	${10}^7$	${10}^9$	2.5	$1.5\times {10}^6$	0.32	0.065	$1.40\times {10}^4$	10.4	281
Dwarf 4	${10}^8$	${10}^{10}$	4.0	$2.4\times {10}^7$	0.7	0.14	$2.22\times {10}^5$	164	$4.44\times {10}^3$

Equations

$$T_{\mathrm{SF}}(M_{\mathrm{b}})=T_{\mathrm{SF}}(M_{\mathrm{MW}}){\left(\frac{M_{\mathrm{b}}}{M_{\mathrm{MW}}}\right)}^{\alpha }$$

$$f_{\mathrm{SFR}}(t,M_{\mathrm{b}})=\{\begin{array}{ll}A{M}_{\mathrm{b}}\exp \left(-\frac{t}{\tau (M_{\mathrm{b}})}\right)& \mathrm{if}t<T_{\mathrm{SF}}\\ 0& \mathrm{otherwise}.\end{array}$$

$$M_{\ast }(t,M_{\mathrm{b}})=\{\begin{array}{ll}A{M}_{\mathrm{b}}{T}_{\mathrm{SF}}\left(1-e^{-2t/T_{\mathrm{SF}}}\right)/2& \mathrm{if}t<T_{\mathrm{SF}},\\ A{M}_{\mathrm{b}}{T}_{\mathrm{SF}}\left(1-e^{-2}\right)/2& \mathrm{otherwise}.\end{array}$$

$${\rho }_d(R,z)=\frac{M_{\mathrm{b}}}{4\pi R_d^2{z}_d}\exp \left(-\frac{R}{R_d}\right){\mathrm{sech}}^2\left(\frac{z}{z_d}\right)$$

$$R_d$$

$$z_d$$

$${\rho }_{\mathrm{DM}}(r)=\frac{{\rho }_0}{\frac{r}{r_h}{\left(1+\frac{r}{r_h}\right)}^2},$$

$${\rho }_0=\frac{M_{DM}}{4\pi r_h^3\left(\log (1+C)-\frac{C}{1+C}\right)},$$

$$R_h=32.9{\left(\frac{M_b}{10^8{M}_{\odot }}\right)}^{1/3}kpc.$$

$${\Phi }_{\mathrm{DM}}(r)=-\frac{4\pi {\rho }_{0}G{r}_h^3}{r}\ln \left(1+\frac{r}{r_h}\right),$$

$${\mathbf{a}}_{\mathrm{DM}}(\mathbf{r})=-\nabla {\Phi }_{\mathrm{DM}}(r).$$

$$N_{\mathrm{NSB}}=N_{\mathrm{CCSN}}/x$$

$$p_l(V_{\mathrm{kick},2})=\frac{1}{\sqrt{2\pi }{V}_{\mathrm{kick},2}{\sigma }_{\ln V,l}}\exp \left(\frac{\ln (V_{\mathrm{kick},2}/{\bar{V}}_l{)}^2}{2{\sigma }_{\ln V,l}^2}\right),$$

$$p_{l,h}(\Delta M)\frac{1}{\sqrt{2\pi }\Delta M{\sigma }_{\ln \Delta M,l,h}}\exp \left(\frac{\ln (\Delta M/{\overline{\Delta M}}_{l,h}{)}^2}{2{\sigma }_{\ln \Delta M,l,h}^2}\right),$$

$${\mathbf{V}}_{\mathrm{CM}}$$

$$\frac{a_0}{a}$$

$$1-e^2$$

$$T_{\mathrm{gw}}(a,e,m_1,m_2)=\frac{12c_0^4}{19\beta }\mathcal{G}(e)$$

$$\beta$$

$$c_0$$

$$\mathcal{G}(e)$$

$${\mathbf{V}}_{\mathrm{CM},\mathrm{i}}={\mathbf{V}}_{\mathrm{CM}}+{\mathbf{V}}_c,$$

$${\rho }_{\mathrm{med}}=\max ({\rho }_d({\mathbf{x}}_f),n_0{m}_p{\mu }_a),$$

$$\{\begin{array}{l}{R}_{\mathrm{rem}}(t)=v_{\mathrm{ej}}t0\le t\le t_1,\\ v_{\mathrm{rem}}(t)=v_{\mathrm{ej}}0\le t\le t_1.\end{array}$$

$$\{\begin{array}{l}{R}_{\mathrm{rem}}(t)=R_{\mathrm{rem}}(t_1)+\xi {\left(\frac{E_{\mathrm{rem}}}{{\rho }_0}{\left(t-t_1\right)}^2\right)}^{1/5}{t}_1<t\le t_2,\\ v_{\mathrm{rem}}(t)=\frac{2}{5}\xi {\left(\frac{E_{\mathrm{SN}}}{{\rho }_0}\right)}^{1/5}{\left(t-t_1\right)}^{-3/5}{t}_1<t\le t_2,\end{array}$$

$$t_3={\left(\frac{2}{7v_{\lim }}\frac{R_{\mathrm{rem}}(t_2)}{t_2^{2/7}}\right)}^{7/5},$$

$$\{\begin{array}{l}{R}_{\mathrm{rem}}(t)=R_{\mathrm{rem}}(t_2){\left(\frac{t}{t_2}\right)}^{2/7}{t}_2<t\le t_3,\\ v_{\mathrm{rem}}(t)=\frac{2}{7}\frac{R_{\mathrm{rem}}(t_2)}{t_2}{\left(\frac{t}{t_2}\right)}^{-5/7}{t}_2<t\le t_3.\end{array}$$

$$m_{\mathrm{ret}}=\text{max}\left({\int }_{\Sigma \cap \mathcal{V}}\left(\frac{\mathrm{d}{m}_{\mathrm{ej}}}{\mathrm{d}\Sigma }\right)\mathrm{d}\Sigma ,0.5m_{\mathrm{ej},\mathrm{tot}}\right),$$

$$\sqrt{x_f^2+y_f^2}-R_{\mathrm{rem}}(t_3)$$

$$|z_f-R_{\mathrm{rem}}(t_3)|$$

$$m_{\mathrm{ret}}={\int }_{{\Omega }_f}\left(\frac{\mathrm{d}{m}_{\mathrm{ej}}}{\mathrm{d}\Omega }\right)\mathrm{d}\Omega ,$$

$$X(E_i)=\frac{m_{\mathrm{ej},\mathrm{ret}}(E_i)}{m_{\mathrm{ej},\mathrm{ret}}}=\frac{{\left(Y(E_i)\right)}_{\odot }⟨A_i⟩}{{\sum }_k{\left(Y(E_k)\right)}_{\odot }⟨A_k⟩},$$

$$[{\mathrm{E}}_{\mathrm{i}}/{\mathrm{E}}_{\mathrm{j}}]={\log }_{10}\left(\frac{Y(E_i)}{Y(E_j)}\right)-{\log }_{10}{\left(\frac{Y(E_i)}{Y(E_j)}\right)}_{\odot },$$

Full text

††pubyear: 2019††pagerange: Neutron star binary orbits in their host potential: effect on early r-process enrichment–Neutron star binary orbits in their host potential: effect on early r-process enrichment

1 Introduction

The merger of compact binaries comprising at least one neutron star (NS) has long been thought to be the site for the production of a significant fraction of the heavy elements above the iron group via the so called $r$-process nucleosynthesis (e.g. Lattimer & Schramm, 1974; Eichler et al., 1989; Freiburghaus et al., 1999). The recent detection of a kilonova transient (AT2017gfo) associated with the gravitational wave (GW) signal produced by two NSs in the late phases of their inspiral (GW170817, e.g. Abbott et al., 2017b, a; Drout et al., 2017; Tanaka et al., 2017; Pian et al., 2017; Tanvir et al., 2017; Kasen et al., 2017; Nicholl et al., 2017; Chornock et al., 2017) has finally provided strong observational support to these ideas and confirmed that NSB mergers are one of the major (if not the main) production site for $r$-process nucleosynthesis elements (see e.g. Thielemann et al., 2017; Rosswog et al., 2018; Hotokezaka et al., 2018; Côté et al., 2018).

The unprecedented quality of the kilonova detection has provided a glimpse of the potential variety associated with this new class of transients. The optical and infrared electromagnetic data are well explained by the radioactive decay of $\sim 0.05~{}\rm M_{\large\odot}$ of material (e.g. Cowperthwaite et al., 2017; Rosswog et al., 2018; Tanaka et al., 2017; Tanvir et al., 2017; Smartt et al., 2017; Kasen et al., 2017; Perego et al., 2017). The presence of different peaks in the light curves and the spectral evolution can be modelled by different components in the outflows with different compositions and, possibly, different physical origins. Modelling of the matter outflow and of its properties is presently accomplished by numerical simulations of the merger and of its aftermath. Matter expelled within the first milliseconds after the NS collision is dubbed dynamical ejecta (e.g. Korobkin et al., 2012; Bauswein et al., 2013; Hotokezaka et al., 2013; Wanajo et al., 2014; Sekiguchi et al., 2015; Radice et al., 2016; Bovard et al., 2017; Radice et al., 2018b, c). On longer time scales (a few hundreds milliseconds) a significant fraction of the ejecta can be expelled as baryonic winds that originate from the merger remnant (e.g. Fernández & Metzger, 2013; Metzger & Fernández, 2014; Perego et al., 2014; Siegel et al., 2014; Just et al., 2015; Martin et al., 2015; Siegel & Metzger, 2018; Radice et al., 2018a). This ejecta can come both from the disc that forms around the central remnant (either a massive neutron star or a black hole) and from the central remnant itself, as long as it does not collapse to a black hole.

Observations of $r$-process elements in the atmosphere of metal-poor stars in our galaxy and in nearby classical dwarf galaxies hint at the occurrence of $r$-process nucleosynthesis also in metal-poor environment, corresponding to the very early stages of the galaxy evolution (McWilliam, 1998; Sneden et al., 2003; Shetrone et al., 2003; Honda et al., 2006; François et al., 2007; Sneden et al., 2008; Roederer et al., 2014; Ural et al., 2015; Jablonka et al., 2015; Hill et al., 2018). Moreover, the abundance of europium (an element synthesised mainly by $r$-process nucleosynthesis) presents a large scatter at very low metallicities (${\rm[Fe/H]}<-3$), suggesting that $r$-process elements must be synthesized in rare and isolated events that inject a significant amount of heavy elements into a relatively small amount of gas (e.g. Sneden et al., 2008, and references therein). This rare-event/high-yield scenario is also corroborated by the comparison of iron and plutonium abundances in deep-sea sediments (Hotokezaka et al., 2015). This hypothesis was also tested in stochastic and inhomogeneous chemical evolution models (Cescutti et al., 2015a; Wehmeyer et al., 2015), which succeed to explain the chemical spread of r-process abundances in halo stars, but only assuming very short delay for the NSB merger.

Among the smallest dwarfs, called ultra-faint dwarf galaxies (UFD), Reticulum II (Ji et al., 2016a; Roederer et al., 2016; Ji et al., 2016b) and very likely also Tucana III (Drlica-Wagner et al., 2015; Simon et al., 2017; Hansen et al., 2017; Marshall et al., 2018) show a large excess of $r$-process elements, while for all the other UFDs robust upper limits on their $r$-process element abundances have been set (Frebel et al., 2010b; Simon et al., 2010; Koch et al., 2013; François et al., 2016). The presence of $r$-process elements in UFDs poses serious challenges to NSB mergers as origin of the $r$-process nucleosynthesis elements. First, these galaxies have low escape velocities (Walker et al., 2015). The kick imparted to any newly born NS by its CCSN explosion could potentially eject the NS from the galaxy, even when the stellar binary system survives the second supernova explosion. Second, their old stellar population is thought to be the results of a fast star formation episode, that ends within the first Gyr of the galaxy evolution after the first CCSN explosions expel a significant fraction of baryons (Brown et al., 2014b; Weisz et al., 2015).

A firm understanding of the $r$-process enrichment in dwarf galaxies is also essential to understand metal poor stars in the Milky Way halo. Indeed both observations (Frebel et al., 2010a; Ivezić et al., 2012) and theoretical models (Helmi, 2008; Griffen et al., 2016) point to the fact that the dwarf satellite galaxies that we observe nowadays around our Galaxy are the remnants of a large population of dwarfs that long ago merged with it to form the galactic halo stellar population. Moreover, recent results coming from GAIA (Gaia Collaboration et al., 2016) point to the fact that a large fraction of the Milky Way halo was actually formed through a merger with a single and relative massive satellite (Haywood et al., 2018; Helmi et al., 2018).

Assuming that fast mergers of compact binary systems require a high natal kick, Bramante & Linden (2016) disfavoured NSB mergers as the source of the $r$-process material observed in Reticulum II. In contrast, further studies concerning the observed distributions of the orbital parameters of double NS systems (Beniamini & Piran, 2016) and the formation channels of such systems (Tauris et al., 2017) suggest that a large fraction of NSB systems could have received a rather small kick and have ejected a small amount of mass as a consequence of the second CCSN explosion. These conclusions imply that a large fraction (up to 60%) of double NS systems could be retained even by UFD galaxies and to merge within the first Gyr of galaxy evolution (Beniamini et al., 2016a). Also Safarzadeh et al. (2019) reached opposite conclusions with respect to Bramante & Linden (2016) by considering the possibility that high-kick NSBs are either on highly eccentric orbits or form with very short separations due to an additional mass-transfer between the first-born neutron star and a naked helium star, progenitor of the second neutron star.

Due to the low stellar content of UFD galaxies and the subsequent small number of NSB systems expected in these galaxies, objects like Reticulum II and Tucana III should be a minority. Moreover, in this scenario these UDFs should have hosted a single NSB merger that was able to significantly enrich them in $r$-process material. These first analyses was later refined by Beniamini et al. (2018), who considered that some of the $r$-process material synthesised during the mergers could still escape from the galaxy, thanks to the large kinetic energy that characterizes NSB merger ejecta ($\sim 10^{50}$ erg). By performing a more detailed analysis of the enrichment in iron from CCSNe and in $r$-process elements from NSB mergers in UFD galaxies, they confirmed the compatibility between the abundances observed and an $r$-process enrichment due to the ejection of material from rare events taking place inside dwarf galaxies with low escape velocities.

In many of the above mentioned studies, it was assumed that a NSB bound to its host coalesces always well inside the galaxy. The potential relevance of the merger location, and thus of the imparted kick, in explaining the abundances observed in UFDs was first underlined by Safarzadeh & Scannapieco (2017). In particular, if the merger time is not extremely short (Bonetti et al., 2018; Safarzadeh et al., 2019), the binary will start orbiting the galaxy and there is a high chance that the merger will still happen far from the regions of the galaxy where the next generation of stars will form. In this work, we systematically explore the potential effect of the dilution of the ejecta on enrichment of the $r$-process material due to the location of NSB mergers relative to the host galaxy (see also Safarzadeh & Côté, 2017, for a similar analysis focused on the Milky-Way). Under a common set of minimal assumptions about the properties of the host galaxy, the star formation rate inside it, the NSB birth and coalescence as well as the ejecta properties, we investigate a wide sample of galaxy masses and we compute the fraction of the ejecta material retained by the galaxy.

The paper is structured as follows: in sections 2-5 we present the model we have adopted for our calculations and we detail all its components. The results we have obtained are presented in section 6, while in section 7 we discuss our results and compare them with observations. We finally conclude in section 8.

2 Overview of the chemical enrichment model

Modelling the evolution of galaxies and of their chemical enrichment is an extremely complex task, since it requires to follow many different processes that span a huge range of scales (both in space and in time). In this section, we first present a summary of the (simplified) model adopted in this work to study the enrichment in $r$-process material due to NSB mergers.

We consider disc galaxies with a different content of baryonic matter $M_{\rm b}$, ranging from $10^{5}M_{\odot}$ up to $10^{8}M_{\odot}$, as well as a model with $5\times 10^{10}M_{\odot}$. The former interval is expected to correspond to the initial gas content of ultra-faint and classical dwarf galaxies while the latter value is the one of a MW-like galaxy. During the cosmic history, gas is converted into stars with a certain star formation rate ($f_{\rm SFR}$) starting from the galaxy formation ($t=0$) up to $t\approx T_{\rm SF}$. For a MW-like galaxy, $T_{\rm SF}\sim 14~{}{\rm Gyr}$ and the final stellar content is thought to be comparable to the initial gas mass. For smaller galaxies SN feedback and/or environmental processes (e.g. tidal perturbations and ram pressure stripping) are expected to quench the conversion of gas into stars on shorter timescales and to remove a fraction of the gas from the galaxy (Revaz & Jablonka, 2018). In the case of UFD, we assume $T_{\rm SF}\lesssim 1~{}{\rm Gyr}$, while for classical dwarfs $T_{\rm SF}\sim 3-5~{}{\rm Gyr}$. The SFR has a non-trivial dependence on the cosmic time and on the individual history of each galaxy. Since we are not interested in the detailed time evolution of the metal content nor in the reproduction of a specific galaxy model, we adopt an exponential dependence on time, possibly dependent on the initial gas mass. The stellar content at any time $t$, $M_{*}(t)$ is then computed as the integral of $f_{\rm SFR}$ over time.

For each galaxy, we generate a pool of $N$ initial conditions for NSBs forming from stellar binaries as a consequence of a double CCSN explosion. For each NSB, the relevant initial conditions after the second SN are the initial position in the galaxy, the NS masses ($m_{1},m_{2}$), the semi-major axis and eccentricity ($a,e$), and the center of mass (CoM) velocity of the NSB with respect to galaxy frame ($\mathbf{V}_{\rm CM}$). For all the dwarf galaxies in the our sample, the expected number of merger is $\lesssim 10^{4}$, therefore we choose $N\gg N_{\rm merg}$ to properly span the parameter space. For the MW-like galaxy, the number of mergers could potentially exceed a few millions. We expect such a high number of configurations to be large enough to significantly cover the whole NSB parameter space. Thus we consider a pool of $N\sim 2\times 10^{6}$ NSB configurations. We evolve each initial condition by integrating the trajectory of the CoM of the NSB for the GW driven coalescence time $T_{\rm gw}(a,e)$, unequivocally determined by the initial binary parameters. Thus, we can associate to each NSB in our pools the time at which the merger happens and the corresponding location.

To compute the amount of $r$-process material produced by NSB mergers that has enriched a specific galaxy at a time $t$ we proceed as follows:

• •

starting from the galaxy stellar mass at $t$, $M_{*}(t)$, we estimate the amount of CCSN explosions, $N_{\rm CCSN}$, and from that the expected number of stellar binaries surviving a double CCSN explosion and producing a NSB system within $t$, $N_{\rm merg}$;

• •

we sample the actual number of NSB mergers $N_{r}$ from a Poisson distribution with average equal to $N_{\rm merg}$, and we randomly choose $N_{r}$ cases from the galaxy NSB pool. We stress here that a (possibly large) fraction of these binaries could merger on a time $T_{\rm gw}\lesssim t$. In the case of the MW-like galaxy, if $N_{r}>N_{r,{\rm max}}=10^{5}$ we select $N_{r,{\rm max}}$ NSB from the pool of $N\sim 2\times 10^{6}$ elements. In this case, at the end of the analysis, the true values of the total and retained masses are obtained by rescaling the computed quantities by a factor $N_{r}/N_{r,{\rm max}}$. This approach is motivated by the need of reducing the necessary computations and the amount of data. However, it is justified by the large number of expected mergers with respect to their intrinsic variability;

• •

for each sampled binary merging within $t$, we model the properties of the material ejected during the coalescence, and its expansion due to the interaction with the diffuse gaseous halo of the host in order to determine the fraction and composition of the ejecta that gets injected in the galaxy disc. We consider that the ejected mass can be in form of wind ejecta (isotropically distributed) as well as dynamical ejecta, characterized by both a polar and an equatorial component.

Finally, in order to estimate both the average and the scatter in the distribution of the total ejected and retained material, for each galaxy model we perform 200 different realizations of the expected NSB populations. For each realization we assume full mixing of the ejecta. Despite being a simplification, this hypothesis is not a limitation because we will compare with estimates of the total galactic abundances. More detailed studies assuming inhomogeneous mixing and stochastic stellar enrichment are planned for the future.

In the following we fully detail the procedure adopted in our model. In section 3, we start with the modeling of the host galaxy and with the determination of the number of stars, SNae and merging NSBs for each galaxy model. We also discuss the modelling of the host galaxy potential, which determines the velocity of the NSB progenitor binary system, its CoM orbital evolution as well as the escape velocity. We then discuss the initialization of the intrinsic parameters of each NSB (section 4.1): the masses of the two NSs ($m_{1}$ and $m_{2}$), the semi-major axis and eccentricity of the the binary ($a$ and $e$), and the kick that the NSB gets due to the two SNae explosions ($\mathbf{V}_{\rm CM}$). The initialization of the initial position and velocity of the NSB CoM is described in section 4.2. The ejecta properties immediately after the NSB coalescence are described in section 4.3, while their evolution within the host halo is described in section 4.4.

3 Galaxy models

We consider five different disc galaxies spanning a wide range of possible masses, $M_{\rm b}$: four dwarf galaxy models with total baryonic mass of $10^{5}~{}\rm M_{\large\odot}$, $10^{6}~{}\rm M_{\large\odot}$, $10^{7}~{}\rm M_{\large\odot}$, and $10^{8}~{}\rm M_{\large\odot}$ respectively, and a MW-like host modelled with five dynamical baryonic components (galactic bulge + thin and thick stellar disks + HI and H2 disks, see next and Barros et al., 2016, for full details) surrounded by a dark halo.

The stellar population of UFD galaxies seem to be dominated by very old ($\sim 12~{}{\rm Gyr}$) stars (Brown et al., 2014a, and references therein), pointing to a rapid SF, $\lesssim 1~{}{\rm Gyr}$. Old stars ($\lesssim 10~{}{\rm Gyr}$) dominate also dwarf galaxies (Grebel, 1997), indicating that the bulk of the SF in them happens within 2-3 Gyr. To approximately catch the dependence of the duration of the star formation on the initial baryonic mass $M_{\rm b}$, we assume a power-law dependence normalized to the MW-like case:

[TABLE]

with $\alpha=0.2$ and $T_{\rm SF}(M_{\rm MW})=14~{}{\rm Gyr}$. For the star formation rate, we adopt the following exponential dependence:

[TABLE]

where $A$ is a constant fixed by the requirement that the MW-like model reproduces the presently observed star formation rate in the MW, i.e. $f_{\rm SFR}(10~{}{\rm Gyr},M_{\rm MW})=1.65~{}~{}\rm M_{\large\odot}~{}{\rm yr^{-1}}$ (Licquia & Newman, 2015). For the timescale appearing inside the exponential factor, we choose $\tau=T_{\rm SF}/2$. This choice is broadly compatible with (simple) models of the MW (e.g. Snaith et al., 2014) and of classical dwarf galaxies (North et al., 2012). Finally, the gas mass converted into stars as a function of time and of initial baryonic mass can be easily computed as

[TABLE]

3.1 Dwarf galaxies

We assume that the potential of the host is well described by only two components: a baryonic disc (where we assume that the stars and gas follow the same profile) and a dark matter spherical halo.

The baryonic density profile of the galaxy is modelled as an exponential disk:

[TABLE]

where $R$ and $z$ are the cylindrical radial and vertical coordinates, while the length scales $R_{d}$ and $z_{d}$ are (Mo et al., 1998)

[TABLE]

An analytic form for the potential (and consequently for the acceleration) cannot be obtained, we therefore employ a numerical sampling of the disk density profile with $\approx 8\times 10^{5}$ tracers, and splitting every tracer in 8 sub-tracers by changing the sign of 1, 2 or all the 3 coordinates, in order to preserve the symmetry of the potential. The disc acceleration $\mathbf{a}_{\rm disc}(\mathbf{r})$ is then evaluated through direct summation over all the sampled particles, where we use a gravitational softening of $\epsilon_{\rm soft}=0.01\left(M_{b}/10^{8}~{}\rm M_{\large\odot}\right)^{1/3}$ kpc to avoid spurious strong scattering due to the finite number of tracers used (see e.g. Monaghan & Lattanzio, 1985). Such a kind of exponential disk profiles111Or equivalently Sérsic profiles with a very small Sérsic index fits well the stellar brightness profile of many observed dwarf galaxies, including dwarf ellipticals and dwarf spheroidals, (see e.g. Faber & Lin, 1983; Binggeli et al., 1984; Kormendy, 1985; Graham, 2002; Graham & Guzmán, 2003). Furthermore, it implies that the stellar dynamics is dominated by rotation. While this seems to be the case for the vast majority of dwarf galaxies (Kerr et al., 1954; Kerr & de Vaucouleurs, 1955; Swaters et al., 2009), for dwarf spheroidal galaxies with less rotational support and embedded in virialized galaxy clusters it is reasonable to assume that they were more rotationally supported in the past (during their star formation epoch) and lost their coherence due to galaxy harassment (Moore et al., 1996).

The dark matter density is assumed to follow a NFW profile:

[TABLE]

where the normalization $\rho_{0}$ is given by

[TABLE]

with $M_{\rm DM}=100M_{b}$,222Dwarf galaxies are among the structures with the smallest baryonic to dark matter ratios. Detailed studies (e.g. McGaugh et al., 2010; Chan, 2019) suggest a weak dependence on the baryonic mass, $\sim M_{b}^{0.23}$. For simplicity in our study we assume one single value broadly compatible with the expected ratio in dwarf galaxies. $C=9.4$ and the scale radius is given by:

[TABLE]

The potential generated by such distribution is analytic

[TABLE]

which allow us to compute the acceleration directly from

[TABLE]

The total acceleration is then computed as $\mathbf{a}=\mathbf{a}_{\rm DM}+\mathbf{a}_{\rm disc}$.

3.2 MW-like galaxies

Following Barros et al. (2016), we model the MW by considering a dark halo, a galactic bulge plus two stellar disks (thin and thick) and two gaseous disks (HI and H2). The dark halo is modelled as a logarithmic potential, while the bulge is assumed spherically symmetric and described by an Hernquist profile (Hernquist, 1990). The description of the baryonic disks requires instead additional modelling. In fact, despite from an observational point of view these disks can be accurately fitted assuming exponential disks, as pointed out in section 3.1, they do not admit any analytical form for the gravitational potential nor for the acceleration. This does not represent a real issue for the dwarf galaxy case because the global number of NSB that we have to simulate is rather small. On the contrary, for a MW-like galaxy this number can exceed several millions, with a severe impact on the performance of our calculation. This motivated our choice to describe the baryonic disks with analytical potential-density pairs. In particular, in Barros et al. (2016), each of the four disk is modelled with a superposition of three Miyamoto-Nagai (MN) disks (Miyamoto & Nagai, 1975), where one of the three mass parameter is usually negative in order to mimic the sharp decrease that characterise the density profile of exponential disks. We address the interested reader to section 2 of Barros et al. (2016) for the detailed description of the model, and to table 3-4 in the same paper for the values of the parameter that best describe the MW and that we implemented in our model.

3.3 Average number of NSB mergers for a given galaxy

Given the large uncertainties that affect the determination of this value, we adopt a fairly simple but physically motivated approach. We assume a standard stellar initial mass function (IMF), identical for all models (Kroupa, 2001). For a given galaxy model and a given time, we compute the number of stars that have exploded as CCSNae ($N_{\rm CCSN}$) as the number of stars with mass greater than $8~{}\rm M_{\large\odot}$. In doing that, we are implicitly assuming that the evolution timescale of massive stars is much shorter than any time we are going to explore. We then assume that the number of NSB systems that form represents a fraction of the total number of SNe that have exploded. In particular, we parametrise $N_{\rm NSB}$ as

[TABLE]

where for $x$ we choose four values logarithmically distributed between an optimistic and a very pessimistic estimate, i.e. $x=[50,150,450,1350]$. This broad interval covers present uncertainties in the determination of the ratio between the number of exploding CCSNe and the number of forming NSB systems, as obtained in detailed population synthesis models (Giacobbo & Mapelli, 2018).

In Table 1 we summarize the properties of the different dwarf galaxy models, alongside their names used in the following.

4 Binary NS and their ejecta

4.1 NSB internal parameters and kick velocity

Both the masses and the semi-major axes of each binary are sampled from observationally constrained distributions (e.g. Tauris et al., 2017, and references therein). The masses of the two NSs $m_{1}$ and $m_{2}$ are randomly sampled from normal distribution with $\mu_{1}=1.4~{}\rm M_{\large\odot}$, $\mu_{2}=1.34~{}\rm M_{\large\odot}$ and $\sigma=0.14~{}\rm M_{\large\odot}$. The semi-major axis $a$ is evaluated from a log-uniform distribution where upper and lower limit are set selecting twice the maximum and half the minimum of “observed” NS binary orbital periods, i.e. $P\in[0.05,100]$ days.333Assuming a total mass of $\sim 2.8~{}\rm M_{\large\odot}$ these limits correspond to a minimum and maximum semi-major axis of about $\sim 0.8$ and $\sim 130$ solar radii.

Tidal circularization acting after the first SN explosion does not allow to constrain the kick velocity experienced by the first-born NS. As a consequence, we sample the first kick experienced by the binary CoM ($\mathbf{V}_{\rm CM,1}$) assuming isotropy and a uniform magnitude distribution between 10 and 20 km/s, as constrained by the observed velocities of High-Mass X-Ray Binaries (see e.g. Coleiro & Chaty, 2013).444For this reason, we do not need to consider the effect of the mass lost during the first SN.

The kick velocity associated to the second SN ($\mathbf{V}_{\rm kick,2}$) is isotropically generated in the $m_{2}$ rest frame. The kick magnitude ($V_{\rm kick,2}$) is distributed as follows: following Beniamini & Piran (2016), we assume that the distribution is bimodal, with 60-70$\%$ of the cases following a log-normal "low-kick distribution"

[TABLE]

with $\bar{V}_{l}=9$ km/s and $\sigma_{\ln{V},l}=0.8$, and the remaining 40-30$\%$ of the cases are sampled from a second log-normal distribution $p_{h}(V_{\rm kick,2})$ with the corresponding parameters $\bar{V}_{h}=158$ km/s and $\sigma_{\ln{V},h}=0.5$.555The value of the chosen parameters are selected in order to broadly reproduce the merger time and the eccentricity distribution of observed NSB.

The mass loss associated to the second SN is also sampled from a log-normal distribution:

[TABLE]

with $\bar{\Delta M}_{l}=0.3~{}\rm M_{\large\odot}$, $\bar{\Delta M}_{h}=1~{}\rm M_{\large\odot}$, and $\sigma_{\ln{\Delta M},l}=\sigma_{\ln{\Delta M},h}=0.5$. Conservation of linear momentum implies that the velocity acquired by the CoM is:

[TABLE]

where the second and third terms on the right hand side of the equation correspond to the contributions to the NSB CoM velocity from the kick experienced by the second remnant and from the (assumed instantaneous) mass change of the binary system (Postnov & Yungelson, 2014).

The last parameter we estimate is the eccentricity $e$ of the NSB immediately after the second SN. $e$ is obtained assuming the conservation of energy and angular momentum after the second SN:

[TABLE]

where $\chi=((m_{2}+\Delta M)+m_{1})/(m_{2}+m_{1})\geq 1$ is the fractional change in mass, and $V_{\rm kep}=\sqrt{G((m_{2}+\Delta M)+m_{1})/a_{0}}$.666Note that the assumption here on the direction of the secondary velocity immediately before the second SN is only made for clarity reasons. In the actual sampling $\mathbf{V}_{\rm kep}$ is assumed isotropic..

The time to coalescence due to GW emission $T_{\rm gw}$ is then evaluated from $(a,e,m_{1},m_{2})$ by performing the integration (Peters, 1964)

[TABLE]

where

[TABLE]

In Fig. 1 we present the probability distribution of the merger times, $t_{\rm m}$ computed as equation 18, for a large set ($10^{4}$) of NSBs whose properties have been generated using the above mentioned distributions for the internal parameters and kick velocities. In particular, we consider the case $p_{\rm l}=0.7$. In 42% of the cases the NSB merges within the Hubble time. This percentage decreases to 27% and 14% if we limit $t_{\rm m}$ to 1 Gyr and 0.1 Gyr, respectively. We thus notice that our choice of a log-uniform distribution in the binary semi-major axis down to $a\lesssim R_{\odot}$ and the inclusion of highly eccentric cases (mainly for the high kick population) introduce a significant population of fast merging binaries (see e.g. Belczynski et al., 2002; O’Shaughnessy et al., 2008; Beniamini & Piran, 2016). However, this percentage is possibly smaller (by a factor $\lesssim 2$) than the one implied by the observed distribution of Galactic NSBs (Beniamini & Piran, 2019).

4.2 Binary NS external parameters and dynamical evolution

We sample the position of the CoM of the NSB at the time of the second SN explosion from the stellar disk profile (i.e. the exponential disk for dwarf galaxies or the superposition of MN for MW case), and we assume an initial CoM velocity

[TABLE]

where $\mathbf{V}_{\rm CM}$ is the velocity acquired due to the SN kicks (see equation 15), while $V_{c}=\sqrt{R\,{\rm d}\Phi/{\rm d}r}$ is the circular velocity in the galactic plane (i.e. at $z=0$) at the radius at which the NSB is initialized. We then integrate the motion of the CoM of the NSB in the chosen galactic system by numerically solving the equations of motion with an 8th order adaptive stepsize Dormand-Prince-Runge-Kutta algorithm (DOPRI853). We stop the integration at $T_{\rm gw}$ and we denote the NSB position at that time, $\mathbf{x}_{f}$, as the merger location. As a final step, we sample the direction of the NSB angular momentum from an isotropic distribution in order to infer the orientation of the binary orbital plane with respect to the galactic reference frame.

4.3 Binary NS ejecta: initial conditions and nucleosynthesis yields

The NSB is assumed to eject a certain amount of mass during and immediately after the merger. We consider two kinds of ejecta: dynamical and disc wind ejecta.

4.3.1 Dynamical ejecta

Dynamical ejecta is expelled on a timescale of a few milliseconds due to tidal torques and hydrodynamics shocks. To model the properties of the dynamical ejecta, we use the parametric fit reported in Radice et al. (2018b) for the total ejecta mass and average speed. To compute the NS compactness parameters associated with the two NS masses $m_{1,2}$ and required by these fits, we choose a nuclear equation of state compatible with all present nuclear and astrophysical constraints (SLy, Haensel & Potekhin, 2004). Numerical simulations of NSB mergers show that this mass ejection happens preferentially along the equatorial plane. We thus adopt a $\sin^{2}{\theta}$ dependence of the mass spatial distribution on the polar angle measured with respect to the NSB rotational axis (Perego et al., 2017; Radice et al., 2018b). Moreover, neutrino irradiation, more intense inside the polar funnels, increases the electron fraction above 0.25 for $\theta\lesssim\pi/4$ and $\theta\gtrsim 3\pi/4$. As a consequence, we assume that polar ejecta produces $r$-process nucleosynthesis yields between the first and second $r$-process peaks, for which we consider an atomic mass number interval $70\leq A\leq 120$. On the other hand, the more neutron rich equatorial ejecta produces elements between the second and third $r$-process peaks for which we assume an $r$-process nucleosynthesis interval $120\leq A\leq 232$.

4.3.2 Disc wind ejecta

The ejection of matter in the form of disc winds happens on timescales longer than the ones of the dynamical ejecta (up to a few hundreds milliseconds). It is due to neutrino absorption, magnetic and viscous processes inside the remnant. To estimate the total mass contained inside this ejecta, we assume that a fixed fraction of the disc $\xi_{\rm wind}$ becomes unbound and is launched with an average velocity $v_{\rm disk}$. In this work we assume $\xi_{\rm wind}=0.2$ and $v_{\rm disk}=0.08{\rm c}$ (see e.g. Just et al., 2015; Fahlman & Fernández, 2018). The disc mass is determined through the fitting formula reported in Radice et al. (2018b). As in the case of the dynamical ejecta, the calculation of the dimensionless tidal coefficient of the merging NSB required by this fitting formula is computing assuming the same SLy nuclear equation of state used for the dynamical ejecta. This ejecta is expected to be more isotropic, both in terms of mass and electron fraction distribution. However, neutrinos could also affect its nucleosynthesis (Perego et al., 2014; Martin et al., 2015; Lippuner et al., 2017). In particular, if the merger results in a long-lived massive NS, the ejecta electron fraction could be systematically shifted above 0.25 such that the production of heavy $r$-process elements is prevented and the nucleosynthesis produces only nuclei between the first and second $r$-process peaks, i.e. for $70\leq A\leq 120$. Otherwise, if a black hole forms promptly or on a timescale smaller than the disc viscous timescale, the production of all $r$-process elements is foreseen. In this case we assume a mass number interval $80\leq A\leq 232$ for the $r$-process nucleosynthesis. To distinguish between the first and the second case we use an empirical threshold value suggested by the CoRe NSB merger database (Dietrich et al., 2018): if $(m_{1}+m_{2})<1.3M_{\rm NS,max}$ (where $M_{\rm NS,max}\approx 2.05~{}\rm M_{\large\odot}$ is the maximum cold NS mass predicted by the SLy nuclear equation of state) then the central remnant does not collapses to a BH before the wind ejecta is expelled and the production of all $r$-process elements in the disk ejecta is prevented.

4.4 Binary NS ejecta: evolution

The ejecta expanding in the intergalactic (IGM) or interstellar (ISM) medium will slow down, forming a NSB merger remnant similar to the remnant produced by SN explosions (see e.g. Montes et al., 2016). This will eventually mix with the surrounding medium, enriching it with its nucleosynthesis yields. We assume the IGM/ISM to follow the same density profile as the model galaxy down to a floor number density $n_{0}$, and to be formed by hydrogen atoms with a temperature of $T\sim 10^{4}{\rm K}$. The relation between the mass and the number floor densities is given by $\rho_{0}=n_{0}m_{p}\mu_{a}$, where $m_{p}$ is the proton mass and $\mu_{a}$ the mean atomic weight, $\mu_{a}=1.27$. Values of $n_{0}$ around galaxies are largely unknown and strongly dependent on the environment and cosmological epoch. We consider a fiducial value $n_{0}=10^{-4}{\rm cm^{-3}}$, larger but still comparable to the present average density of the Universe ($\sim 2\times 10^{5}{\rm cm^{-3}}$). To explore the possible impact of larger densities, expected for example in the early Universe, we investigate also $n_{0}=10^{-2}{\rm cm^{-3}}$. We notice that the resulting interval $10^{-4}{\rm cm^{-3}}\lesssim n_{0}\lesssim 10^{-2}{\rm cm^{-3}}$ is compatible with the values inferred by the GRB afterglow emission of GW170817 (Margutti et al., 2017; Ghirlanda et al., 2019).

To model the evolution of the remnant, we assume that the two kinds of ejecta (characterized by different initial kinetic energies) will produce two remnants, that we treat independently. Moreover, for simplicity, we consider the remnant expansion to happen inside a uniform medium of density

[TABLE]

where $\mathbf{x}_{f}=(R_{f},z_{f})$ represents the merger location in the galactic reference frame with $\rho_{d}(\mathbf{x}_{f})$ denoting the baryonic density (cf. section 3.1) at such point. If $R_{\rm rem}$ and $v_{\rm rem}$ denote the radius and the speed, respectively, of the NSB remnant forward front, we model their time evolution as in the case of SN remnants and we follow an approach similar to the one described in Haid et al. (2016) and in Beniamini et al. (2018). More specifically, we consider the following expansion phases:

• •

free expansion: during this phase, the ejecta expands with constant velocity $v_{\rm ej}$ equal to either $v_{\rm dyn}$ or $v_{\rm disk}$, depending on the nature of the ejecta considered:

[TABLE]

This phase lasts up to $t_{1}=\frac{1}{v_{\rm ej}}\left(\frac{3m_{\rm ej}}{4\pi\rho_{0}}\right)^{1/3}$, i.e. the point when the ejecta has swept-up a mass comparable to its total mass.

• •

Sedov-Taylor expansion: this phase lasts up to the point where radiative cooling becomes relevant. This transition usually happens at $t\sim t_{\rm TR}=4\times 10^{4}\,{\rm yr}\,n_{0}^{-0.53}$ (Haid et al., 2016). The phase $t_{\rm TR}\leq t<C_{\rm TR}t_{\rm TR}$, with $C_{\rm TR}\approx 1.83$, is an intermediate phase in which the expansion cannot be expressed as a self-similar solution. In this first study we neglect this additional complication and assume that the self-similar Sedov-Taylor expansion continues up to $t_{2}=C_{\rm TR}t_{\rm TR}$

[TABLE]

where $\xi=25/(4\pi)$ and $E_{\rm rem}\approx m_{\rm ej}v_{\rm ej}^{2}/2$ is the initial kinetic energy of the ejecta.

• •

snowplow expansion: in this phase, the ejecta has produced a thin shell, containing most of the ejecta mass, that expands driven by the hot interior (pressure-driven snowplow phase). This phase is characterized again by a self-similar solution, and continues up to the point where the ejecta velocity equals the sound speed of the IGM/ISM, $v_{\rm lim}\approx 10~{}{\rm km/s}$. If we denote this time as,

[TABLE]

then

[TABLE]

The maximum remnant expansion is assumed to occur at $t_{3}$, since after that the ejecta dissolves into the IGM/ISM.

5 Retained ejecta and elemental abundances

5.1 Fraction of retained ejecta

Let us consider a NSB that merges at a certain location $\mathbf{x}_{f}=(x_{f},y_{f},z_{f})=(R_{f},z_{f})$. If the resulting remnant is large enough (or alternatively the merger point is not too far from the host galaxy), a certain fraction of the ejecta can be retained by the galaxy. Moreover, if the NSB inspiral and the remnant evolution are fast enough (i.e. $T_{\rm GW}+t_{3}\lesssim T_{\rm SF}$), the retained mass will pollute the gas that will form a new generation of stars. To compute the fraction of retained ejecta, the first step is to define a suitable contour for the galaxy. We define it as the iso-density contour in the baryonic density profile $\rho_{\rm lim}=\rho_{d}(\xi_{\rm lim}R_{\rm d},0)$ and we consider the surface $\mathcal{S}$ (enclosing the volume $\mathcal{V}$, i.e. the galaxy) defined by all points $(R_{p},z_{p})$ that satisfy the condition $\rho_{d}(R_{p},z_{p})=\rho_{\rm lim}$. We set $\xi_{\rm lim}=3$ and $\xi_{\rm lim}=5$ in order to explore the impact of the chosen threshold on the retained mass.777Here we comment on the dwarf galaxy case only (i.e. $\rho_{d}$ is given by the exponential disk), with the understanding that the procedure followed for the MW is conceptually the same. If $\rho_{\rm med}$ denotes the density inside which the remnant expands (see previous subsection), these are the possible scenarios (see Fig. 2 for a cartoon sketch):

• •

if $\rho_{d}(\mathbf{x}_{f})>\rho_{\rm lim}$, the coalescence happens within the host. When the ejecta reaches the maximum expansion radius, $R_{\rm rem}(t_{3})$, we evaluate the intersection of the corresponding spherical shell $\Sigma$ with the volume $\mathcal{V}$ and we estimate the amount of the retained mass as the fraction of the spherical surface that lies inside the galaxy times the total ejected mass, i.e.

[TABLE]

where ${\rm d}m_{\rm ej}/{\rm d}\Sigma$ is the ejecta angular distribution on the remnant sphere. The maximum with $0.5m_{\rm ej,tot}$ ensures that at least half of the ejected mass remains inside the galaxy. This is due to the fact that the increasing density gradient moving from $\mathbf{x}_{f}$ towards the galactic disc plane will certainly slow down the ejecta moving in that direction.

• •

if $\rho_{d}(\mathbf{x}_{f})<\rho_{\rm lim}$ and $\Sigma$ does not expand across the galaxy volume $\mathcal{V}$, then none of the ejecta can enrich the galaxy. Practically, the situation of non-intersection verifies when one of conditions below is satisfied

[TABLE]

where $\xi_{\rm lim}R_{d}$ is the largest disk contained inside $\mathcal{V}$, while the interval $[-z_{\rho_{\rm lim}},z_{\rho_{\rm lim}}]$ represents the maximum $z$ extent of the galaxy at a radial distance $R_{f}=\sqrt{x_{f}^{2}+y_{f}^{2}}$.

• •

if $\rho_{d}(\mathbf{x}_{f})<\rho_{\rm lim}$ and the remnant sphere extends inside or beyond the galaxy, then a certain amount of the ejected mass can pollute the galaxy. In order to estimate this mass, we identify the fraction of solid angle that intersects the galaxy, $\Omega_{f}$, using a Monte Carlo approach.888To this end, we generate $\sim 4\times 10^{6}$ propagation directions starting from $\mathbf{x}_{f}$ and spanning the whole solid angle. The retained mass is then computed as

[TABLE]

where ${\rm d}m_{\rm ej}/{\rm d}\Omega$ is the ejecta angular distribution.

5.2 r-process nucleosynthesis abundances

As stated in section 4.3, NSB merger ejecta is characterized by a heterogeneous composition of heavy $r$-process elements potentially produced in different components. In order to asses the contribution of each element to the mass retained by the host galaxy, we consider that the relative abundances of heavy elements closely follows the Solar System (SS) abundances obtained from Lodders (2003) and decomposed in its $s-$ and $r-$process contributions following Sneden et al. (2008). If $A_{\rm min}$ and $A_{\rm max}$ are the minimum and maximum mass number produced by the $r$-process nucleosynthesis in a certain ejecta component, the mass fraction of each element Ei (whose isotopes are such that $A_{\rm min}\leq A_{i}\leq A_{\rm max}$) inside the retained mass of each component can be expressed as

[TABLE]

where $Y({E_{i}})$ is the elemental abundance of Ei, $\langle A_{i}\rangle$ its mean mass number, as obtained by the SS abundances, and the sum at denominator ranges over all produced elements. 999Throughout the literature, the SS abundance of an element Ei is usually quoted in terms of the astronomical logarithmic scale $A(\textrm{E}_{i})=\log_{10}(n_{i}/n_{H})+12$, with $n_{H}$ the number density of Hydrogen atoms or according to the cosmochemical scale, which instead normalises the number of Silica atoms to $10^{6}$. The relation between the two scale is given by $A(\textrm{E}_{i})=1.54+\log_{10}N_{\rm cosmo}(E_{i})$. Once $A_{\rm min}$ and $A_{\rm max}$ are given, $X(E_{i})$ is unambiguously determined and allows us to obtain the mass composition of the $r$-process nucleosynthesis retained material in each component.

The total retained mass of each element is then obtained as the sum of the element retained masses in all components. Finally, the abundance ratio of two elements, defined as

[TABLE]

can be equivalently computed as

[TABLE]

where the last term is a constant determined by the SS abundances.

In addition to the retained mass in $r$-process elements, we also compute the iron mass produced by CCSNe. According to observations, 2/3 of CCSNe explode as type II SN and produce, on average, 0.02 $~{}\rm M_{\large\odot}$ of iron; 1/3 are SNIbc and produce, on average, 0.2 $~{}\rm M_{\large\odot}$ of iron (Li et al., 2011; Drout et al., 2011). For each galaxy realization, if $N_{\rm CCSN}$ is the number of CCSNe occurred within a time $t$, we sample $N_{\rm CCSNII}$ type II SNe from a Poisson distribution with average equal to $2N_{\rm CCSN}/3$, and compute $N_{\rm CCSNIbc}=N_{\rm CCSN}-N_{\rm CCSNII}$. Since CCSNe are expected to explode inside the galaxy, we further assume that all the iron is retained inside the galaxy and we estimate its amount from the average quantity produced per event. More detailed studies taking into account the kinetic energy of the ejecta as well as the elemental mixing (Beniamini et al., 2018; Emerick et al., 2018), showed that the amount of retained elements could sensitively depend on the source. Thus, our values provide upper limits to the iron enrichment due to CCSNe. For example, Beniamini et al. (2018) estimated a retain factor of 0.2-0.9 for ejecta products exploded within a dwarf galaxy.

6 Results

In this section we present the result obtained when all the ingredients presented in sections 2-5 are combined. We first focus on a fiducial case characterised by $n_{0}=10^{-4}{\rm cm^{-3}},\xi_{\rm lim}=3,x=150,p_{l}=0.7$ in Sec. 6.1. We then discuss the impact of parameters by varying each of them one at a time in Sec.6.2.

6.1 Model exploration

A first qualitative interpretation of our results is given by Fig. 3, where we report the $x-z$ projection of the NSB merger locations (black dots in the figure) as well as the expansion sphere (circles in the projections) of the ejected $r-$process material. In this plot and for the other that follow, we consider, from left to right, three galaxy models with total baryonic mass of $10^{5},10^{7},\sim 5\times 10^{10}~{}\rm M_{\large\odot}$, indicatively representative of UFDs, classical dwarf and MW-like galaxies. From the figure it is evident that two factors primary determine the level of enrichment of the progenitor galaxy (dark blue contour in the figure): the number of NSBs that form in the galaxy and merge within $T_{\rm SF}$, and the distance from the galaxy at which each NSB travels prior to coalescence. Both factors crucially depend on the baryonic mass of the galaxy. In very low mass galaxies (e.g. left panel of Fig. 3) the enrichment level is quite low since the average number of formed NSB is usually small, around unity or fractions of it. In addition, the potential well provided by the galaxy is rather shallow, therefore the SN kicks can imprint enough velocity such that a NSB system can substantially recede from the parent galaxy. Even if a NSB is gravitationally bound to the galaxy, it often spends most of its orbital time rather far from the galactic disc. On the contrary, for massive galaxies (e.g. MW-like ones, right panel of Fig. 3) the level of enrichment is much higher given the increased number of massive stars that can produce NSB systems, the longer $T_{\rm SF}$, as well as the higher escape velocity that prevent NSB to travel far away from the galaxy (the energetics of the SNae is reasonably assumed to be the same irrespective of the galaxy mass). Between the above situations, galaxies with intermediate masses continuously connect the two extremes, as visible for example in the central panel of Fig. 3.

This can be more quantitatively inferred from Fig. 4, in which we report the mass distribution of the whole amount of produced $r-$process material (dashed magenta lines) compared to the $r-$process mass actually retained by the host galaxy (solid black lines). For each galaxy we also consider only NSBs that coalesce within the time available for star formation ($T_{\rm SF}$), since a galaxy with no more star formation activity cannot form stars enriched with heavy elements.

In the least massive galaxy case ($M_{\rm b}=10^{5}~{}\rm M_{\large\odot}$, left panel), only a small fraction of the realizations ($<10\%$) presents at least one NSB merger happening within $T_{\rm SF}=1~{}{\rm Gyr}$. The small probability that the corresponding ejecta intersects the galactic disc further reduces the probability that a galaxy realization displays a $r$-process material enrichment ($\sim 2\%$). For more massive galaxies the distributions are more peaked with an higher number of realisations providing some enrichment, meaning that in these systems NSB mergers are more likely and generally closer to the galaxy. Moreover, from the figure it is particularly clear that there exists a dilution factor between the $r-$process mass produced by NSB mergers and the actual amount captured by the host galaxy. This dilution factor is simply due to the weakening of the galactic gravitational field and increases with decreasing galaxy mass. For instance, in the least massive galaxy case the amount of retained material is diminished by at least a factor 10 with respect to the produced one. Overall, there is a potential discrepancy of more than three order of magnitude for the lightest galaxies between the produced and the retained $r$-process material. For larger galaxy masses, as of a consequence of the larger numbers of merging NSBs and of the deeper gravitational potential of the galaxy, the amount of $r$-process material increases and the discrepancy between the produced and the retained material decreases. A global dilution factor of the order of ten is still visible for the $M_{\rm b}=10^{7}~{}\rm M_{\large\odot}$ model (central panel). Finally, no significant discrepancy is visible for the most massive case (right panel), where a large number of NSB merger occur within $t=10~{}{\rm Gyr}<T_{\rm SF}$, most of them pollute the galaxy, and no significant dilution factor applies. In the following, unless differently specified, we will always refer to $r-$process mass actually retained by galaxies.

A deeper insight into the specific origin of $r-$process material is given in Fig. 5 and Fig. 6, for the CCSN to NSB ratio $x=$ 150 and 50, respectively. Upper panels show $r-$process mass generated from dynamical ejecta, both from equatorial (solid red lines) and polar components (dashed blue lines). Middle panels report instead the mass coming from disk wind ejecta, labelled as “wind 1” for the wind coming from more massive NSB that do form a BH before the disk wind emerges (solid orange lines) and “wind 2” for the ejecta coming from NSB merger that form a long-lived NS (dashed green lines). We recall that the latter cases do not produce a significant amount of Eu. Finally, lower panels show the distributions of total retained $r-$process mass (solid black lines) as well as retained europium mass (dashed cyan lines). For the dynamical ejecta, the $\sin^{2}\theta$ mass distribution, as well as the larger solid angle, favours pollution from the equatorial component, inside which europium is produced. For the disk winds, lighter NSB mergers produce more massive disc whose ejecta can pollute the galaxy more significantly, but without producing significant amount of heavy $r$-process elements, including europium. Summing up of the relevant contributions, the amount of retained europium is thus usually three orders of magnitudes smaller than the total amount of $r$-process material. Results for the more optimistic case ($x=50$, Fig. 6) are qualitatively very similar to the $x=150$ ones, but showing an higher number of successful events and larger amount of $r$-process material and europium masses. Still for the least massive case the ratio of retained vs produced europium mass remains usually very low, i.e. around $10^{-6}-10^{-5}$. Once again, larger galaxy masses reduce fluctuations in the distributions, leading to narrow histograms in the MW-like cases.

6.2 Parameter exploration

We now turn to explore the dependence of our results on some of the key parameters of the model.

In Fig. 7 we present the distribution of $r$-process material for different $x$, i.e. the ratio between CCSNe and formed NSBs. Also for the least massive galaxy the increased number of potential NSB mergers leads to an increment of the retained mass, but still revelling a high scatter, meaning that the enrichment process is dominated by small number statistic. A clear trend is instead observable in the other two cases (central and right panels of Fig. 7) showing an increase of nearly two orders of magnitude in the peak of retained mass when going from $x=1350$ (very pessimistic) to $x=50$ (optimistic).

Another crucial parameter that can affect the amount of $r-$process enrichment, especially at low galactic masses, is the IGM density. If NSBs merge outside from the galaxy, a higher IGM density determines a smaller expansion of the ejecta bubble, preventing the ejected material to fall back on the parent galaxy. Such trend is shown in Fig. 8, in which we compare our standard case with an IGM number density of $n_{0}=10^{-4}~{}{\rm cm^{-3}}$ (solid blue lines) to a situation with $n_{0}$ one hundred times higher. No sensible effect is seen in the MW case, while for Dwarf model 3 we witness a slightly decrease in the number of successful enrichment realisations and, at the same time, a decrease in the efficiency of the $r$-process enrichment. Very different is instead the case of Dwarf model 1 where assuming $n_{0}=10^{-2}~{}{\rm cm^{-3}}$ no $r-$process enrichment is verified (with 200 galaxy realisations).

In addition, for the europium mass only, we verify the dependence on our assumptions about the the $r-$process nucleosynthesis in the various ejecta components. In particular in Fig. 9 we compare our standard case, in which no europium is produced in the disk wind emerging from a remnant characterized by a long-lived massive NS, to a situation in which all disc winds produce europium, i.e. we assume the same nucleosynthesis with $80\leq A\leq 232$ for all winds. Again for the smallest galaxy case allowing this extra production of europium can increase up to one order of magnitude the quantity of Eu, but we stress that the result is heavily affected by small number fluctuations, with essentially only two realisations (out of 200) showing a significant increase in the retained europium mass. For more massive galaxies the impact is instead milder, determining at most an increase of a factor of a few.

Finally, we have tested the impact of other input parameters of the models, including the extension of the galaxy in computing the retained mass and the probability of receiving of low kick, i.e. $\xi$ and $p_{l}$. In both cases, the amount of retained material scales as expected in the case of low mass galaxies: an increase from $\xi=3$ to $\xi=5$ translates in an increase of $\sim(5/3)^{2}$ in the retained mass, due to a more extended disc surface.101010Nevertheless, we expect our standard choice ($\xi=3$, corresponding to three times the radial scale of the exponential profile) to cover a significant fraction of the galaxy in which stars are produced. For the latter parameter, a decrease in the probability from 0.7 to 0.6 produces essentially no noticeable differences in the amount of retained mass, except for the least massive galaxy case, where a slightly decrease of the captured mass arises. Any dependence on both $\xi$ and $p_{l}$ becomes less and less relevant for large galaxy masses and, in particular, for a MW-like galaxy.

7 Comparison with dwarf galaxy observations

Fig. 10 shows the amount of retained europium mass as a function of the iron mass produced by CCSNe explosions. We stress that the nucleosynthesis contribution of SNIa is neglected in our work. This assumption is valid as long as star formation happens significantly before SNIa start to explode (as in the case of UFD galaxies). Otherwise, the amount of iron computed provides a lower limit. In the figure we report the comparison between our results (small dots and empty circles) and observational data points referring to UFDs (red triangles, square for Ret II and diamond for Tuc III) and classical dwarfs (blue square), respectively. Each small dot (empty circle) represents a successful realisation in which europium is effectively captured by the host galaxy assuming a CCSN/NSB ratio of $x=150$ ($x=50$). The latter values have been chosen by comparing the amount of europium retained in MW-like galaxy models with estimate of europium mass in our Galaxy. Assuming a mass of $5000~{}\rm M_{\large\odot}$ for $r$-process material with $A\geq 90$ (see e.g. Hotokezaka et al., 2018, and references therein) and an europium mass fraction of $m_{\rm Eu}/m_{r-\rm proc,A\geq 90}=0.0042$ (Lodders, 2003), we found $50\lesssim x\lesssim 150$, depending on the detailed parameter choice. Observational data points are taken from figure 1 and table 1 of Beniamini et al. (2016b). We refer to their section 2 for the presentation of the UFD and classical dwarf samples and for the relevant references. Concerning Tucana III, we followed the same procedure as in Beniamini et al. (2016b), but from stellar abundances reported in Hansen et al. (2017).

From left to right, each group of realisations clusters around specific values of the iron mass, which is proportional to the explored baryonic masses (in the interval $10^{5}-10^{8}~{}\rm M_{\large\odot}$), and shows a small dispersion. The europium mass instead shows a much larger scatter, especially for Dwarf model 1 and 2, where variations span three orders of magnitude. On the contrary, at increasing galaxy mass the scatter decrease. Again this result is a combination of the small number of events and the geometric dilution factor. From Fig. 10 an increasing trend with galaxy mass for the europium mass is clearly visible, but when compared to the observational data points, our procedure systematically underestimate $m_{\rm Eu}$ for all galaxy models with baryonic mass in the range $10^{5}-10^{7}~{}\rm M_{\large\odot}$. To test potential systematic uncertainties in our model, in the right panel of Fig. 10, we report the same quantities of the left panel but we consider that europium production is active also in the wind emerging from long-lived remnant. This cause a shift of at most 0.5 dex for $m_{\rm Eu}$, easing the tension, but without definitively solving it. Therefore, despite the large uncertainties that affect the followed procedure, our results even in the most optimistic case (right panel of Fig. 10, $x=50$ case), might imply a possible tension with the scenario in which europium enrichment is only generated by NSB mergers.

In Fig. 11 we provide histograms of $[{\rm Eu/Fe}]$ for all our four dwarf models. Once again, results for models 1,2, and 3 are systematically smaller than the values required to explain abundances in Reticulum II and in classical dwarfs, as reported in table 1 of Beniamini et al. (2016b) by at least one dex, even in the most optimistic case $x=50$.

It is still possible to argue that the effective amount of $r$-process material produced by NSB mergers is not so well constrained and a larger production could solve the problem. However, there is at least another prediction of the model which we found more difficult to reconcile with observations. As previously pointed out, in our model dwarf galaxies show a significant spread of Eu enrichment compared to a relative fixed enrichment of iron. This spread is dependent on the galaxy mass and increasing for less massive objects. Intuitively, this is connected to the shallower potential well that the NSB encounter after the second SN kick in a less massive host galaxy. In the model, we do not follow any possible stochastic enrichment inside the galaxy due to finite dimension of the pollution by our NSMs; the spread obtained is between the average Eu in single galaxies. Therefore, the observational expectation is that classical dwarf galaxies can differ in their average Eu/Fe ratio from 1 dex for the more massive one to more than 3 dex for the lightest ones. However, according to the data up to now collected, this is not the case. Observed dwarf spheroidal galaxies with similar final stellar masses show a relative small scatter in iron (McConnachie, 2012), which is compatible with our model. On the other hand, the enrichment in Eu for classical dwarf is, in most of the cases, proportional to that of iron, differently from the outcome of our numerical modelling. For the classical dwarf, the only hint of a substantial variation of the Eu enrichment is the [Eu/Fe] ratio measured in stars of Sagittarius dwarf galaxy (McWilliam et al., 2013). The situation is different for UFD galaxies. At the moment most of them seem to have an extremely low enrichment of neutron capture elements with - at present - the only exception Reticulum II (Ji et al., 2016a) and possibly Tucana III (Hansen et al., 2017). According to our model, NSB mergers can enrich of the order of 1-2% galaxies in the mass regime of Reticulum II. Therefore, we could have simply randomly detected this object, although the chance are relatively low considering the dozen of UFD galaxies with measured stellar abundances. However, in the case the measurements of europium in four additional stars of Tucana III will be confirmed (Marshall et al., 2018), then the random probability will be certainly too low and a clear tension with the prediction would be confirmed also from this prospective. We should underline that the debate is still on whether UFD galaxies used to be isolated galaxies or are just fractions of larger tidal disrupted objects. This would relax the constrain from this side.

If the stars that we observe in UF and classical dwarfs formed at high redshift $z$, the tension with our models could potentially increase. In fact, we expect the average particle density in the Universe to increase as a function of $z$ as $(1+z)^{3}$. This implies that for dwarf galaxies forming stars within the first 3.3 Gyr after the Big-Bang (i.e. $z\sim 2$) the IGM could be 10 larger than the present Universe average and likely larger than our standard case ($n_{0}=10^{-4}{\rm cm^{-3}}$). This increase could limit the extension of the merger remnant, further reducing the dilution factor and, ultimately, the amount of matter enriching the host galaxy.

Our results depends on the chosen binary semi-major axis distribution. Neglecting eccentricity effects, a uniform distribution of $\log_{10}(a)$ translates in a $t^{-1}$ delay time distribution for NSB mergers. Beniamini & Piran (2019) recently showed that the observed distribution of Galactic NSBs might suggest the presence of a significant population (at least 40%) of fast merging ($\lesssim 1~{}{\rm Gyr}$) binaries, larger than the fraction implied by our distribution ($\sim$ 27%). A closer inspection of the NSB samples used in our models reveals that a significant increase in the enrichment fraction (at least at 10% level) in most of enriching binaries would require a much larger fraction of binaries merging not only within 1 Gry, but within a few times $10^{7}$ yr and $10^{8}$ yr in the $10^{5}~{}M_{\odot}$ and $10^{7}~{}M_{\odot}$ dwarf galaxy model, respectively.

8 Conclusions

NSB mergers are nowadays clearly recognized as one of the major source of $r$-process nucleosynthesis in the Universe and a key player in galactic chemical evolution. The event rate and the mass per event needed to explain UFD enrichment appears to be consistent with those in the Milky Way and with the NSB merger properties obtained by the analysis of GW170817 (see e.g. Hotokezaka et al., 2018). In this paper, we have explored the impact of the orbital motion of binary systems of NSs around galaxies prior to merger on the $r$-process enrichment. We have found that for low mass systems, i.e. disc galaxies with a baryonic mass $M_{\rm b}$ ranging from $10^{5}$ up to $10^{8}~{}\rm M_{\large\odot}$, the motion of the binary due to the kicks imparted by the two SN explosions determines a merger location potentially detached from the disc plane, even for gravitationally bound systems.

The immediate consequence is a dilution of the amount of $r$-process material retained by the galaxy within its star forming age $T_{\rm SF}$ (and thus potentially available for the next generation of stars). This effect is more severe for low mass disc galaxies. Assuming a log-flat distribution for the semi-major axis, realistic distributions of NSB parameters, a production rate of one double NS system every 150 CCSNe, and a rather dilute IGM (with a density of $n_{0}=10^{-4}{\rm cm^{-3}}$), in the least massive case we have explored ($M_{\rm b}=10^{5}~{}\rm M_{\large\odot}$) a galaxy has a $\sim 10\%$ probability of producing a merging NSB within ($T_{\rm SF}\lesssim 1~{}{\rm Gyr}$), and a significantly lower probability of retaining $r$-process elements from this single event ($\sim 2\%$). Since the merger happens at distances comparable or larger than the galactic disc size and outside from the disc plane, a fraction ranging between 0.1 and 0.001 of the ejected mass is actually retained (corresponding to $10^{-5}-10^{-4}~{}\rm M_{\large\odot}$) in spite of the ejection of a few $10^{-2}-10^{-3}~{}\rm M_{\large\odot}$ of $r$-process material, most of which stops and mixes with the IGM. This dilution effects is also present in more massive galaxies, but it becomes less and less relevant as $M_{\rm b}$ increases, and it has practically no relevance for a MW-like galaxy.

We have also estimated the amount of Eu retained by the galaxy. We considered that mass ejection in NSB mergers happens through different channels and each channel has a potentially different nucleosynthesis, based on the influence of neutrino irradiation. We have found that the amount of retained Eu is usually $\sim 10^{-3}$ the amount of retained $r$-process material. This is due to the presence of a significant fraction of binaries that produce a long-lived massive NS. For these remnants, discs are usually more massive and the persistent neutrino irradiation suppresses the synthesis of elements between the second and the third $r$-process peaks. The relative fraction of Eu is largely unaffected by the galaxy mass, however the paucity of enrichment events in the low mass case produces larger fluctuations.

We have tested the robustness of our results with respect to the parameters entering the model. The most relevant parameters are the fraction of double NS systems with respect to the number of CCSNe ($x$) and the IGM density ($n_{0}$). For the former, even considering the most optimistic case (i.e. one double NS system every 50 CCSNe), the results stay within a factor of a few our standard case. For the latter, a significantly larger density (which is what we expect in the early Universe) determines less extended NSB remnants and prevents $r$-process enrichment in the least massive case.

The precise modelling of the composition of NSB ejecta is still affected by uncertainties, as well as our knowledge of the fraction of systems that forms a long-lived remnant. Thus, we have repeated our analysis assuming that all disc winds produce Eu in the same ratio. We verified that our standard results are qualitatively robust, since the amount of retained Eu increases only by a factor of few. Calculations of the retained $r$-process material allows us to compare with observed abundances of iron and europium from UFD and classical dwarfs. Due to the dilution effect on the retained $r$-process material, all our realizations show a systematic deficiency in the europium abundance compared with observed abundances. Additionally, for a fixed galaxy model and especially for low mass galaxies, the large intrinsic variability introduced by the dilution process and small number statistics is not visible in the observations. Since the dilution is practically negligible in the MW-like galaxy model, in that case we have estimated the CCSN to NSB rate within our assumptions and we found $50\lesssim x\lesssim 150$, with a preference on the low value side.

Both the discrepancies on the absolute values and on their spread could potentially imply a tension in the explanation of the europium production in dwarf galaxies as produced by NSB mergers only. This conclusion could depend on observational uncertainties in the estimate of the elemental abundances or it could imply a significant change in one or more of the canonical assumed parameters. However, to reconcile our models with observations a systematically optimistic choice of parameters is required, together with a significant revision of the observed elemental abundances.

On the one hand, a possible solution to the discrepancy is to assume the existence of a fast merger population (e.g. Beniamini & Piran, 2019), such that a significant fraction of NSBs merge within $10^{8}$ yr, or even a few times $10^{7}$ yr. Under these conditions the merger does not happen too far from the galaxy and the dilution factor becomes $\gtrsim 0.1$. If, on the other hand, this tension is confirmed, the existence of additional sites for the production of $r$-process elements is necessary. A similar conclusion was also obtained by other authors and motivated by the difficulty to inject r-process elements early enough to explain the Eu abundances in metal-poor stars (Matteucci et al., 2014; Cescutti et al., 2015b; Wehmeyer et al., 2015; Haynes & Kobayashi, 2019). Moreover, also Côté et al. (2019) and Simonetti et al. (2019) further suggested an extra production site of europium as a possible way of reproducing the decreasing trend of [Eu/Fe] in the Galactic disk, a different feature that is the result of $\sim 12$ Gyr of chemical evolution.

BH-NS mergers are also possible sources of $r$-process elements (e.g. Shibata & Taniguchi, 2011). The lack of a MNS in the remnant favours the production of the heaviest $r$-process elements (e.g. Roberts et al., 2017; Lippuner et al., 2017). This could enhance the amount of retained europium by a factor of a few. However, even if BH-NS merger rates are highly unconstrained, we expect them to be significantly lower than NSB merger rates (Abbott et al., 2018). A possible solution would thus require qualitatively different binary parameters, for example much lower kick velocities or much smaller initial separations, resulting from a possibly different binary evolution. This could however be in tension with population synthesis results, considering that a significant amount of ejecta from BH-NS merger requires a not too large mass ratio ($\lesssim 5$), for moderately high BH spins.

An alternative solution is represented by special classes of CCSNe, able to produce $r$-process elements with an amount comparable to the one of NSB mergers and with a similar rate. Possible examples include magnetically-driven CCSNe (Fujimoto et al., 2008; Winteler et al., 2012; Nishimura et al., 2015; Mösta et al., 2018) and disk ejecta from collapsar models (e.g. Malkus et al., 2012; Siegel et al., 2019). In these cases the ejection of $r$-process material happens still on a sufficiently short time scale and inside the galaxy, such that no significant dilution factor affects the enrichment (Beniamini et al., 2018). Magnetically-driven CCSNe have also the advantage that they could also explain the chemical enrichment of Galactic halo (Cescutti & Chiappini, 2014). Further multi-physics studies combining galactic (chemical) evolution, binary population synthesis, NSB mergers, and taking into account nuclear uncertainties are required to address all the open issues in the field. In this respect, a crucial role is represented by the forthcoming determination of more precise compact binary merger rates by the Advanced LIGO and Virgo detectors (Aasi et al., 2015; Acernese et al., 2015; Abbott et al., 2018).

Acknowledgements

MB, AP and MD acknowledge CINECA, under the TEONGRAV initiative, for the availability of high performance computing resources and support. G. Cescutti acknowledges financial support from the European Union Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No. 664931.