The chemical evolution of r-process elements from neutron star mergers: the role of a 2-phase interstellar medium

TL;DR

This paper demonstrates that a two-phase interstellar medium model, with specific cooling times and yield injection ratios, successfully explains observed europium trends from neutron star mergers in the Milky Way.

Contribution

It introduces a two-phase ISM model with differential yield injection to better match r-process element observations, resolving previous discrepancies in chemical evolution models.

Findings

Two-phase ISM model reproduces observed Eu trends.

A merger timescale of about 150 Myr fits data.

Europium likely originates from neutron star mergers or similar sources.

Abstract

Neutron star mergers (NM) are a plausible source of heavy r-process elements such as Europium, but previous chemical evolution models have either failed to reproduce the observed Europium trends for Milky Way thick disc stars (with [Fe/H] ~ -1) or have done so only by adopting unrealistically short merger timescales. Using analytic arguments and numerical simulations, we demonstrate that models with a single-phase interstellar medium (ISM) and metallicity-independent yields cannot reproduce observations showing [Eu/alpha] > 0 or [Eu/Fe] > [alpha/Fe] for alpha-elements such as Mg and Si. However, this problem is easily resolved if we allow for a 2-phase ISM, with hot-phase cooling times \tau_{cool} of order 1 Gyr and a larger fraction of NM yields injected directly into the cold star-forming phase relative to alpha-element yields from core collapse supernovae (ccSNe). We find good agreement with observations in models with a cold phase injection ratio f_{c,NM}/f_{c,ccSN} of order 2, and a characteristic merger timescale \tau_NM=150 Myr. We show that the observed super-solar [Eu/alpha] at intermediate metallicities implies that a significant fraction of Eu originates from NM or another source besides ccSNe, and that these non-ccSN yields are preferentially deposited in the star-forming phase of the ISM at early times.

Equations

$${[Mg/Fe]}_{\mathrm{pl}}=\log \frac{y_{\mathrm{Mg},\mathrm{ccSN}}}{y_{\mathrm{Fe},\mathrm{ccSN}}}-\log \frac{M_{\mathrm{Mg},\odot }}{M_{\mathrm{Fe},\odot }}.$$

$${[Mg/Fe]}_{\mathrm{pl}}\sim \log \left(\frac{y_{\mathrm{Fe},\mathrm{ccSN}}+y_{\mathrm{Fe},\mathrm{SNIa}}}{y_{\mathrm{Fe},\mathrm{SNIa}}}\right)=\log \left(1+\frac{y_{\mathrm{Fe},\mathrm{SNIa}}}{y_{\mathrm{Fe},\mathrm{ccSN}}}\right),$$

$$\hbox{[Eu/Fe]}_{\rm pl}\leq\hbox{[Mg/Fe]}_{\rm pl}$,$$

$

$$[Eu/Fe]\le [Mg/Fe].$$

$${[Eu/Si]}_{\mathrm{pl}}\le [\mathrm{Eu}/{\mathrm{Si}}_{\mathrm{ccSN}}]+0.1$$

$$\hbox{[Eu/Fe]}_{\rm pl}\leq\hbox{[Si/Fe]}_{\rm pl}+0.1$.$$

$

$${[Eu/Mg]}_{\mathrm{pl}}\le {\log }_{10}\frac{{\mathrm{f}}_{\mathrm{c},\mathrm{NM}}}{{\mathrm{f}}_{\mathrm{c},\mathrm{ccSN}}}$$

$${[Eu/Fe]}_{\mathrm{pl}}\le {[Mg/Fe]}_{\mathrm{pl}}+{\log }_{10}\frac{{\mathrm{f}}_{\mathrm{c},\mathrm{NM}}}{{\mathrm{f}}_{\mathrm{c},\mathrm{ccSN}}}.$$

$$Z_{\mathrm{inst}}(t)=\frac{y_{\mathrm{inst}}}{1+\eta -r}\left[1-e^{-t/{\tau }_{\mathrm{d}}}\right]$$

$$Z_{\mathrm{del}}(t)=\frac{y_{\mathrm{del}}}{1+\eta -r}\left[1-e^{-\frac{\Delta t}{{\tau }_{\mathrm{d}}}}-\frac{{\tau }_{\mathrm{del}}}{{\tau }_{\mathrm{del}}-{\tau }_{\mathrm{d}}}(e^{-\frac{\Delta t}{{\tau }_{\mathrm{del}}}}-e^{-\frac{\Delta t}{{\tau }_{\mathrm{d}}}})\right],$$

$${Z_{\rm Eu}(t)\over Z_{\rm Mg}(t)}\approx\left({y_{\rm Eu}\over y_{\rm Mg}}\right){\rm f_{\rm c,NM}F(t,\tau_{\rm d})+(1-\rm f_{\rm c,NM})F(t,\tau_{\rm cool})\over\rm f_{\rm c,ccSN}F(t,\tau_{\rm d})+(1-\rm f_{\rm c,ccSN})F(t,\tau_{\rm cool})}$,$$

$

$$F(t,\tau)=1-e^{-t/\tau}$.$$

$

$$\frac{Z_{\mathrm{Eu}}(t)}{Z_{\mathrm{Mg}}(t)}\approx \frac{1}{2}\left(\frac{y_{\mathrm{Eu}}}{y_{\mathrm{Mg}}}\right)\left(\frac{{\mathrm{f}}_{\mathrm{c},\mathrm{NM}}}{{\mathrm{f}}_{\mathrm{c},\mathrm{ccSN}}}\right)\left(\frac{t}{{\tau }_{\mathrm{NM}}}\right)\left(1-{\left(\frac{t_{\min ,\mathrm{NM}}}{t}\right)}^2\right)$$

Full text

††pubyear: 2019††pagerange: The chemical evolution of r-process elements from neutron star mergers: the role of a 2-phase interstellar medium–A2 Further model tests

1 Introduction

With the latest direct observation of a neutron star merger, as well as previous studies from Swift GRB counts (Wanderman & Piran, 2015), it has become fully evident that neutron star mergers (NM) are numerous enough and have a sufficient heavy element yield (Smartt et al., 2017) to play a significant role in shaping the chemical evolution of r-process elements like Europium (Eu). However, there has been a long tradition of chemical evolution papers challenging the impact of neutron star mergers (see below). The central issue is the delay time distribution (DTD), i.e. the average distribution in time of events, such as Type Ia supernovae (SNeIa) or NM, after the birth of a stellar population. The delay time distribution is best characterised by an onset time ($t_{\rm min}$, short and not so important for neutron star mergers) and a characteristic timescale, e.g. the decay time $\tau$ of an exponential. Recent studies have approached the chemical evolution problem with different modelling strategies (see further below), but agree widely on the outcome: they have either assumed reasonable DTDs for neutron star mergers and concluded that they cannot explain measured abundance patterns (Argast et al., 2004; Wanajo & Ishimaru, 2006), or claimed that NM work while adopting implausibly short characteristic timescales of order $1\,{\rm Myr}$ (Matteucci et al., 2014), which are not feasible in population synthesis.

In this paper we will show that neither argument applies. We will demonstrate that classical chemical evolution models that assume a 1-phase interstellar medium (ISM) and metallicity-independent yields cannot possibly reproduce the observed stellar abundances. Specifically, such models will always predict too low europium/r-process abundances of metal-poor stars relative to those at solar metallicity. Conversely, our models that account for differential injection into the warm/hot gas phase reproduce the high [Eu/Fe] ratios observed for stars with $\hbox{[Fe/H]}\sim-1$. Our arguments will further show that a significant fraction of r-process elements must come from a different origin than core-collapse supernovae (ccSNe).

Two concise and informative reviews of the problem of neutron star mergers and the chemical evolution of r-process elements can be found in Hotokezaka et al. (2018) and Cowan et al. (2019). Evidence from the observed bolometric luminosity evolution confirms a low-$Y_{e}$ (electron to nucleon ratio) nucleosynthesis with sufficient heavy element production, which accounts for the origin of Lanthanides (Rosswog et al., 2018). Ever since the discovery of neutron stars, it has been predicted (see e.g. Lattimer & Schramm, 1974) that r-process elements can be generated and injected to the star-forming ISM by mergers of a neutron star with another neutron star or a black hole, when parts of the neutron star decompress rapidly and get expelled either in a tidal tail (likely favouring heavier r-process nuclei, and in excellent agreement with the solar elemental composition, Freiburghaus et al., 1999), or in a polar outflow (likely favouring less massive nuclei, Perego et al., 2014; Martin et al., 2015). Thielemann et al. (2017) show that the abundance pattern of r-process elements can be fully modelled by the nucleosynthesis in decompressing neutron star matter during mergers. The alternative site for r-process element generation would be pockets in ccSNe which have a large free neutron flux and low electron fraction (Witti et al., 1994; Takahashi et al., 1994). The most promising ccSN site are collapsar models, where jets/outflows from the black hole accretion disks formed in magneto-rotational core collapse hit the stellar material (Bisnovatyi-Kogan, 1970; LeBlanc & Wilson, 1970; Nishimura et al., 2017; Halevi & Mösta, 2018). In a recent discussion, Siegel et al. (2019) argue that these could be the dominant r-process production site, as their direct association with ccSNe guarantees a minimal time delay.

On the observational side, Guiglion et al. (2018) show data that hint at slightly different trends of the r-process elements [Gd/Fe] and [Dy/Fe] vs. [Fe/H] than [Eu/Fe], which would argue for some difference in their nucleosynthesis sources, but it remains to be seen if this trend remains robust when accounting for NLTE or 3D effects (departures from the assumed local thermodynamic equilibrium or the assumptions of 1-D atmospheres vs. a full 3D atmosphere model causing bias in the abundances inferred from stellar spectra). While there were doubts if the rate of neutron star merger events would suffice to supply sufficient r-process elements, this is no concern any more after the recent observation of the GW170817 event by LIGO (Hotokezaka et al., 2018; Côté et al., 2018). What remains unsolved is the central problem of the detailed enrichment history and the resulting abundance patterns.

Galactic chemical evolution models can assess if r-process elements predominantly originate in special environments in ccSNe or in neutron star mergers by exploiting their different timescales: If r-process elements form in ccSNe, there should be no difference in enrichment trends between r-process elements and alpha elements (e.g. O, Mg, Si), while r-process enrichment by neutron star mergers will show the signatures of significant time delay, i.e. under-abundant r-process elements relative to alpha elements for stars formed early in the history of the system. We will discuss the relevant DTD timescales in Section 3. Chemical evolution studies usually use one of the following approaches: i) modelling very metal-poor stars with stochastic evolution of various flavours that model many separate domains or allow for local deviations (e.g. Argast et al., 2004; Cescutti & Chiappini, 2014; Wehmeyer et al., 2015), an approach that is vital if one attempts to fit these stars ($\hbox{[Fe/H]}\lesssim-2$), ii) running a classical chemical evolution model (e.g. Matteucci et al., 2014), i.e. one phase of ISM and one or multiple zones of the Galactic disc(s), or iii) inserting chemical evolution recipes into hydro simulations, like van de Voort et al. (2015) or Haynes & Kobayashi (2019).

Matteucci et al. (2014) formulated successful r-process enrichment models only by assuming unfeasibly short delay time distributions, with characteristic delay timescales of only $1\,{\rm Myr}$. Other models, like Cescutti & Chiappini (2014) and Cescutti et al. (2015), have been quite successful in modelling halo stars with stochastic chemical evolution, but again with rather short characteristic delay timescales of $\tau_{{\rm NM}}\sim 10\,{\rm Myr}$. Successful models of r-process abundances in halo stars (Hirai et al., 2015) have been restricted to $\hbox{[Fe/H]}<-1$, and the continuation to higher metallicities appears problematic. On the other side of the same medal, recipes in hydro simulations, like van de Voort et al. (2015), get the mean [Eu/Fe] abundance wrong ($\hbox{[Eu/Fe]}<0$ instead of $\hbox{[Eu/Fe]}>0$) along the entire metal-poor track. The real challenge (as will be shown in our discussion below) is to explain high r-process abundances in thick disc stars, while still reaching solar values at later stages, a point also noted by Côté et al. (2019). One way out could be to go for metallicity-dependent yields, setting much larger r-process yields at low metallicity (see e.g. Mennekens & Vanbeveren, 2016), but models of this sort fit the data even worse. Solving the problem this way would also imply very high [Eu/Fe] values and anomalously high [Eu/$\alpha$] in dwarf galaxies, in tension with observations (Geisler et al., 2005; Lanfranchi et al., 2008), and it would incline the [Eu/Fe] vs. [Fe/H] plateau of the thick disc downwards (higher [Eu/Fe] at lower [Fe/H]), again in tension with the data. The collapsar model (Siegel et al., 2019) would yield a significant change of r-process element production near solar metallicity and merits further investigation of its chemical evolution predictions, in particular a comparison with the current $\hbox{[Eu/Fe]}(R)$ profile in the Milky Way.

We will argue in this paper that the essential ingredient for resolving this problem is proper accounting for enrichment into different phases of the ISM. While discussions of the ISM often distinguish three or more phases – e.g., cold, warm, and hot – for our purposes we need only distinguish (cold) gas that can immediately form stars from (hot) gas that must first cool on a $\sim\,{\rm Gyr}$ timescale before entering the star-forming phase.

Our paper is structured as follows: we start with a general discussion of the basic picture and chemical evolution models in Sections 2 and 3, followed by a description of our adopted chemical evolution model in Section 4, detailing the use of different gas phases. In Section 5, we present some general analytic arguments and constraints for chemical evolution. The models to solve the r-process abundance problem are step by step developed from a more classical chemical evolution picture in Section 6. In Section 7 we summarize the conclusions from this analysis. Appendix A1 presents an analytic model, based on the formalism of Weinberg et al. (2017), that reproduces the findings of Section 6. Appendix A2 demonstrates the insensitivity of our models to the assumed SNIa DTD and varies $\tau_{\rm NM}$ for a 2-phase model.

2 The abundance plane

The purpose of this Section is to give a quick overview of how to interpret the structures in abundance planes, like the [$\alpha$/Fe]-[Fe/H] plane, and to introduce some of the terms used throughout this paper. Fig. 1 1 displays data in the $\hbox{[Mg/Fe]}-\hbox{[Fe/H]}$ and $\hbox{[Eu/Fe]}-\hbox{[Fe/H]}$ abundance planes taken from the SAGA database (Suda et al., 2008, 2017) and the Bensby & Battistini sample (Bensby et al., 2014; Battistini & Bensby, 2015, 2016). The experienced reader may wish to inspect our definition of domains in Fig. 1 and proceed to the next section.

Understanding the abundance planes is our only route to explore the sources of early chemical evolution. Testing characteristic delay times of order $100\,{\rm Myr}$ with stellar ages, would demand systematics-proof age measurements with accuracies and precision better than $\sim 100\,{\rm Myr}$ for stars that are more than $10\,{\rm Gyr}$ old; a hopeless prospect at current stage, even with the accuracy of Gaia parallaxes. However, relative abundances naturally resolve this and serve as a finely resolved clock at early times: changes in relative abundances happen when new progenitor classes with different DTDs start to contribute. However, note that metallicity is usually a poor indicator of a star’s age. Accreted stars from dwarf galaxies with their own clocks will overlap with low-metallicity Milky Way stars, which again have at each time a position-dependent metallicity.

Metal-poor stars are dominated by stochastic chemical evolution, i.e. Poisson noise due to a low number of contributing ccSNe, SNeIa and NM, as well as direct stochasticity, e.g. by the question where and how close a star-forming cloud is to such an event. This stochasticity drops rapidly towards larger [Fe/H] (we remind that an increase by $1\,{\rm dex}$ in [Fe/H] equals a factor $10$ in its abundance), and becomes negligible near $\hbox{[Fe/H]}\gtrsim-2\,{\rm dex}$.

The features that we are most interested in are the values $\hbox{[Mg/Fe]}_{\rm pl}$ and $\hbox{[Eu/Fe]}_{\rm pl}$ of the high-[Mg/Fe] ridge, which is commonly the chemically defined thick disc of the Milky Way and at low [Fe/H] shares its location with old halo stars. This ridge is particularly well defined in the more precise Bensby et al. data (blue points), as well as the location of the thick disc knee, which marks the set-in of SNeIa in the inner Galactic disc, and which is nearly guaranteed to be free of any contamination by halo stars or stochastic chemical evolution effects. Lower [Mg/Fe] values imply that SNeIa have contributed a significant fraction of the iron budget (i.e. at least a few Gyrs after the formation of the system) and thus designate younger objects: the thin disc on the high [Fe/H] end separated from the younger halo star populations at $\hbox{[Fe/H]}\lesssim-0.7\,{\rm dex}$. We note that this [$\alpha$/Fe] clock is again not an absolute time - systems have different star formation histories, and [$\alpha$/Fe] values may depend on galactocentric radius even in the Milky Way.

3 Basic chemical evolution: sourcing the elements

To understand the chemical evolution processes shaping the distribution of stars in abundance space, it is useful to draw the gas balances. Fig. 2 shows a graph of gas balances in a 1-phase (top) vs. a 2-phase model accounting for a hot ISM (bottom). To understand the problem of r-process abundances, we only need to account for the three progenitor classes that produce iron (Fe), alpha elements, and r-process elements: ccSNe, SNeIa, and NM. Each of these classes will, with its own delay time distribution, yield its own mix of elements back to the ISM, from where new stars are formed.

The functional shape of DTDs is frequently taken to be $t^{-1}$ for SNeIa; this shape is suggested by population synthesis models compared with data (Hachisu et al., 2008; Totani et al., 2008). Similarly, the DTDs for neutron star mergers is frequently described by a $t^{-1}$ or $t^{-1.5}$ law. From a practical point of view, this functional shape is unattractive, since it is not integrable to at least one side, thus requiring a cut-off and normalisation conditions sensitive to this cut-off. Here, we will adopt a simple exponential: ${\rm DTD}(t)=\exp(-t/\tau)$, where $\tau$ is the characteristic timescale. In Appendix A2 we show that using a double exponential law for SNeIa, which fits the power law more closely, does not affect our results. The difference for NM is even less significant due to the shorter timescale.

Let us outline these three progenitor classes:

• •

Massive stars exploding (usually as ccSNe): most massive stars attain an onion-like structure just prior to collapse, with an Fe/Ni core at the centre, surrounded by shells dominated by $\alpha$ elements of decreasing mass number. In the explosion, most of the shells are expelled into the ISM, providing the bulk of today’s inventory of $\alpha$ elements, like Si, or Mg, or Ca. The mass cut where material can still escape the emerging neutron star determines the early production of Fe. Theoretical models of ccSNe do not predict the mass cut robustly. Some empirical indications come from estimates of ${}^{56}{\rm Ni}$ mass in ccSNe. The picture of a simple mass cut is only a naive approximation to real explosions, where convection and anisotropy (evidenced by neutron star kicks and remnant observations) allow for some expulsion of material below the average mass cut. However, the average ${\rm Fe}$ yield from ccSNe is constrained by the \hbox{[\alpha/Fe]}_{\rm pl} value of the thick disc (see Fig. 1). The value of \hbox{[\alpha/Fe]}_{\rm pl} is close to $\log_{10}{2}=0.3$, which shows that ccSNe account for roughly half of today’s ${\rm Fe}$ inventory. The massive progenitors of ccSNe ($M\gtrsim 10\,{\rm M}_{\odot}$) explode after lifetimes of a few Myrs to less than $100\,{\rm Myr}$.

• •

Exploding white dwarfs (SNeIa): we know that the remaining half of today’s ${\rm Fe}$ must be produced in SNIa explosions. There is still considerable debate about which systems account for which fraction of SNeIa, e.g. accretion onto a white dwarf to reach the Chandrasekhar mass ${\rm M_{ch}}$, or triggered explosions below ${\rm M_{ch}}$ by a pressure wave from surface He burning, or white dwarf mergers in a double degenerate system. The common denominator is the explosion of a C/O white dwarf, which will burn about half of its mass directly to ${}^{56}{\rm Ni}$ decaying into Fe. We use a characteristic delay timescale of $\tau_{\rm SNIa}=1.5\,{\rm Gyr}$.

• •

The source of r-process elements: either neutron star mergers or pockets with particular conditions (jets, neutrino driven winds, collapsar scenario) in ccSNe. In the latter case, the delay time distribution is nearly identical with ccSNe. For NM the DTD timescale should be faster than that of SNeIa, since the progenitor systems are more massive, arise on the ccSN timescale, and at least the single degenerate scenario for SNeIa depends on the evolution of the secondary star. Also, the Swift gamma ray burst counts show that $\sim 40\%$ of the neutron star mergers have no measurable delay relative to star formation (Wanderman & Piran, 2015; Fong et al., 2017). The frequently used $t^{-1}$ or $t^{-1.5}$ DTD (Hotokezaka et al., 2018; Côté et al., 2018) combined with their shorter minimum delay time relative to SNeIa implies an equivalent characteristic timescale $\tau_{{\rm NM}}$ that is also about an order of magnitude shorter than for SNeIa. We therefore take $150\,{\rm Myr}$ as a best guess for $\tau_{{\rm NM}}$, but we will investigate a range of values.

Fig. 3 shows the tracks of cold ISM abundance through the [Mg/Fe] vs. [Fe/H] abundance plane. All models (described further in Section 4) start near metal-free, and move from left to right. The initial abundance ratios (y-axis) on the left-hand side are always set by the fastest channel of enrichment. At $\hbox{[Fe/H]}>-2.5\,{\rm dex}$, they are set by all sources that contribute on a timescale of about $100\,{\rm Myr}$. The time spacing between model points at each of the three plotted galactocentric radii is $15\,{\rm Myr}$. The onset of SNIa explosions begins to drive down [Eu/Fe] and [Mg/Fe] around $\hbox{[Fe/H]}\gtrsim-1\,{\rm dex}$. We will later see that accounting for a hot phase, which can lock up yields for of order $1\,{\rm Gyr}$, can significantly alter the effective DTD for yields from the two faster progenitor classes (ccSNe, NM), or even reverse their temporal order, if their cold channel fractions $f_{c,i}$ differ.

4 Galaxy model

4.1 General description

For this paper we use the Schönrich & Binney (2009) model with the parameters of Schönrich & McMillan (2017). A full description is provided in those two papers, so we will limit ourselves to just some basic facts. Similar to classical chemical evolution models Chiappini et al. (1997) the model divides the galactic disc into concentric rings, with a spacing of $\sim 0.2\,{\rm kpc}$, using a timestep of $15\,{\rm Myr}$. We experimented with higher resolutions, but found that they do not matter for modelling chemical evolution in general, as long as we stay out of the extremely metal-poor regime, where results are anyway dominated by stochastic chemical evolution. The model fully accounts for time-dependent yields from stars, radial migration in both stars and gas, and stellar orbital motions (which further broaden the distribution in stellar radii beyond the radial migration). The model has an inflow of fresh gas to sustain star formation (Chiosi, 1980), which by angular momentum conservation drives radial flows (Lacey & Fall, 1985) into the disc according to the prescription of Bilitewski & Schönrich (2012). The relative angular momentum of infall (about $0.75$ of the disc angular momentum) is fitted to match the radial abundance profile of the Galactic disc as measured with Cepheids.

Yields for Mg, Si, and Fe in the ccSNe are taken from Chieffi & Limongi (2004), and we use the SNIa yields for all three elements from Iwamoto et al. (1999). The results are compared to the solar abundance scale from Grevesse et al. (2007). The solar abundances for Mg and Fe are well-measured. One can argue a little about the yields, but ultimately this comes down to fixing the iron yield from ccSNe to fit the thick disc plateau and then adapting the SNIa rates to reach solar abundances; thus, using different sources would not significantly affect our conclusions. The yields for r-process elements are still very uncertain, and attempts to compile such data would only distract from our main aim to investigate the range of possible delay-time distributions. Thus, we have implemented r-process yields as a generic element and normalise them such that at late times the standard model achieves $\hbox{[r/Fe]}=0$.

An open question is the loss rate of metals to the IGM. If one assumes a Salpeter IMF, only minimal losses of order $10\%$ are needed to reach near-solar metallicity in the local thin disc. However, the now favoured Chabrier or Chabrier-like (e.g. Kroupa) IMF has far less low-mass dwarf stars to lock up metals, and so implies a higher loss-rate of $\sim 0.6$ (only $40\%$ of metals are retained). Yields from stellar populations and thus required loss rates could also be lower if many massive stars collapse directly into black holes without a significant explosion (Adams et al., 2017).

In the numerical chemical evolution models applied here, metal loss is treated by a direct loss from the galaxy without evaporating a significant amount of the ambient cold ISM. The opposite extreme is ejection of gas at the ISM metallicity characterised by a mass loading parameter $\eta=\dot{M}_{\rm out}/\dot{M}_{\star}$. In the prescription of Andrews et al. (2017), which does not implement radial flows and thus needs a higher loss rate, we require $\eta=2.5$. We show in Appendix A1 that the choice of prescription affects the chemical evolution of [Eu/Fe] and [Fe/H] at very low metallicity, $\hbox{[Fe/H]}\lesssim-2\,{\rm dex}$, but does not affect our conclusions based on the thick disc metallicity regime. The SNIa rates in the model have been adjusted to match the [Mg/Fe] plane, and are consistent with expectations for the Milky Way. As discussed in §5, we believe that [Mg/Fe] ratios may be artificially boosted by NLTE effects, making Si a more reliable reference element in measurements. However, our SNIa rates are calibrated to match observed [Mg/Fe] rather than [Si/Fe], and rather than readjust the rates when comparing to Si observations, we simply reduce model ${\rm Si}$ values by $0.07\,{\rm dex}$ (and thus decrease [Si/Fe] and increase [Eu/Si] by the same amount).

The delay time distribution for SNeIa is set as a single exponential with timescale $\tau_{\rm SNIa}=1.5\,{\rm Gyr}$, with a minimum delay time of $150\,{\rm Myr}$. We show in Appendix A2 that the precise formulation of this law does not affect our results significantly. Analogously, we use a minimum time delay of $15\,{\rm Myr}$ and several choices of timescale $\tau_{{\rm NM}}$ for neutron star mergers.

4.2 Distribution of elements to different gas phases

The main benefit of using the Schönrich & Binney (2009) model is that it accounts with a simple 2-phase ISM for the presence of hot gas, which cannot directly partake in the star formation. But which elements and which yields go where? One could be tempted to think that all yields are first locked up in the hot gas phase, and then cool back into the star-forming ISM. On the timescales we are interested in, however, a few observations show that this is not the case; for example, the short-lived isotopes identified in Earth’s crust and solar system meteorites Wasserburg et al. (2006) prove that there must be a direct channel from stellar yields to the cold star-forming ISM. Further, we know for sure that at least the ejecta from novae (Amari et al., 2001) and even from ccSNe (e.g. Travaglio et al., 1999) can directly form dust grains, which made their way directly to the meteorite formation sites of the nascent solar system.

Due to their long delay time, SNeIa will commonly explode far from the star forming regions, so that their ejecta are not trapped by a dense surrounding ISM, and they are also more energetic than ccSNe. Neutron star mergers will be different from ccSNe, and are likely affected by the different timing: the early events should still take place predominantly in dense gas regions, but the later ones might not. Most importantly, NM produce relatively dense ejecta with very heavy chemical elements. This matters, since a “normal” ISM can achieve the well-known metastable hot state, where a few million Kelvin hot plamsa is fully ionised and thus has much lower cooling rates than at, say, $10^{5}\,{\rm K}$. However, the ionisation temperature goes with the first ionisation energy, which in turn depends on $Z^{2}$. So, while hydrogen has a typical ionisation temperature of $10^{4}\,{\rm K}$, this temperature is of order $10^{8}\,{\rm K}$ for uranium, i.e. there is basically no astrophysical system where heavy r-process elements are fully ionised, not even when the near-relativistic ejecta shock with the surrounding ISM. Therefore, we expect a larger fraction of NM yields to cool quickly.

Given the above, we can make the following assumptions

• •

We put the bulk of SNIa yields ($0.99$) into the hot gas phase, i.e. ${\rm f_{\rm c,SNIa}}=0.01$. The choice hardly affects chemical evolution (see below).

• •

About $\rm f_{\rm c,ccSN}=0.25$ of ccSN yields is fed to the cold gas phase. We test different values.

• •

The fraction of neutron star merger material directly fed the cold phase is set as $\rm f_{\rm c,NM}=\rm f_{\rm c,ccSN}$ for our initial models and later varied independently.

We note that ${\rm f_{\rm c,SNIa}}$ is not important to our discussion, since the cooling timescale is still even a bit shorter than the characteristic SNIa delay time, i.e. increasing ${\rm f_{\rm c,SNIa}}$ would not significantly alter our results. The cooling timescale (again exponential) is set to $1\,{\rm Gyr}$ as in SB09, in concordance with the mass of warm/hot ISM found in the Milky Way.

5 Some analytical considerations about abundance levels

The key observational results that will drive our conclusions in this paper are that disk stars in the metallicity range $\hbox{[Fe/H]}\sim-1.5$ to $-0.5$ have typical $\hbox{[Eu/Mg]}>0$, typical $\hbox{[Eu/Si]}\sim 0.1-0.2$, and typical $\hbox{[Eu/Fe]}\ga\hbox{[Mg/Fe]}$ (see Figures 3 and 4 below).

Before proceeding to detailed comparison of numerical models to observations, it is useful to derive some basic expectations and implications from analytic arguments. Here we focus on limiting cases, with Appendix A1 presenting a more complete analytic model that confirms these limits and supports our numerical findings.

We begin by making the conventional assumption of a 1-phase, well mixed ISM at all times. If the IMF averaged ccSN yields of Mg and Fe are independent of metallicity, then at early times before SNIa enrichment is important the ratio of these elements in the ISM must simply equal the ratio of their ccSN yields, implying a plateau ratio

[TABLE]

The flatness of the observed [Mg/Fe] plateau supports the assumption of metallicity independent yields, which is also consistent with theoretical expectations for $\alpha$ elements. At later times, SNIa enrichment adds Fe without associated Mg, thus depressing [Mg/Fe]. The solar ratio $\hbox{[Mg/Fe]}=0$ reflects this latter condition, with roughly equal Fe contributions from ccSNe and SNeIa. As a rough approximation, we expect

[TABLE]

though in detail the drop from the plateau to the late-time equilibrium ratio depends on the star formation history (see, e.g., figure 3 of Weinberg et al., 2017).

Now let us assume that the $r$-process yield from NM is also independent of metallicity, and that the time delay for NM is much shorter than for SNeIa. In this case the Eu/Mg rises quickly from zero to a ratio that reflects the ratio of yields, from NM and ccSNe, respectively. Provided the yields remain metallicity independent up to solar [Fe/H], this ratio must correspond to the solar value $\hbox{[Eu/Mg]}=0$. Since $\hbox{[Eu/Fe]}=\hbox{[Eu/Mg]}+\hbox{[Mg/Fe]}$, we have

[TABLE]

which is an inequality because Eu production always lags at least slightly behind Mg production. If the NM delay time is not sufficiently short, then SNIa enrichment can depress [Mg/Fe] before [Eu/Mg] rises to the yield ratio, in which case there is not a flat plateau of [Eu/Mg], but at all times

[TABLE]

The above inequalities approach equality for very short NM characteristic delay times. If the dominant source of Eu is ccSNe rather than NM, and yields are metallicity independent, then we simply have $\hbox{[Eu/Mg]}=0$ and $\hbox{[Eu/Fe]}=\hbox{[Mg/Fe]}$ at all times. Yields of alpha elements from SNeIa allow for a small violation of the above inequalities (see below).

These arguments tell us that a 1-phase model with metallicity independent yields cannot explain observations showing $\hbox{[Eu/Mg]}_{\rm pl}>0$ or $\hbox{[Eu/Fe]}_{\rm pl}>\hbox{[Mg/Fe]}_{\rm pl}$, regardless of whether ccSNe or NM are the main source of Eu.

Spectroscopic measurements of [Mg/Fe] are subject to significant non-LTE corrections (Bergemann et al., 2017), which become stronger towards lower metallicity. We thus feel compelled to compare to Si, as a second alpha element, which should have different biases. However, Si likely has a much stronger yield contribution from SNeIa. This effect is anticipated in theoretical yield models from Iwamoto et al. (1999), which have an ${\rm Fe}/{\rm Mg}$ compared to the Sun that is an order of magnitude larger than for ${\rm Fe}/{\rm Si}$. The resulting change between ${\rm Mg}$ and ${\rm Si}$ in the SM17 chemical evolution model still underpredicts what is inferred from empirically calibrated yield models (Rybizki et al., 2017) and from [Si/Mg] ratios in APOGEE (Weinberg et al., 2019), which both indicate that $\sim 20\%$ of the solar ${\rm Si}$ content are contributed by SNeIa. A contribution of this magnitude allows a violation of the above inequalities by a factor $1.25$ or about $0.1\,{\rm dex}$, i.e.

[TABLE]

and

[TABLE]

For the yields in our numerical model, the allowed violation is only about $0.05\,{\rm dex}$.

Things change if we allow for a hot ISM phase and assume that only a fraction $\rm f_{\rm c,ccSN}$ of ccSN yields and $\rm f_{\rm c,NM}$ of NM yields go directly into the cold, star-forming ISM. Before there is time for cooling from the hot phase, this model has the same effect as reducing yields by a factor $f_{c}$. However, these factors do not enter the normalization to solar abundance ratios provided the ISM cooling timescale is short relative to the time needed to evolve to solar [Fe/H]. Our previous results therefore adjust to

[TABLE]

and

[TABLE]

For Si in place of Mg, the same inequalities hold but with $+0.1$ on the right hand sides. We thus see that the observational constraints can potentially be satisfied if $\rm f_{\rm c,NM}>\rm f_{\rm c,ccSN}$, as we have argued is plausible in Section 4.2. For example, with $\rm f_{\rm c,ccSN}=0.25$ and $\rm f_{\rm c,NM}=0.5$, values of [Eu/Fe] can rise $\sim 0.3\,{\rm dex}$ over the [Mg/Fe] plateau or $\sim 0.4\,{\rm dex}$ over the [Si/Fe] plateau.

One could try to employ 1-phase models to explain the observed Eu evolution by invoking enhanced production at low metallicity, thus boosting [Eu/Mg] and [Eu/Fe] in the plateau regime relative to solar abundance. Metallicity dependence of NM yields could arise from changing binary populations or neutron star production rates, as suggested by Mennekens & Vanbeveren (2016). However, their resulting chemical evolution models do not match the observations they compared to. In general, such a model would at least require fine-tuning, and potentially impossible parameter choices; it naturally predicts enhanced [Eu/Fe] in low metallicity dwarf galaxies, which does not appear supported by observations (Geisler et al., 2005; Lanfranchi et al., 2008), and it appears in contrast to the apparent flatness of $\hbox{[Eu/Fe]}_{\rm pl}$ against [Fe/H] in the SAGA database at low metallicities (see Fig. 1). By contrast, we will show that the 2-phase ISM reproduces observed tracks with plausible ISM physics and parameter choices.

One could mimic the results of our 2-phase model by invoking differential outflow efficiencies for ccSN and NM products at early times that go away at late times, since this effectively changes the relative $\alpha$-element and $r$-process yields between sub-solar and solar metallicities. However, this scenario again requires contrived parameter tuning relative to the natural results of the 2-phase model.

6 Models and comparisons in abundance space

After these general discussions we perform some numerical experiments using the Schönrich & McMillan (2017) model and overplotting those on the collected data in the SAGA database (Suda et al., 2008, 2017). The bottom panel of Fig. 3 shows the evolution of [Mg/Fe] vs. [Fe/H] for three different radii ($4,6,10\,{\rm kpc}$) in the chemodynamical model. For these first plots, we have removed the hot gas phase, i.e. all stellar yields are either lost to the IGM or directly enter the star-forming gas. The evolution follows the usual picture: all trajectories run from left to right at $\hbox{[Mg/Fe]}_{\rm pl}$, forming the thick disc ridge. Once the SNeIa set in, the trajectories tend towards solar [Mg/Fe], while their different positions in [Fe/H] reflect the Galactic radial metallicity gradient. The top panel shows the modelled r-process abundance from neutron star mergers using different timescales vs. the [Eu/Fe] observations compiled in the SAGA database. The qualitative picture is straight-forwardly explained. The initial rise in [r/Fe] derives from the relative time-delay of NM vs. ccSNe. Trajectories can reach a high [r/Fe] plateau when both i) $\tau_{\rm NM}\ll\tau_{\rm SF}$, where $\tau_{\rm SF}=\left(d\ln(\rm SFR)/dt\right)^{-1}$, i.e. the star formation rate must not rise too quickly, ii) $\tau_{\rm NM}\ll\langle\tau_{\rm pop}\rangle$, where $\langle\tau_{\rm pop}\rangle$ is the average population age, and iii) $\tau_{\rm NM}\ll\tau_{\rm SNIa}$, i.e. the SNeIa do not significantly contribute on the NM timescale. Thus, only the shorter timescales with $\tau_{\rm NM}<150\,{\rm Myr}$ can form a high [r/Fe] plateau. In accordance with our approximations in Section 5, no matter how short we choose $\tau_{\rm NM}$, $\hbox{[Eu/Fe]}_{\rm pl}$ just approaches the value of $\hbox{[Mg/Fe]}_{\rm pl}$, but does not reach or exceed it. We will in later figures fix $\tau_{\rm NM}=150\,{\rm Myr}$, but this choice has little impact on our results as long as we keep $\tau_{\rm NM}\ll\tau_{\rm SNIa}$.

The tension between 1-phase models and data becomes even more evident in Fig. 4, which shows the behaviour of [Eu/Mg] (top left) and [Eu/Si] (top right) vs. [Fe/H]. For assessment of the calibration we show the respective [$\alpha$/Fe] vs. [Fe/H] trends in the bottom row. Here we use the merged catalogue from Bensby et al. (2014), Battistini & Bensby (2015) and Battistini & Bensby (2016), so that we have only data from one single, consistently analysed disc sample. While this does not remove all systematic problems, it reduces possible trends arising from heterogeneous data sets. In the lower panels $\hbox{[Si/Fe]}_{\rm pl}$ shows the slightly lower value expected for an element that has a slight production also from SNeIa. As argued in Section 5, the model [Eu/Mg] values never rise significantly above $\hbox{[Eu/Mg]}>0$, nor above $\hbox{[Eu/Si]}>0.05$, where the slight violation is permitted by the small SNIa contribution to Si yields. This is in stark contradiction to the data, which feature $\hbox{[Eu/Si]}\gg 0.05\,{\rm dex}$ and $\hbox{[Eu/Mg]}\gg 0\,{\rm dex}$, despite the fact that non-LTE analysis (Bergemann et al., 2017) shows that in the data presented here, $\hbox{[Mg/Fe]}_{\rm pl}$ and consequently the observed [Eu/Mg] values are biased by $\sim-0.1\,{\rm dex}$.

With Fig. 5 we approach the theme of this paper. The top and bottom panels show again the evolution of [Eu/Fe] and [Mg/Fe]. In comparison to the previous results we now introduce (model shown with red points) the hot phase of the ISM, and assume that most yields are fed into the hot ISM from which they have to cool down back into the cold gas phase. Compared to the 1-phase model, all trajectories shift to the left, since lock-up in the hot phase significantly lowers the metallicity in the cold star-forming gas at any given time $t$ that is comparable to or shorter than the cooling time $\tau_{\rm cool}$. In contrast to the [Eu/Fe]-[Fe/H] plane, the trajectories in the [Mg/Fe]-[Fe/H] plane are only weakly affected, since the SNIa timescale, which delays the additional Fe yields, is larger than $\tau_{\rm cool}$. This again confirms that our choice of a low value of ${\rm f_{\rm c,SNIa}}$ has no practical consequence to this paper. For ${\rm Eu}$ the situation is different: the neutron star merger timescale fulfills $\tau_{{\rm NM}}\ll\tau_{\rm cool}$, so here the lock-up becomes important. This shifts the initial (roughly linear, see Appendix A1, eq. 13) rise of the model trajectories in the [Eu/Fe]-[Fe/H] plane to the left, deeper into the regime dominated by stochastic chemical evolution. Further, along the thick disc plateau, there is also a very mild increase in [Eu/Fe]. It derives from the change in star formation rates, the different effect of NM delay time with the hot gas phase, and a small upward nudge at late times due to the lower ${\rm f_{\rm c,SNIa}}$. However, this effect is still far too small to enable this prescription to break the hard limits discussed in Section 5.

The possible solution is shown in Fig. 6. While keeping all other parameters constant in the 2-phase model presented in Fig. 5, we now allow a larger fraction of neutron star merger yields to enter the cold ISM than the yields from ccSNe, i.e.: $\rm f_{\rm c,NM}>\rm f_{\rm c,ccSN}$. Just as predicted in the simple analytical picture of equation (7), we now obtain plateau values fully in agreement with the [Eu/Fe] observations.

We stress that the goal is not to match the observations at low metallicities: our models predict very few thick disc stars formed in the linearly rising section of the model trajectories at $\hbox{[Fe/H]}\lesssim-1.5\,{\rm dex}$. The stellar populations with $\hbox{[Fe/H]}\lesssim-1.5\,{\rm dex}$ are dominated by accreted objects, i.e. stars from dwarf galaxies. At first order, simply imagine the presented model to be shifted to the left, as those smaller galaxies have lower star formation efficiencies and larger loss rates to the IGM. Many of them will also cease to form stars before evolving to low [Mg/Fe] and [Eu/Fe] values. We thus can expect that the high [Eu/Fe] plateau is extended to $\hbox{[Fe/H]}\lesssim-1.5\,{\rm dex}$ by a superposition of analogous dwarf galaxy models. Fitting this region would also not add new information to our case: the low-[Fe/H] regime is not only composed of a large number of different systems, but the different flavours of stochastical chemical evolution, and uncertainties e.g. about the incidence of magneto-rotational ccSNe, add a large space of additional parameters that do not pertain to our question, which is answered on the more tightly constrained thick disc stars.

In Fig. 7 we present an extension of Fig. 6. We now present for [Eu/Fe] (top row) and [Mg/Fe] (bottom row) a range of values for $\rm f_{\rm c,NM}$ for two possible values of $\rm f_{\rm c,ccSN}=0.25$ (left-hand column) and $\rm f_{\rm c,ccSN}=0.5$ (right-hand column). The right-hand side should be considered only as an instructive comparison, as we do not favour such a large cold ISM channel from ccSNe. We can see by the comparisons that eq. (7) correctly predicts the behaviour. Multiplying $\rm f_{\rm c,ccSN}$ by a factor $2$ shifts $\hbox{[Eu/Fe]}_{\rm pl}$ downwards by $0.3\,{\rm dex}$. Similarly at fixed $\rm f_{\rm c,NM}$ the plateau would shift upwards by $0.3\,{\rm dex}$ compared to the left-hand column of Fig. 7 if we chose $\rm f_{\rm c,ccSN}=0.125$. The same shift applies for every factor $2$ in $\rm f_{\rm c,NM}$. A slight caution needs to be taken as we can see from the high-metallicity/low-[Eu/Fe] end of the trajectories. For large values $\rm f_{\rm c,NM}$, less r-process elements stay locked in the hot phase, which would require a mild correction to our adopted normalisation to bring the trajectories back to $\hbox{[Eu/Fe]}=0$. However, this is a small fraction of the change in the plateau, and readers can easily convince themselves that the difference between the trajectory-endpoints and the plateau rises strongly with $\rm f_{\rm c,NM}$ irrespective of the correction.

In analogy to Fig. 4, Fig. 8 analyses model trajectories in the [Eu/$\alpha$]-[Fe/H] plane when we vary $\rm f_{\rm c,NM}$, affirming the points made above. When we allow for a larger fraction of neutron star merger yields to be contributed directly to the cold gas phase, the model values of [Eu/Mg] and [Eu/Si] finally push significantly into positive territory. ${\rm Si}$ requires a somewhat larger effect than ${\rm Mg}$, which we again ascribe to the problem of [Mg/H] values being slightly overestimated in a 1D-LTE spectroscopic analysis. The figure also underlines again that different yields to the hot/cold gas phases readily account for the observed trends.

Fig. 9 tests different cooling times, $\tau_{\rm cool}$, for the hot ISM, for an otherwise standard model setting with $\rm f_{\rm c,ccSN}=0.25$ and $\rm f_{\rm c,NM}=0.61$, showing the comparison with both the SAGA database in [Eu/Fe] (top) and the Battistini & Bensby sample in [Eu/Si] (bottom). We consider timescales below $\sim 500\,{\rm Myr}$ and significantly above $2\,{\rm Gyr}$ unlikely, since they would either approach values shorter than the orbital timescale or predict too much mass in the warm/hot ISM, but it is instructive to compare the model values. For very short cooling timescales the model of course reverts back to the behaviour of a 1-phase model as shown in Fig. 3. For very long timescales the model approaches the upper limit given in equation 7, and the changes nearly saturate. We are not re-calibrating the model’s loss rates and r-process yields between the different parameters – one can see that the larger lock-up in models with large $\tau_{\rm cool}$ results in slightly lower end-point [Fe/H] and larger [Eu/Fe]. If we re-gauged the parameters, models with $\tau_{\rm cool}>1\,{\rm Gyr}$ would look nearly identical, though requiring slightly smaller yield losses and smaller r-process yield for larger $\tau_{\rm cool}$. When comparing to the SAGA database (top panel) we plot lower [Fe/H] values. It is clear that $\tau_{\rm cool}$ is less important for the very low [Fe/H] regime, where $\tau_{\rm cool}$ is large compared to the age of the stellar population. Most importantly, our assumption on $\tau_{\rm cool}$ impacts the estimated ratio of $\rm f_{\rm c,NM}/{\rm f_{\rm c,SNIa}}$, which is insensitive against changes in the region $\tau_{\rm cool}>1\,{\rm Gyr}$, but would have to be increased for shorter $\tau_{\rm cool}$, e.g. to about $4$ if one assumed $\tau_{\rm cool}=500\,{\rm Myr}$.

Readers could now argue that we do not know what the relative contributions to the respective gas phases of each yield class are, and they would be right. The point is that we can a priori expect differences, i.e. $\rm f_{\rm c,NM}\neq\rm f_{\rm c,ccSN}$, but we have no way to quantify them at present, apart from looking at the abundance plane. Thus, to complement this discussion, Fig. 10 compares the chemical evolution in a model with both gas phases and treating the neutron-star mergers like ccSNe, i.e. $\rm f_{\rm c,NM}=\rm f_{\rm c,ccSN}=0.25$ to a model with both gas phases in which neutron-star merger yields go entirely into the hot phase. Obviously the latter model strongly exacerbates the problem that neutron-star mergers contribute their yields later than the ccSNe, and the low [Eu/Fe] trajectories strongly disfavour this scenario, even if we could find other ways (e.g. metallicity dependent yields) to help raise the [Eu/Fe] ratios at early times.

7 Conclusions

We should take three main points from this analysis:

• •

Modelling the different phases of the interstellar medium (hot vs. cold star forming gas) is vital to understand chemical evolution on timescales smaller than $\sim 1\,{\rm Gyr}$.

• •

In contrast to 1-phase chemical evolution models, we find that neutron star mergers as source of r-process elements with reasonable delay time distributions (delay times of order $100\,{\rm Myr}$) can explain the observed abundance patterns, provided that the fraction of NM yields delivered directly to the cold star-forming phase of the ISM is higher than that of ccSN yields.

• •

The only other significant sources of r-process elements are (possibly) ccSNe, but we can only explain the $\hbox{[Eu/Si]}>0$ and $\hbox{[Eu/Mg]}>0$ values in the thick disc if there is a source of r-process elements that differs from a constant ccSN contribution. This implies a significant source of r-process elements besides ccSNe - most naturally neutron star mergers.

The main problem is not as much explaining r-process abundances at the very metal-poor end, where stochastic chemical evolution and the superposition of stellar populations from different accreted dwarf galaxies complicate the picture (and provide freedom in parameter choice). Instead, the challenge is the relatively high-metallicity edge of the high [Eu/Fe] sequence, near the “knee” of the [Mg/Fe] vs. [Fe/H] distribution.

A cautionary point that would allow us to lower the NM contribution to the r-process budget is the possibility of metallicity dependent yields from ccSNe/collapsars. We did not have space here for an in-depth discussion. However, these models require a lot more fine-tuning and face three observational challenges that are not yet convincingly resolved: a significant increase of r-process yields at low metallicities would have to i) still conform with normal [Eu/Fe] values observed in dwarf galaxies, and to ii) be fine-tuned in order not to incline the thick disc plateau. Further, a collapsar model with a metallicity dependence on the metal-rich end, as suggested in Siegel et al. (2019) would have to still produce the radial [Eu/Fe] profile measured in the Milky Way.

We have shown with simple analytical considerations in Section 5, a more comprehensive analytical model in Appendix A1, and with the chemodynamical models of Schönrich & McMillan (2017) in Section 6, that \hbox{[Eu/\alpha]}>0 or \hbox{[Eu/Fe]}>\hbox{[\alpha/Fe]}, cannot be achieved in a 1-phase chemical evolution model in general. However, they are readily reproduced with a simple chemical evolution model as soon as we account for a 2-phase ISM and allow yields from neutron star mergers to enter the cold, star forming gas at a mildly higher fraction than the ccSN yields, $\rm f_{\rm c,NM}>\rm f_{\rm c,ccSN}$. The factor $\rm f_{\rm c,NM}/\rm f_{\rm c,ccSN}$ required here can be directly estimated from the [Eu/$\alpha$] value of the thick disc knee to be about $0.2\,{\rm dex}$ or $\rm f_{\rm c,NM}/\rm f_{\rm c,ccSN}\approx 2$.

We find that the this factor is quite insensitive to cooling timescales $\tau_{\rm cool}\gtrsim 1\,{\rm Gyr}$ and would need to be increased for shorter $\tau_{\rm cool}$. There is hope to derive better constraints for $\tau_{\rm cool}$ from the chemical evolution of elements produced on different timescales and by different sources, e.g. s-process elements.

The slightly different behaviour of [Si/Fe] and [Mg/Fe] in observations serves as a reminder that one should not blindly trust observed abundance trends. In particular, trends vs. [Fe/H] are vulnerable to increasing deviations from the frequently assumed LTE, when with lower metallicity the importance of collisions in the photosphere’s plasma decreases. Similarly, we calibrated this model on [Mg/Fe] and thus have a mismatch by $\sim~{}0.07\,{\rm dex}$ with [Si/Fe]. This has no impact on our qualitative results, but it puts systematic uncertainty on the exact value needed for $\rm f_{\rm c,NM}/\rm f_{\rm c,ccSN}$. Further, the qualitatively consistent behaviour of \hbox{[Eu/\alpha]}>0 and \hbox{[Eu/Fe]}>\hbox{[\alpha/Fe]} for several different $\alpha$ elements makes it unlikely that the problem is a mere result of our still incomplete understanding of stellar spectra. Non-LTE analyses of ${\rm Eu}$ are very rare, though from Zhao et al. (2016) it seems that ${\rm Eu}$ and [Eu/Fe] measurements are quite robust and if at all increase vs. LTE analysis. More systematic studies are quite sorely needed.

The re-distribution of r-process elements between different gas phases and the resulting changes to abundance trends might at first hand seem like an annoying complication, but it may be turned into an important/useful diagnostic if simulations have an appropriate recipe for hot/cold gas dynamics and feedback. In this light, the strong discrepancy in papers like van de Voort et al. (2015) between the modelled and observed trends of [Eu/Fe] vs. [Fe/H] in the abundance plane is a chance to improve on the ISM modelling in hydrodynamical simulations, and our paper may provide an analytic framework to understand the differences between these and other hydro simulations that matched the data better (e.g. Shen et al., 2015).

We point out that the 2-phase ISM assumed here is a simplified picture of a real galaxy, which will have a whole range of delay times that will vary with how far from the star-forming regions material is expelled and how long it takes for yields to re-enter the star-forming phase. So some of the difference in $\rm f_{\rm c,NM}$ vs. $\rm f_{\rm c,ccSN}$ may be due to ccSN yields being expelled further from the disc.

More generally, we conclude that observed r-process abundance trends provide empirical evidence for both an important NM contribution to r-process enrichment and differential distribution of yields to the cold and hot phases of the ISM.

Acknowledgements

We thank our referee F. Thielemann for enlightening comments. It is a pleasure to thank Eliot Quataert for pointing out the problem of neutron star chemical evolution. We would like to thank M. Bergemann, J. Binney, J. Johnson, J. Magorrian, B. Metzger, N. Prantzos, D. Siegel, and S. Smartt for helpful comments. RS is supported by a Royal Society University Research Fellowship. This work was performed using the Cambridge Service for Data Driven Discovery (CSD3), part of which is operated by the University of Cambridge Research Computing on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). The DiRAC component of CSD3 was funded by BEIS capital funding via STFC capital grants ST/P002307/1 and ST/R002452/1 and STFC operations grant ST/R00689X/1. DiRAC is part of the National e-Infrastructure.

Appendix

A1 Analytical model

In this section we develop a more sophisticated analytical model for Eu vs. Mg enhancement. As an approximate guide to the behaviour of [Eu/Mg] abundance ratios with a 2-phase ISM, we can adapt equations (34) and (37) of Weinberg et al. (2017, hereafter WAF), which describe the abundance evolution for an element with instantaneous enrichment and an element whose production follows an exponential delay time distribution with a minimum time delay $t_{\rm min}$.

These equations assume a constant SFR, a constant star formation efficiency timescale $\tau_{*}\equiv M_{g}/\dot{M}_{*}$, and outflow with a constant mass-loading factor $\eta\equiv\dot{M}_{\rm out}/\dot{M}_{*}$. With some notational changes the enrichment in the instantaneous and delayed cases can be written

[TABLE]

and

[TABLE]

where $r\approx 0.4$ is the recycling fraction, $\Delta t=t-t_{\rm min}$, and $\tau_{\rm d}=\tau_{*}/(1+\eta-r)$ is the timescale on which the gas supply would be depleted (in the absence of accretion) by the combination of star formation and outflow. For present purposes, we use equation (9) multiplied by $\rm f_{\rm c,ccSN}$ to represent the ccSN Mg enrichment direct to the cold phase, and we use equation (10) multiplied by $(1-\rm f_{\rm c,ccSN})$ to represent ccSN Mg enrichment that passes through the hot phase, with $\tau_{\rm del}=\tau_{\rm cool}$. We use equation (10) multiplied by $\rm f_{\rm c,NM}$ to represent NM Eu enrichment direct to the cold phase, with $\tau_{\rm del}=\tau_{{\rm NM}}$ and a minimum delay time $t_{\rm min,NM}$, and we use equation (10) multiplied by $(1-\rm f_{\rm c,NM})$ and with timescale $\tau_{\rm del}=\tau_{{\rm NM}}+\tau_{\rm cool}$ to represent NM Eu enrichment that passes through the hot phase. In similar fashion, we can compute the ccSN and SNIa contributions to iron enrichment.

Although the assumptions and parameters of this approximate analytic model differ from those of our numerical model, we find that it reproduces the qualitative behaviour of the numerical [Eu/Mg] and [Eu/Fe] results for different choices of $\tau_{{\rm NM}}$ and $\rm f_{\rm c,NM}$.

For typical yield datasets (Andrews et al., 2017, or the yields used here), an outflow mass-loading $\eta\approx 2.5$ is needed to obtain approximately solar abundances at late times. This factor would be moderately modified in the presence of radial flows. Figure 11a shows the [Mg/Fe] evolution (dotted curve) and [Eu/Fe] evolution (solid curves) for a model with $\eta=2.5$, $\tau_{*}=1\,{\rm Gyr}$, $\rm f_{\rm c,NM}=\rm f_{\rm c,ccSN}=0.25$, and a range of choices for $\tau_{{\rm NM}}$. We have adjusted the ccSN iron yield to give a plateau at $\hbox{[Mg/Fe]}_{\rm pl}\approx 0.45$ similar to the numerically solved chemical evolution in Section 6, and we have used the same SNIa parameters (${\rm f_{\rm c,SNIa}}=0.01$, $\tau_{\rm SNIa}=1.5\,{\rm Gyr}$, minimum time delay $0.15\,{\rm Gyr}$), cooling time $\tau_{c}=1\,{\rm Gyr}$, and minimum NM delay ($0.015\,{\rm Gyr}$) as the numerical models. Comparison to Figure 3 shows similar trends with $\tau_{{\rm NM}}$, similar peak values of [Eu/Fe], and similar evolution for $\hbox{[Fe/H]}\ga-1.5$. At times $t\gg\tau_{{\rm NM}}$, ccSN and NM enrichment proceed at equal rates, so all [Eu/Fe] curves asymptote to match the [Mg/Fe] curve, the shape of which is driven by SNIa iron enrichment.

At early times and low [Fe/H], the analytic predictions of [Eu/Fe] in Figure 11a are significantly lower than those in Figure 3, primarily owing to the different treatment of outflows. The model presented in the main text assumes a fixed loss of stellar yields ($60\%$) vs. the mass loss ($\eta$) formalism applied in this section, which drives an outflow at the ISM metallicity. In the analytic model, less metals are lost at early times, when the ISM metallicity is low. The ccSN enrichment of the $\eta=2.5$ model therefore proceeds faster at early times, moving the initial rise in [Eu/Mg] and [Eu/Fe] to the right. We replicate the difference between the models by multiplying all yields by a factor of 0.4 and setting $\eta=0$ in Figure 11b, which as expected shifts the early rise in [Eu/Fe] back to the left, approaching Figure 3, though if we changed to a single-phase model ($\rm f_{\rm c,ccSN}=\rm f_{\rm c,NM}={\rm f_{\rm c,SNIa}}=1$) then the early rise would still be slower than in Figure 3. The remaining differences may be explained by the instantaneous return of locked up metals in the analytic approach of this Section vs. the correctly delayed return allowed by the numerical solution. Radial gas flows in the model of Section 6 depress the local metallicity, so that its net mass loss rate is smaller. Most important for our purposes, the maximum values of [Eu/Fe] are nearly identical for both outflow implementations.

Figures 11c and 11d show varying values of $\rm f_{\rm c,NM}$ with $\tau_{{\rm NM}}$ fixed to $0.15\,{\rm Gyr}$, for the same $\eta$ and yield combinations in panels (a) and (b), respectively. Here we adopt $\rm f_{\rm c,ccSN}=0.25$, as in the left column of Figure 7. The behaviour of the analytic models resembles that of the numerical models, though the values of [Eu/Fe] are again lower at low [Fe/H], and the maximum values of [Eu/Fe] are slightly lower (by $\sim 0.05$ dex). The final panels adopt $\rm f_{\rm c,NM}=0.50$ and $\tau_{{\rm NM}}=0.15\,{\rm Gyr}$ and vary the outflow mass-loading $\eta$ or the hot-phase cooling time $\tau_{c}$. Changing $\eta$ alters the final [Fe/H] but not the maximum of [Eu/Fe]. Shortening $\tau_{c}$ to $0.1\,{\rm Gyr}$ allows hot-phase enrichment to occur much more quickly, so the model has already evolved to $\hbox{[Fe/H]}\approx-0.5$ before [Eu/Fe] reaches its maximum. However, models with $\tau_{c}=0.5$ or $2.0\,{\rm Gyr}$ are only slightly different from our standard case of $1.0\,{\rm Gyr}$. The results in Fig. 11f qualitatively resemble those in Fig. 9.

We can use the analytic model to gain some insights into behaviour in different regimes. For a star formation efficiency timescale $\tau_{*}=1\,{\rm Gyr}$ and $\eta=2.5$, the depletion timescale is $\tau_{d}\approx 0.3\,{\rm Gyr}$. For $\tau_{c}\approx 1\,{\rm Gyr}$, we expect $\tau_{c}>\tau_{d}>\tau_{{\rm NM}}$. If we make the approximation $\tau_{\rm cool}\gg\tau_{\rm d}\gg\tau_{{\rm NM}}$, implying $\tau_{\rm cool}/(\tau_{\rm cool}-\tau_{\rm d})\approx 1$ and $\tau_{{\rm NM}}/(\tau_{{\rm NM}}-\tau_{\rm d})\approx 0$, we get

[TABLE]

where

[TABLE]

When $t\ll\tau_{\rm cool}$ so that the first terms dominate, the abundance ratio is equal to the yield ratio multiplied by $\rm f_{\rm c,NM}/\rm f_{\rm c,ccSN}$. When $t\gg\tau_{\rm cool}$, the abundance ratio is simply equal to the yield ratio, which should be similar to the solar Eu/Mg ratio if the model is to reproduce observations.

The approximation that leads to equation (11) breaks down when $t\ll\tau_{{\rm NM}}$. In this regime, one can instead ignore the hot phase return and use a Taylor expansion of equations (9) and (10) to find

[TABLE]

(for $t>t_{\rm min,NM}$). Thus, the Eu/Mg ratio grows approximately linearly at early times, leveling off as $t$ approaches $\tau_{{\rm NM}}$. By this time, the Mg abundance has reached a fraction $\rm f_{\rm c,ccSN}(1-e^{-\tau_{{\rm NM}}/\tau_{\rm d}})\approx\rm f_{\rm c,ccSN}(\tau_{{\rm NM}}/\tau_{\rm d})$ of its final equilibrium value, assuming that hot phase return is still negligible. From these analytic arguments, we conclude that $\rm f_{\rm c,NM}>\rm f_{\rm c,ccSN}$ is required to produce [Eu/Mg]$>0$ at any time and that the condition for reaching [Eu/Mg]$>0$ by the time [Mg/H]$\approx-1$ is approximately $\rm f_{\rm c,ccSN}(\tau_{{\rm NM}}/\tau_{\rm d})<0.1$. Both conditions are physically plausible. Reaching [Eu/Mg]$>0$ while [Mg/Fe] is significantly super-solar requires $\tau_{{\rm NM}}<\tau_{\rm Ia}$, which is again physically plausible.

A2 Further model tests

In Fig. 12 we present models with a 1-phase ISM (left-hand column) and a 2-phase ISM (right-hand side) for three different SNIa DTDs, as discussed in Weinberg et al. (2017). These double-exponential DTDs are designed to more closely match the frequently claimed $t^{-1}$ time-dependence, while they preserve the straight-forward normalisation and integrability of exponential functions.

While there might be some targeted applications where this matters, the short message for our work is: it does not significantly impact our results.

Finally, Fig. 13 shows different characteristic delay timescales $\tau_{\rm NM}$ for models with our standard cooling timescale of $\tau_{\rm cool}=1\,{\rm Gyr}$ and more yields from neutron star mergers directly entering the cold gas phase ($\rm f_{\rm c,NM}=0.61$ vs. $\rm f_{\rm c,ccSN}=0.25$. The behaviour is as expected, with the models with shorter $\tau_{\rm NM}$ showing faster ${\rm Eu}$ enrichment and thus reaching a larger peak [Eu/Fe] value, before SNeIa and cooling of less r-process rich material from the hot phase set in.