Metallicity-Suppressed Collapsars Cannot be the Dominant r-Process Source in the Milky Way

We develop a high-performance analytical model of Galactic Chemical Evolution, which accounts for delay time distributions and lock-up of stellar yields in a thermal-phased ISM. The model is capable of searching, for the first time, through the high-dimensional parameter space associated with the r-process enrichment of the Milky Way by its possible sources: Neutron Star Mergers and Collapsar events. Their differing formation mechanisms give these two processes different time dependencies, a property which has frequently been used to argue in favour of collapsars as the dominant r-process source. However, we show that even with large degrees of freedom in the allowed thermal, structural, and chemical properties of the galaxy, large regions of parameter space are in strong tension with the data. In particular, whilst we are able to find models in which neutron star mergers produce the majority of r-process material, the data rule out all models with dominant collapsar yields. With no other identified source, we conclude that Neutron Star Mergers must be the dominant contributors to the modern Milky Way r-process budget.


INTRODUCTION
Since the landmark paper of Burbidge et al. (1957), it has been widely accepted that, in order to explain the abundance distribution of chemical elements observed in the universe, we require a number of distinct nucleosynthesis channels, operating in unison.
Primordial nucleosynthesis sourced the lightest elements in the universe (H and He, Alpher et al. 1948), while the heavier elements are created in processes such as shell burning in the stellar interior (the elements, including O, Si and Mg, Hoyle 1954), explosive nucleosynthesis during supernovae events ( elements, plus the iron peak elements: Fe, Ni, Co, Arnett & Clayton 1970) or cosmic ray spallation (Li, Be, B, Reeves et al. 1970).
The majority of the elements heavier than iron, however, are sourced from a variety of neutron-capture processes: the slow (s), intermediate (i) and rapid (r) neutron capture processes. Whilst the origin of the s-process is well understood (Clayton et al. 1961), and the i-process thought to contribute significantly to only a handful of isotopes (Côté et al. 2018), the debate about the astrophysical sites that can lead to sufficient r-process synthesis has remained an open and enduring question for many years, with several possible sites for the r-process being identified: Woosley et al. (1994) argue that neutrino heating in Core Collapse Supernovae (CCSN) creates favourable conditions such that normal CCSN can provide a site for the r-process. This model, however, is plagued by overproduction of certain elements, and the high entropy conditions required have been questioned in more recent studies (i.e. Fischer et al. 2012). In addition, the existence of ultra metal poor, but highly r-process enriched stars (such as that found in Sneden et al. 1996) indicates that the source of r-process nucleosynthesis must be a rare, high yield event. CCSN are therefore not considered a viable source of r-process enrichment.
The disruption of a neutron star by tidal interactions during a merger with a black hole (Lattimer & Schramm 1974), or by a binary collision between two neutron stars (Symbalisty & Schramm 1982;Freiburghaus et al. 1999;Rosswog et al. 1999) are also candidates for the r-process. The detection of the combined gravitational wave GW170817 and GRB event GRB170817A, confirmed to arise from a NS-NS merger event (Abbott et al. 2017), and the subsequent detection of r-process material in the ejecta (Chornock et al. 2017;Rosswog et al. 2018) provided direct observational evidence that Neutron Star Mergers (NSM) produce r-process material.
Though the existence of NSM as an active r-process pathway is rarely called into question, it is seen as concerning that time-delayed nature of NSM formation would naïvely predict entirely different enrichment pathways in [Eu/Fe]-[Fe/H] space than is observed, leading to either the conclusion that NSM cannot be dominant r-process sources (Argast et al. 2004;, or the invocation of neutron star properties incompatible with their understood behaviour (Matteucci et al. 2014).
In an alternative approach, Fujimoto et al. (2006) combined an MHD jet method with the collapsar models of Woosley (1993), and demonstrated that this can produce a significant amount of r-process enrichment. Collapsars occur when the core of a collapsing star is rotating sufficiently fast to delay radial infall, resulting in MHD jets driven by accretion onto a compact engine, and are thought to be the source of Long Gamma Ray Bursts. Although LGRB are welldocumented events and often tied to unusual forms of CCSN due to their formation in regions of rapid star formation (Bloom et al. 2002) and several closely tied observations of supernovae associated with LGRBs (Kulkarni et al. 1998;Mazzali et al. 2003;Sollerman et al. 2006), there is no direct evidence linking their formation with rprocess synthesis. However, the model is widely favoured, since the high progenitor masses imply a short lifetime and thus allow for very early r-process enrichment.
Other sources for r-process material have also been studied. For example, neutrino-driven winds (Wanajo et al. 2001) and electron-capture supernovae (Wanajo et al. 2011), however Haynes & Kobayashi (2018) showed that these did not produce sufficient quantities of r-process material. Whilst in the case of the intermediary (i) process-producing White Dwarf binaries (Cowan & Rose 1977;Denissenkov et al. 2017), it was predicted by Côté et al. (2018) that only specific isotopes are produced in significant quantities, with 45 per cent of solar Mo predicted to be i-process in origin, but less than 10 −2 of solar Eu. For the sake of clarity and simplicity, we will therefore neglect these sources.
Confusing matters further, there is also evidence of an incomplete (or weak) r-process (Honda et al. 2006), in which the lighter r-process elements are synthesised, but not the heavier second and third peak elements. Magneto-rotational supernova discussed in Nishimura et al. (2017) (also referred to as 'hypernovae', though this phenomenological term can refer to collapsars) or the recently proposed Quark-Deconfinement Supernovae of Fischer et al. (2020) are thought to be good candidates. For our purposes, we use the term 'r-process synthesis' to refer to the complete r-process, in which all r-process material up to the third peak is synthesised.
It might feel natural to assume that multiple pathways actively contribute to the r-process enrichment, with collapsars providing the early time yield, and then NSM coming in later. However, this bears two problems: i) the enrichment profiles for [Eu/Fe] are poorly replicated in simple GCE models whenever NSM are significant contributors, and ii) Sneden et al. (1996) demonstrated a common r-process fingerprint: a remarkable consistency of relative r-process abundances, with Sneden et al. (2008) (henceforth Sn08) extending this relationship. The relative abundances of r-process material in the metal-poor but highly r-process enriched star CS 22892-052 match those of the solar system, despite the fact that the high enrichment indicates very early enrichment from the unmixed ejecta of a single r-process event. The common 'fingerprint' with the solar system implies that the material is produced in the same ratio throughout galactic history that dominates today's abundances, a tension if one would like to assume that the dominant 1 channel switched from collapsars early on to NSMs today.
In this paper, we will put quantitative limits on the relative contributions from both neutron star mergers and collapsars by comparing chemical evolution models with observed stellar abundances. This will show that under reasonable assumptions the relative contribution from collapsars is highly limited.
Section 2 will introduce some of the key aspects of collapsar formation as it pertains to chemical evolution, with section 3 detailing the observational data that we will attempt to replicate in our models. Section 4 introduces the analytical SACEM model, and briefly discusses the full evolutionary simulation RaMiCES, whilst Section 5 outlines our attempt to eliminate regions parameter space, and §6 and §7 discuss our findings regarding the excluded regions of parameter space, and the required properties of our models.

MODELLING COLLAPSARS & THEIR YIELDS
Collapsars are a corollary to the existence of the 'failed supernovae' proposed in Bodenheimer & Woosley (1983) and since observed by 1 Throughout this paper we use the nomenclature that a dominant source is one which can be assumed to be the sole source, neglecting all others. A non-dominant source which produces more than 50 per cent of production is a majority source. Adams et al. (2017). 'Failed Supernovae' are the fate of stars which are so massive that the usual supernova mechanism is insufficient to prevent runaway radial collapse. However, progenitors with large angular momentum cannot collapse spherically. Instead, they collapse into a compact accretion disc around the growing black hole: a potential site for r-process nucleosynthesis and LGRBs.

Dependence on Initial Metallicity
The nature of collapsars necessitate a mass and core angular momentum cutoff for their formation: the star must be massive enough to defy normal supernova mechanisms, but have enough internal angular momentum to stave off radial collapse. There are also strong indications of a metallicity dependence on collapsar rates, both from models of how metallicity impacts the stellar interior (through variations in the opacity and angular momentum transport efficiencies), and from observational LGRB counts such as Perley et al. (2016). Theoretical models of stellar evolution quote a strict cutoff in metallicity around c = 0.1 ⊙ (MacFadyen & Woosley 1999;Yoon & Langer 2005;Woosley & Heger 2006), up to c = 0.3 ⊙ (Yoon et al. 2006). In these models, stars with > c undergo core-braking to below the critical threshold, and hence cannot be progenitors of collapsar events.
In contrast, whilst observational constraints of GRB events points to find suppressed collapsar activity at higher , they instead find a much larger cutoff value, as high as = ⊙ (Wolf & Podsiadlowski 2007), though they do find suppression beginning beforehand. Perley et al. (2016) propose that this discrepancy arises because there are both single-progenitor and double-progenitor pathways for GRB events, with the single progenitor pathway (isolated collapsars) dying away fastest.
Some proposed double-progenitor LGRB models (i.e. Podsiadlowski et al. 2010), however, source their energy from explosive detonation inside a common envelope rather than accretion disc/central engine interactions, making a significant contribution to r-process nucleosynthesis questionable. In addition, some observed LGRB events have no accompanying supernova detected (Fynbo et al. 2006), further indicating a diversity of origins for LGRBs beyond r-process-producing collapsars.
In summary, current observational evidence indicates that rprocess producing collapsar events form a non-constant subset of all LGRB events. This is concerning as many studies (i.e. Siegel et al. 2019, henceforth Si19) assume that the r-process nucleosynthesis tracks the GRB rate perfectly: such models thus might source rprocess elements from collapsars long after they have stopped contributing in reality.

r-Process Yields
In principle, one may therefore continue to search for conditions which would lead to collapsar formation, and couple these with rotational, metallicity and stellar-mass distribution functions and yield tables to simulate the 'true' collapsar contribution. However, this approach is unfeasible for a number of reasons: • The models of i.e. Heger et al. (2000) and Limongi & Chieffi (2018) show a complex/non-linear relationship between initial rotation and final core angular momentum, barring simple predictions for the likelihood of a star going collapsar from initial conditions.
• The physical processes occuring in the interior of massive stars are highly inaccessible. As a result, exact numerics of the yields from a given system must be considered poorly unconstrained.
Without rigorous observational or theoretical constraints, our models will therefore need to utilise an approximation for coll ( , , zams ). We choose the simplest approximation: a constant yield which is suppressed at a metallicity = c (Eq.C14). We justify this approach in §5.3 and conclude that it has neglgible impact on the strength of our conclusions.

PATTERNS IN THE R PROCESS ABUNDANCES
In this section we delineate the main features of the r-process abundances, both in the [Eu/Fe]-[Fe/H] plane, and as a function of time. The use of Eu as a proxy for the total r-process enrichment is justified in Appendix A.

The Abundance Planes
As stellar ages are hard to obtain accurately, we follow Tinsley (1979) by studying the chemical history of the Galaxy in the [X/Fe]-[Fe/H] plane, where: X is the number density of species X, and ⊙ denotes solar values. Figure 1 shows a sample of 965 stars drawn from the SAGA database 2 which possess both upper and lower bounds for all of the elements of interest, and have [Fe/H] > -2.5, since the low metallicity end of the distribution is dominated by stochastic processes, and by stars in the galactic halo which formed within dwarf galaxies accreted during the growth of the Milky Way. The low metallicity end therefore likely represents a superposition of the chemical histories of these dwarf galaxies (Ojima et al. 2018), rather than the in-situ history of the Milky Way. Due to the low mass of such accreted objects, their final impact on galactic chemistry is negligible, and this is therefore beyond the scope of this paper to discuss.
The chemical data in Fig. 1 are sourced from a variety of surveys, and encompasses studies of disc FGK dwarfs (Mishenina et al. 2013;Reddy et al. 2006), dwarfs in both the disc and the halo (Fulbright 2000), as well as studies of both giants and dwarfs in the galactic halo (Sakari et al. 2018;Allen et al. 2012), thereby providing us with a wide sample of the enrichment of the galaxy both in physical and temporal space. Following the prescriptions of Bergemann et al. (2017) and Zhao et al. (2016), we performed a minor NLTE correction for both Mg and Eu for stars with [Fe/H] < -2, though this affected our mean trends by less than 0.02 dex, and so did not meaningfully alter our results.
The  Fig. 1 may find the discussion in Appendix B helpful.

Time Dependence
We leverage the data from Bensby et al. (2014) and Battistini & Bensby (2016)  with stellar age. These papers provide age and abundance estimates for 714 G and F dwarfs in the Solar Neighbourhood. A subset of 339 stars has ages and abundances for [Eu/Fe] with both upper and lower limits. These datapoints are plotted in the top panel of Fig. 2.
The dataset also assigns membership probabilities of each star to the Galactic thick and thin discs based on kinematics only. However, due to the large overlap of disc components in kinematics, we instead prefer a chemical cut in the [Mg/Fe]-[Fe/H] plane, while we also tested that our conclusions remain unchanged irrespective of this strategy.
We used a standard Bayesian approach to estimate the best-fit gradient of a given sample of stars, marginalising over unknown errors (the dataset provided only [Fe/H] uncertainties, not [Eu/H]).  Battistini & Bensby (2016). Top: The compiled data showing the inferred stellar ages and [Eu/Fe] abundances. Bottom: a plot of the calculated abundanceage gradients of each temporal subsample as a function of the maximum age permitted into the set. Bins shaded grey have fewer than 10 stars in them and so are liable to large errors.
Using 2-Gyr long sampling periods, we build up a picture of how this gradient changes with the age of the stars, shown in the lower panel of Fig 2. We see that at early times the time derivative d[Eu/Fe]/d strongly positive -around +0.02 dex per Gyr for the combined sample and +0.01 dex per Gyr for the chemical thin disc. A positive gradient with stellar age is equal to a negative gradient with respect to forward time, so this shows that the [Eu/Fe] abundance was decreasing during this period.
Between 7Gyr and the present day, however, gradient has decreased such that the average change in [Eu/Fe] is consistent with zero across this time period, indicating that chemical equilibrium was reached approximately 7 Gyr ago.
We note that the data of Fig. 2 is derived under the assumption of LTE, and so both the age and abundance estimates may change under a full NLTE treatment. However, we note from Fig. 15 of Zhao et al. (2016) that significant NLTE corrections for the elements in question occur (if ever) for stars with [Fe/H] −1, which implies that stars younger than 10Gyr are safe against this bias. Our conclusions about the recent equilibrium of the galaxy are therefore robust against the LTE approximation.

Star Formation Rate
Observational studies show that the Milky Way has sustained a Star Formation Rate until recent times which is no less than an order of magnitude below its maximal value. Though the exact value of the MW star formation rate is contested, the general consensus is that it is in the region 1 − 2 ⊙ yr −1 (Mor et al. 2019;Chomiuk & Povich 2011;Robitaille & Whitney 2010;Aumer & Binney 2009).

Observational Summary
Usually models are fit to the data with the goal of optimizing to a best-fit parameter set, which is then prone to systematic biases. In this work, we walk a different path, where we instead try to falsify classes of models, based on their ability to reproduce a minimum set of constraints, which we draw from the observational evidence.
For each constraint, we present a broad conclusion which can be drawn from the data (in bold), followed by how this would be replicated within a model (in italics): (iv) Galactic chemistry is (almost) in equlibrium (v) Star formation has continued until late times

MODELS
In all our models we make the simplifying approximation that there are only 3 sources of nucleosynthesis for r-process elements: a small contribution from Core Collapse Supernovae 3 , Neutron Star Mergers, and Collapsars. Whilst other sources of r-process may exist and be important for explaining the abundance of individual stars, we expect the majority of such sources to be of subdominant importance in the case of the overall r-process trends in the galaxy. In the case of BH-NS mergers, we note that their formation mechanisms and timescales are so similar to NS-NS mergers that we can consider them part of the same process, though we note from Pannarale & Ohme (2014) that for BH 14 ⊙ , the vast majority of NS companions will be swallowed without tidal disruption, and hence without r-process ejecta.

Simple Argument
Before deriving our analytical chemical evolution models, it is instructive to first discuss the qualitative appearance of Collapsar In agreement with previous works (i.e. Si19), we parametrise rprocess yields from collapsars (see eq. 3) with a constant synthesis rate for metallicity < c . This assumption guarantees that the dominant collapsar pathway produces a plateau in both the [Eu/Fe] and [Eu/Mg] planes. These approximations are justified by comparison with the data in Figure 1, which shows a plateau in both planes at early times.
Since the SNIa channel opens with a time delay and produces mostly iron, the plateau in the [Eu/Fe] plane will be interrupted, producing the sharp dip in the [Eu/Fe] abundance (the knee) as the iron production increases. The [Eu/Mg] abundance, however, will exhibit only minor changes (due to metallicity variations or the subdominant Mg yield from SNIa), matching the observations of Fig. 1. When the background metallicity reaches ∼ c , the Eu synthesis rate will drop as collapsars cease to be formed. Iron and Magnesium production, however, undergo no such transition. This will manifest a 'knee' in both the [Eu/Fe] and [Mg/Fe] planes, with both abundance tracks dropping significantly from their previous values. However, observational evidence for such a feature is lacking. Fig. 1 does not show a second knee in the iron plane, and nor is there evidence of a significant drop in [Eu/Mg] abundance -in fact, [Eu/Mg] remains almost flat across the entire space.
This simple argument already appears to put severe constraints on the dominance of collapsars in the production of r-process elements. Note, however, that whilst these naive constraints hold in general, they can be modified, e.g. by radial migration (Tsujimoto & Baba 2019), the variations in the SFH, and the cooling of hot gas.

Analytical Model
The purpose of this paper is to explore the full possible parameter space of collapsar contributions to examine which regions are ruled out, not to point to some plausible solutions. This necessitates a rapid, streamlined model which can evaluate the chemical histories for billions of combinations of possible parameters.
To this end, we developed a Simple Analytical Chemical Evolution Model (SACEM), which makes a number of simplifying assumptions to allow a fully analytical solution. We use this model to track the enrichment of three elements representative of three key groups: the alpha elements (Mg), the iron peak (Fe) and r-process (Eu).
SACEM features many aspects common to standard GCE models, and we limit the following discussion to the points distinguishing it from other models. More details are found in Appendix C.
SACEM is a single-zone model but features two gas phases: a cold phase capable of forming stars and a hot phase into which the majority of newly synthesised material is ejected. This follows the work of Schönrich & Weinberg (2019) which showed that such thermal splitting had a drastic impact on the chemical histories inferred for r-process material, and allows us to incorporate effects such as the diffuse return of hot gas to the star-forming disc: the 'galactic fountain' of Shapiro & Field (1976).

Modelling r-Process Synthesis
The rate at which collapsars synthesise r-process material, r ( ), is dependent on the properties of the progenitor star (mass, , metallicity, , and ZAMS rotation speed, ), the r-process yield ( , , ) of each event from such a progenitor and the rate at which these progenitors undergo a collapsar event, coll = d coll /d d d , such that the exact synthesis rate is: In §2.2 we identified several problems barring us from formulating and coll , so we must instead opt for simple approximations to Eq. (2). Consistent with the approach of Si19 we assume that the yields from individual events are a constant, and that collapsars are a subset of CCSN events such that a fraction coll ( , c ) of all stars undergoing CCSN meet the criteria to become collapsars. The constant coll accounts for the mass and rotational constraints, whilst ( , c ) is the metallicity suppression function which obeys ( , c ) = 0 for > c . Therefore: Here¯ coll is the constant characteristic yield of a single collapsar event, CCSN ( ) is the rate of CCSN events at the time and ( ) is the metallicity of the star-forming gas at a time . Note that although the ejecta-mass from a single collapsar is assumed to be a constant, due to the metallicity dependence of the collapsar rate, the total synthesis rate is strongly metallicity dependent.
The formulation of reasonable functions ( ) for different processes and ( ) is discussed in detail in Appendix C, whilst the formulation of ( ) is found in §4.2.2.

Metallicity Decoupling
SACEM tracks the abundance of 3 chemical species, but not the overall metallicity. However, it is known that the yields of single elements can show strong and individual metallicity dependencies (Maeder 1992;Chieffi & Limongi 2004), making ( ) a functional of its past evolution which is impossible to evaluate analytically. For simplicity, we therefore fix = Z( ), where Z( ) is an externally imposed, monotonic function of time. As seen in Appendix C, the only explicitly metallicity dependent part of our model is the suppression rate of collapsar yields, which 'turns off' at = c . Utilising the monotonicity of Z( ), we therefore instead suppress collapsar contributions when = coll , such that: We choose Z( ) to be the function recovered from an instantiation the full RaMiCES simulation ( §4.3) which, as per SB09, replicates many aspects of the historic enrichment profile. The evolution of Z( , ) from the best-fit RaMiCES model is shown in Fig. 3. Since our data is mostly from the solar neighbourhood, we evaluate this function at the solar radius, which is well represented for 0 < < 14Gyr by the polynomial:

Deviating from the Metallicity Enrichment History
Varying the properties of our model galaxies will cause the subsequent chemical history to diverge from that used to generate Eq. (5). This is potentially highly problematic for us, as we are now imposing coll as a model parameter, instead of the physically meaningful value c . It is possible, therefore, that a model successfully replicates the r-process enrichment history of the Milky Way, and has a true cutoff metallicity such that < 0.3 ⊙ . However, because we have not updated Z for this model, we find that Z( coll ) > 0.3 ⊙ , and we would reject this model as 'unphysical'.
However, by reference to i.e. Casagrande et al. (2011);Schönrich & McMillan (2017); Haywood (2008), we have constraints on what a physically meaningful ( ) can be. Since RaMiCES (and hence Eq. 5) is calibrated specifically to reproduce this information, the observational evidence constrains how much Z( coll ) and can differ. The 'falsely rejected' models above are explicitly those which lie in tension with this data.
Any model which has coll coll ( c ) can therefore be rejected for either: • Having > 0.3 ⊙ , and therefore failing the test of §2.1 • Requiring a metallicity enrichment history which contradicts observed evidence The converse is also true: we may generate models with ( coll ) < 0.3 ⊙ < . These models would pass the test of §2.1, but are unphysical. However, due to our focus on negative inference, accepting unphysical models limits our conclusions to upper bounds on the contribution of collapsars. Of far greater concern to us is rejecting physical models, which we must ensure we do not do.

Yield Calibration
Due to the parameterisation of the stellar yields, each pathway has an undetermined constant in the form of an effective yield¯ , usually a result of the IMF-weighting of the true yield, as well as effects such as galactic-ejection fractions and remnant-lockup. A true GCE model (e.g., Portinari & Chiosi 2000) would try to derive these prefactors from first principles -for the purposes of SACEM, we instead fix these values by requiring that the curves of interest to us (those for Fe, Eu and Mg) replicate some chosen calibration values. The chosen calibration points, and their values in the nominal model, are shown in Table 1.
We note again that this strategy allows additional degrees of freedom into our model, but that since we are investigating which regions of parameter space are excluded, this in fact strengthens any conclusions we might draw: additional constraints on¯ would serve to reduce the viable parameter space, not expand it.

Action of Delays & Thermal Phasing
Central to understanding our chemical evolution models are the differing Delay-Time Distributions (DTD) of yields, and delays of injection to and freeze-out from the hot gas phase.
For the DTDs, CCSN and Collapsars occur almost immediately as Table 1. The calibration points used in the SACEM to fix the yields of CCSN, SNIa, Collapsar and NSM events. The final constrant, , has no nominal value as it the parameter we wish to investigate: the contribution of collapsars to the total r-process budget.

Observable Symbol Nominal Value
0.02 Collapsar fraction (at = today ) the lifetime of high-mass stars is negligible compared to the timescale of chemical evolution ( 10Gyr). SNIa and NSM, however, rely on the death of a previous population of stars for their formation mechanism, giving a minimum timescale of about 200Myr for SNIa and 10Myr for NSM. With reference to §2.1, we also have that the collapsar channel closes at the time coll , as the metallicity at this time is too large to allow collapsar formation. The left hand panels of Fig. 4 shows the (calibrated) rates at which the three tracked metals are synthesised within a given SACEM initialisation, demonstrating the impact that these differing timescales and cutoffs have on the associated yield functions.
In addition to the absolute rate of synthesis, Schönrich & Weinberg (2019) showed that the rate at which polluted gas becomes available for star formation also has a large impact on GCE models of rprocess synthesis, because we expect diffrent timescales for CCSN, NSM and SNIa gas-availability. For example: the majority of CCSN occur close to regions of high star formation (and hence close to the feedback 'chimneys' of Norman & Ikeuchi 1989), so we expect a large contribution of the gas to be stored into the hot reservoir, or temporarily ejected from the galaxy, only to 'fountain' back in, as per (Shapiro & Field 1976). NSM, however, have a significant time delay, and also experience natal kicks, so can be expected to be far removed both spatially and temporally from feedback chimneys, and their synthesis of high-mass, neutron rich material would lead to strong line cooling of the ejecta, making the material available for star formation much more rapidly. Microphysics such as dust-formation may also have an impact on the availability of gas for star formation.
To model these effects, SACEM places a portion of the material synthesised by process into the hot gas reservoir, with the remainder going into the hot-gas reservoir 4 . We then allow the each process material to cool at a slightly different rate such that cool = hot . Following Schönrich & Weinberg (2019), we adopt ( CCSN , Coll , NSM , SNIa ) = (0.75, 0.75, 0.4, 0.99) and = CCSN = 1Gyr as our nominal model, though we note that their fixing of the non-NSM values was somewhat arbitrary. The right hand panels of Fig. 4 show how the synthesised yields are injected, cooled and consumed by star formation over the course of galactic history.
In this class of models the 'dominant' source of Eu changes drastically with time: even though the final contribution of collapsars to the Eu budget is ∼ 10 , at early times it accounted for more then 90 per cent of the synthesised Eu, of which more than 50 per cent was found within the hot gas phase, and the collapsar-dominant The rate of synthesis of the three species within SACEM and right: the resulting mass of each element stored in the hot and cold gas phases for a calibrated SACEM instantiation chosen to generate r-process contributions ( , , ) = (0.1, 0.02, 0.88) at simulation end, and with a collapsar cutoff time coll = 1.35Gyr, corresponding to c ≈ 0.3 ⊙ . The values of ( , , ) were selected to amplify the signal of the differing behaviours, rather than generate a viable chemical history. This model would be considered 'unsuccessful' by §5.4.1 -though we note that a qualitatively identical (but harder to interpret) plot for ( , , ) = (0.02, 0.02, 0.96) would be considered successful. early time contribution was maintained even as < 0.01. Even if collapsars are negligibly responsible for r-process synthesis at late times, Fig. 4 shows that they may have been important contributors at early times. This model class would be disfavoured due to the Sn08 fingerprint indicating a monolithic r-process source, however, if we suggest that the high-[Eu/Fe] stars sampled to form the 'fingerprint' may not be representative of the bulk of stars at this early time (which did source their Eu from collapsars), this tension is alleviated.

Resolving Contribution Ambiguity
, the collapsar contribution fraction, has two distinct definitions: (i) The mass ratio of collapsar-sourced europium to the total europium within the cold, star forming gas at a time (ii) The total mass ratio of all collapsar-sourced material ever produced to the total amount of europium ever synthesised (i.e. including that which has been subsequently locked up in stars, stored in the hot gas reservoir or folded into black holes) These two definitions can diverge significantly: material produced early on is either lost to the IGM or locked up in stars/stellar remnants, so the definition (i) sets weights much more towards more recent enrichment. With that significant difference, we note that the qualitative picture and our inferred conclusions remains unaltered between the two choices and we will use option (i), defining as the collapsar contribution to the current europium in the cold gas phase. The s-process contribution, and NSM contribution, are defined similarly.

SACEM: A Recipe
With reference to the derivations in Appendix C, we set out a procedure to derive a chemical history as a function of the parameters of a given SACEM model (listed in table 2): (i) Solve differential equations Eq. (C2)-(C5) to generate a SFR: (ii) Use the SFR to generate the (uncalibrated) event rates for the four processes, using I from Eq. (C17): (iii) Solve Eq. (C8)-(C9) to find the (uncalibrated) cold-gas mass in the disc produced by each process as a function of the relevant yield, thermal and lockup parameters.
∈ {ccsn, colls, nsm, snIa} (iv) Construct models for the mass of each element within the cold gas as a function of unknown prefactors: (v) Invoke a function C which calibrates the unknown prefactors against the observed data and model inputs of Table 1: This method generates 4 analytical functions which can be used to plot Tinsley diagrams of the generated chemical history of the galaxy, and forms the core of SACEM.

Full Simulation
SACEM is designed to run quickly with minimal resources, allowing for maximal parameter searches. The penalty for this, however, was a number of potentially unpalatable approximations. In order to ensure that SACEM is not leading us astray, we also make use of a full multizone, multi-phase GCE model which captures much more physicsat the cost of orders of magnitude more computation time.
We use a modified version of of the Radial Migration with Chemical Evolution Simulation (RaMiCES) code developed by Schönrich & Binney (2009a) (SB09), with the updated parameters and inside-out disc growth developed in Schönrich & McMillan (2017) and the dual-phase NSM r-process injections of Schönrich & Weinberg (2019).
A brief description of the code and its functionality can be found in Appendix D, though the interested reader may find the above papers more complete.
For the present work, we modified the base model: we have expanded and updated the RaMiCES chemical yield grid to account for a collapsar contribution (modelled simply as a subset of CCSN events) and a low level s-process contribution to the europium synthesis. In addition, we can now track elements by source -allowing us to distinguish between NSM-origin metals and collapsar-origin.

Varying coll in Analytical Models
As an initial experiment (and to confirm the intuition developed in §4.1), we observe the effects of modifying the collapsar cutoff time, coll on SACEM chemical histories.
The investigated model was calibrated against the iron and magnesium distributions -the only constraints placed on the Europium abundance is that, at simulation end, [Eu/Fe] = 0 and ( , ) = (0.98, 0.02), i.e. be 98 per cent collapsar in origin, with no NSM contribution.
Using the cut − coll relationship of §4.2.2, we generate 6 chemical histories corresponding to cutoff metallicities of 0.1, 0.3, 0.5, 0.7, 1 and 20 ⊙ , with the final cutoff being deliberatly large so as to take place at an infinite time in the future (we refer to this model as ' c ≫ ⊙ ').
The resultant chemical histories, shown by the solid lines in Fig.5

RaMiCES & Radial Structure
This is reinforced by Fig 6, in which we also run the same set of constraints on the RaMiCES simulation, verifying that these initial results hold up in the full multi-zone model. Unlike SACEM, the models used to produce Fig. 6 are calibrated to match the early-time paths, and so we see these curves plunge below [Eu/Fe] = 0, rather than their early-time abundances shooting up.
We see that the RaMiCES models suffer from a greater differential between their early time abundances and their final values than SACEM did. This indicates that the approximations in SACEM make it more lenient with regards to gas depletion timescales, or late-time Eu abundances sourced from the hot gas phase. As a result, SACEM will provide less stringent constraints for the late Eu production than the full simulation. This works in our favour, as it means any regions of parameter space that we are able to exclude is likely to be a lower bound on the size of the excluded region.
As a multi-zoned model, RaMiCES allows us to investigate the guiding-centre radius ( ) behaviour of the abundance distributions, shown in Fig. 7. The rise in [Eu/Fe] seen from = 3 to ∼ 10kpc would be what we expect from the impact of the galactic metallicity gradient: the inner galaxy hits Z = c first, and so rprocess abundances begin to plummet first at lower . The inner radii have been Europium-deprived for longer than mid-disc radii, For = 0.5, 0.7 and 1 ⊙ , this pattern in [Eu/Fe] is interrupted by an 'arch' at > 10kpc -an effect induced by the inflow/accretion prescriptions of RaMiCES. Currently, RaMiCES tethers the composition of infalling material to the composition of the gas found at a certain galactocentric radius. The 'arch' indicates the point at which accreted material becomes the dominant driver of the cold phase metallicity -and since the abundances of the inflowing material cannot yet be reliably determined, this indicates the point at which RaMiCES cannot robustly predict the radial structure.
The behaviour inwards of the solar radius, however, is robust against the IGM prescription chosen and is a necessary consequence of the metallicity gradient of the galaxy. Though a lengthier discussion of the radial behaviour of [Eu/X] is beyond the scope of this paper, we note that the radial distribution of europium, as inferred from cepheids in Luck & Lambert (2011) is inconsistent with a such a drastic increase in [Eu/Fe] with galactic radius. This further emphasises the conclusions drawn in §3.2: there is little evidence for a large change in galactic chemistry in the recent past -any changes to the galactic chemistry must have been far enough in the past for the radial mixing of the galaxy to smooth out the emergent patterns. This could be leveraged to provid a much greater constraint on the success of our models -the models of Fig. 4 undergo their collapsar cutoff during the thick/thin disc transition, and so we would expect to see the resultant radial patterns strongly rule out these model, even if properly calibrated. The restrictions imposed on our SACEM instantiations are far from the strictest we can generate, further emphasising that our conclusions are the upper bound of collapsar contributions.

Impact of the Simple Yield Approximation
In the absence of any observational constraints on the collapsar yield function coll we follow the literature (i.e. Si19) by adopting the simplest possible form: that of Eq. (C14). One could argue that this functional form is inappropriate, or otherwise corrupts any conclu-  sions we might draw from our models. However, several factors justify that approach.
The behaviour of Figs. 5 and 6 shows that even if the yield functions are not well-approximated by a simple function, the overall range in [Eu/Fe] is a necessity that cannot be mended by introducing an additional variation into coll . Once collapsars stop forming, the [Eu/Fe] ratio drops at a fixed rate which depends only on the gas depletion timescale, and independent of any chosen form of the yield function. The only way to ensure that [Eu/Fe] ≈ 0 at simulation end is by adding more Eu before collapsars die off. However, because collapsars are only active for a short time, no matter the functional form of the yield, increasing the Eu abundance pushes [Eu/Fe] above our 'reasonable domain', and the model would be ruled a failure. We demonstrate this principle in Fig. 8, in which we replace our simple yield function with a positive random-walk function +1 = | + | where is a random number in the range [-1,1], allowing for arbitrary functional forms of coll ( ). We see that no matter the functional form, after the collapsar suppresion begins, all of the yields converge rapidly to a single path through [Fe/H]-[Eu/Fe] space. This path is determined solely by the required amount of Eu in the galaxy at the time that collapsars are suppressed, and hence is independent of the collapsar yield function: i.e., the models care how much Eu there is in the galaxy at coll , not how it got there.
The post-turnoff behaviour of these models would therefore be very well represented by a 'simple'-type model which was tuned to produce the same Eu mass at collapsar turnoff as the varying model. For a given set of model parameters if a simple-type yield would cause the model to be ruled incompatible with the observed evidence at any time after collapsar suppression, then so would the corresponding model with a varying yield function.
By eliminating model classes based on their simple-yield approximations, we cannot be falsely eliminating any physically meaningful models, so we conclude that our results are robust against the impacts of the approximation of the simple form of coll .

Expanding the Search
In these initial forays we were altering a single parameter in a single 'best fit' model, so one could argue that the search missed the right combination of parameters that allows for a high collapsar contribution despite a reasonable cut-off metallicity.
The challenge is that such a set of parameters must not simultaneously make, for example, the [Mg/Fe] histories unrealistic: we must ensure that any alterations produce a simultaneously realistic model in multiple chemical planes.

'Success' & Viable Models
Following up on our summary in §3.4 we now define the criteria for a viable chemical evolution model. We emphasise that, in line with our efforts to falsify model classes rather than find a best fit model, the criteria developed here do not imply that a model is fully physical if it fulfils the criteria, but models that breach the criteria are clearly in contradiction to the empirical evidence.
We define a galactic chemical history as 'viable' if the resulting Tinsley curves for [Mg/Fe], [Eu/Fe] and [Eu/Mg] lie between the two black two curves in the corresponding panels of Fig.1. Whilst a good model should also reproduce the observed distributions not just pass within its range, in keeping with our approach of negative inference, we instead choose to have room to be generous and yet still constrain the models. All three planes must simultaneously meet this criterion for the model to count as successful.
We make the intentional choice that the constraints only apply for [Fe/H] > −1.5. Though this weakens our constraints, it gives us numerous advantages: we limit ourselves to the strictly non-stochastic regime, avoid halo contamination and as per 5.3, we avoid eliminating classes of models which have high variability at low metallicities.
In addition, we study the effect of introducing a constraint on the time evolution of the models in concordance with the observational constraint in Fig.2 where this is measured over the final 2 Gyrs of the simulation, and is in agreement with Fig. 2. Finally, a model galaxy must be capable of a sustained rate of star formation, even at late times. In GCE models it is common to constrain this through the star formation efficiency, represented by (up to a factor of order sfr ) the quantity c / * . For the Milky Way this value can be measured (i.e., McKee et al. 2015), indicating that we should constrain this value to ∼ 0.1.
We reemphasise that we have left the acceptance criteria intentionally generous: we want to find which models are excluded by failing to meet even these lax criteria.

Variable Selection for SACEM
We perform a Monte-Carlo exploration of the high-dimension parameter space by randomly and independently selecting the model parameters (except , the collapsar contribution ratio, and coll , the collapsar cutoff time) from between bounds determined by the imposed constraints (see Table 2). For each random set of parameters, we then perform a sweep ofcoll space, such that each instantiation is evaluated uniformly across this space. Note that the enforcement of the temporal gradient constraint (Eq. 19) plays a special role, as we perform all simulation runs both with and without this constraint.
We call the different constraint sets in Table 2 the Unconstrained set (denoted U), Weak (W) and Viable (V). As indicated by the name, in Set U, parameters are allowed to take almost any value, with essentially no constraint or regard for their physical reality. Naturally, the 'unconstrained' set is not truly unconstrained. Instead we mean that the boundaries imposed reflect our desire to search 'interesting' regions of parameter space within finite computation time, and do not meaningfully eliminate any kinds of physically interesting galaxies we might care to consider. We include this set to act as a null hypothesis -that without physical constraints, models can be made to fit to the data -and hence that any region of parameter space that is eliminated is due to the imposition of physical constraints.
The Weak set of parameters (W) adds an amount of physical intuition -the order of magnitude of the galactic mass is determined, cooling timescales are of an order close to what we might expect, and so on. These constraints encapsulate the behaviour of most galaxies, but do not uniquely identify the Milky Way. In contrast, the Viable set imposes the minimum requirements to match the properties of the Milky Way.
A fourth variant, the Mixed model (M), was also studied, which uses the Viable constraints for all parameters except those relating to the SFR, which are treated as Unconstrained. The mixed set forms the basis of the discussion in §7.

SACEM Search
For each of the constraint sets of §5.5 we randomly generated 10 7 parameter sets, which were then each evaluated across a 101×101 grid in − coll space, for a total of ∼ 10 11 models.
Due to the high density of generated models, we may therefore be confident that the regions in which no models could be found (the 'exclusion zone') are genuine forbidden regions of parameter space. This inference is strengthened if the exclusion zone are contiguous: whilst statistical fluctuations may alter the position of the boundary, the contiguous nature implies that regions a significant distance from this boundary are robustly excluded.

Without Temporal Gradient Constraints
We first consider those models which did not use the condition of Eq. (19) to evaluate the success of a model. The density of successful models is shown in Fig. 9, with regions coloured black denoting the points where no successful model could be found.
As expected, for the Unconstrained (U) and Mixed (M) constraints, we are able to find allowed models for all values of the collapsar contribution and the cutoff time coll . However, the density variation for models with large-but low coll (i.e. lower right of each panel) shows that such models are disfavoured even with these lax constraints: in both the Unconstrained and Mixed investigations, NSM-dominated models (left of each panel) were favoured by a factor ∼ 500 over the corresponding high--lowcoll models.
With the introduction of the only Weak constraints, we see the formation of an exclusion zone in the high-, lowcoll region: the Weak constraints fail to find any models which are collapsar dominated (Ω ∼ 1) for coll < 6Gyr, and limits the Milky Way to Ω > 0.65 for coll < 3Gyr. The exclusion zone is not qualitatively changed by the introduction of the Viable constraints, but the size increases such that collapsar dominated models are prohibited for coll < 7 Gyr, Ω > 0.5 is prohibited for coll < 4Gyr and Ω > 0.2 is prohibted for coll < 2Gyr.
This immediately brings us into direct tension with the theoretical models of collapsars. As per §2.1, no theoretical models yet allow for events above = 0.3 ⊙ , but following the discussion in §4.2.2 and Figure 9. The results of the SACEM gridsearch for models where no gradient check was performed, and where the fraction of europium, , was measured from the present-day cold gas reservoir. Black regions are those in which no successful models (as determined by the success conditions of §5.4.1) could be found. Note that the coll axis does not extend exactly to coll = 0 (the region where collapsars are never active), as ≠ 0 is nonsensical if collapsars are inactive for all history. Therefore coll ≥ 0.16Gyr. the work of i.e. Casagrande et al. (2011), we know the Milky Way was rapidly enriched early in its lifetime. In other words, ∼ 2Gyr is the approximate time coordinate that we should associate with the 'maximum theoretical collapsar time'.
These initial results therefore strongly indicate that, under Weak and Viable constraints, the r-process enrichment of the Milky Way cannot be dominated by collapsars events constrained to occur before = 0.3 ⊙ . However, in §4.2.2, we cautioned against inferences which rely on coupling coll and c too tightly, as there is a small amount of flexibility in the relationship due to our decoupling of the background metallicity evolution from the properties of galaxies. The Viable exclusion zone in Fig. 9 has a boundary at ∼ 4Gyr, so between the Monte Carlo search and the decoupling of c and coll , it is plausible that our results might still allow a collapsar dominated galaxy which we failed to detect: more work is yet needed.

Temporally Constrained Models
In Fig. 10 we examine the same set of constraints as Fig. 9, but with the additional constraint that successful models must obey d[Eu/Fe] d ≤ 0.01 dex Gyr −1 across the final 2Gyr of evolution. We note from Fig.  2 that this is still a generous constraint.
We see that, once more, models which use Unconstrained and Mixed constraints excluded no regions of parameter space, though as before they favoured low-or highcoll models by more than a factor of 10 2 and 10 3 respectively.
However, for Weak and Viable models, we see that large regions of space have been declared non-viable. Comparing models W_noGrad and W_Grad (the lower left panels of Figs. 9 and 10 respectively) we see that the addition of the temporal gradient constraint has expanded the exclusion zone. The region > 0.6 (previously excluded for coll <6Gyr) is now excluded for all times coll < 12Gyr. This makes the exclusion zone in this region extremely resistant to both statistical errors and errors arising from the imposed metallicity function, as such errors would need to cause errors on the order of 5Gyr, which we do not consider a reasonable expectation.
Similarly, with the addition of the temporal constraints, all successful Viable models were limited to < 0.4 for coll < 12Gyr and < 0.3 for coll < 4Gyr. This eliminates all models for which collapsars are either the majority or dominant source of r-process nucleosynthesis, as there is substantial theoretical and observational evidence that collapsars cannot have been occuring more recently than 2Gyr ago.
Again, we emphasise that the vertical expansion of the exclusion zone observed in models W_Grad and V_Grad makes these conclusions extremely robust to statistical and − coll induced errors. We are therefore condfident that these regions are truly excluded portions of parameter space.

RaMiCES Search
We performed a similar search on the RaMiCES simulation, though due to computational limitations, the number of models is much smaller. In addition, RaMiCES is more physically motivated than the semi-empirical approach of SACEM, such that there are fewer parameters which we may alter without undermining some other observable property of our galaxy. In particular, the SFR (which is coupled to gas infalls and radial gas motion) is chosen to replicate features of the solar neighbourhood (see SB09). The remaining free parameters which are not fixed to replicate observable properties of the galaxy are therefore: • The hot-gas injection fraction CCSN   Because the RaMiCES model is inherently much more tightly constrained, we widened the success condition to prevent overfitting of the model. The only success condition applied to the RaMiCES models is that, for the solar radius, the final [Eu/Fe] value must lie within the 'viable domain' of Fig. 1  Note that, unlike SACEM, this simulation does use metallicity for the collapsar cutoff, rather than a metallicity-inferred time. We also note that due to the properties of the RaMiCES simulation, we cannot target a specific final collapsar-contribution fraction ; we have to generate models with a given set of europium yields and determine at simulation end. Hence, unlike the models of the previous section, the models are not produced on a uniform grid of , and due to the computational constraints it proved somewhat difficult to even generate a meaningful number of models with ∼ 1 at low c . In particular, we were not able to generate a single model for which c = 0.1 ⊙ but > 0.85. However, given the contiguous (and large) nature of the final exclusion zone, we do not expect this to impact our final conclusion.
We ran 3634 iterations of the RaMiCES simulation, with the subset of variables drawn from the Weak sample shown in Table 2. The results of this search are plotted in Fig. 11. Despite the fact that the random parameters were generated with the Weak model, with no explicitly included gradient consideration and with an even more generous acceptance criteria, RaMiCES produced a grid with an even larger exclusion zone than the temporally-contrained Viable model of SACEM.
It seems that the initial results we inferred from Fig. 6 -that a multi-zone model with more physically motivated parameters suffers a significantly larger drop in [Eu/Fe] than the corresponding SACEM model -continues to hold, thereby indicating that the SACEM results should be interpreted as an upper bound on the maximum collapsar contribution.

PROPERTIES OF SUCCESSFUL MODELS
Figs. 9 and 10 demonstrate that the metallicity dependence of collapsar events necessarily limits them to sub-dominant contributors to the galactic r-process budget, yet it is more instructive to exmaine in more detail those models deemed "successful".
SACEM is unique in its low computational footprint, so we can for the first time systematically scan the full parameter space. This allowed us to entertain non-physical sets of parameters in our search for the minimum required set of constraints: the models generated from the Unconstrained and Mixed sets, in particular.
This was, primarily, an aesthetic choice: we wished to impose as few additional constraints onto our models as possible. However, by examining the behaviour of the global parameters from the lessconstrained models within the Ω − coll plane, we may infer how the successful models fulfill the imposed requirements. This provides three benefits: i) comparing how past studies have been able to match chemical observations, and if this lies in tension with other observables ii) examinations into the limitation to these approaches, i.e. how far we may 'tweak' the parameters of our model whilst remaining physically viable and iii) examining if there are other ways for a model to fit the data beyond the canonical understanding.

Star Formation
The most straightforward way for a model to produce a required chemical history within the specified bounds is with an extremely tightly controlled star-formation rate. Within SACEM, if a model has a highly-peaked early-time SFR, then the large collapsar population could generate a vast amount of Eu over a short period of time. In such Figure 13. As with Fig.12 but with sfr (2)/ sfr (14): a measure of how quickly the star formation in the galaxy dies down. models, a true GCE model would also accumulate an equally vast number of metals during this period due to the high SFR, such that collapsars should shut off. However, because the imposed function Z( ) is insensitive to the particular parameters of our model, the 'true metallicity' of the model diverges significantly from the nominally assumed model, and so collapsars will continue contributing long past the time they should have died out.
This therefore enables models with even tiny values of coll to be considered viable, where a more physically coupled SFR and metallicity evolution would discard these models. We suspect that the vast majority of coll − space accepted under the Mixed regime, but rejected under Weak and Viable models are achieving their success through this method.
Figs. 12 and 13 justify this claim, they respectively depict the ratio of the initial star formation rate and after 2Gyr: sfr (0)/ sfr (2Gyr) and the corresponding ratio between 2Gyr and the end of the simulation, sfr (2Gyr)/ sfr (14Gyr), for two classes of models. The properties of the models in the high--lowcoll region of M_Grad are striking: the search preferentially found models where the initial SFR was more than 2 orders of magnitude greater at simulation start than it was just a few gigayears later, and where the mid-time SFR was on average a factor of 6 greater than the SFR at simulation end, such that there is a total of more than 3 orders of magnitude between the = 0 SFR and the final time SFR. These ratios are orders of magnitude outside the constraints measured for the star formation history of the MW, which typically show a decline in SFR by less or equal to an order of magnitude (Mor et al. 2019;Schönrich & Binney 2009b;Aumer & Binney 2009).
By comparison, the V_Grad model shows no particular bias in either of Figs. 12 or 13, indicating that large-scale constraints such as the total mass of the Milky Way are reasonable methods to eliminate these models.

Hot Gas
The inclusion of a thermal gas phase is a relatively rare feature of GCE models, though recent work has highlighted its importance (Khoperskov et al. 2021;Schönrich & Weinberg 2019). Within SACEM and RaMiCES we implement the thermal phase in a relatively crude fashion: a distinct hot and cold phase, with the hot phase being fed with a fraction of the enriched gas, and then 'cooling' into the cold phase with a decay frequency such that cool = hot . Fig. 14 shows the variation of the parameter NSM in the successful models, which displays an interesting pattern, which we argue both vindicates the inclusion of the hot phase, and is central to our heuristic understanding of the interplay between NSM events and collapsars.
At very low Ω (≪ 0.1), this value is favoured to be low, around 0.6 or below (for comparison, the corresponding CCSN, Collapsar and SNIa values exceed 0.9 almost everywhere) -this was argued for in Schönrich & Weinberg (2019) as a way for the thermal properties to permit NSM to be viable contributors at early times, due to the correspondingly higher immediate contribution to the star forming reservoir. As collapsars become more prominent at early times (higher Ω and lower coll ), it is no longer necessary to invoke the arguments of SW19 to explain the early time [Eu/Fe] values, and hence the need for low values of NSM to meet the chemical criteria at [Fe/H] = -1.5 vanishes, leading to the observed increase in NSM . At even higher values of Ω, the value of NSM decreases again -notably to a value of 0.65, which is the exact midpoint of the permitted range: hence we are seeing ambivalence to the value of NSM when determining if a model is successful. However, for Ω < 0.7 there are clear signs that the thermal properties of NSM are very important for determining the success or failure of a model.
Given that the thermal properties have been shown to be important, we now consider the case where the cooling rate is small. If the hot-gas injection fractions are non trivial, then a large portion of the enriched gas can be secreted away in the hot gas phase, preventing over-enrichment at early times, and providing a 'source' of enriched gas into the star-forming phase long after the gas itself was actually produced. Whilst this is desirable to an extent -it is reasonable to expect exactly this to happen - Fig. 15 shows that the successful models with a large collapsar contribution were almost overwhelmingly those which abused this property to the extreme. Figure 15. The mean ratio cool , the hot gas cooling rate, for sucessful models constrained with M_Grad. Note that the midpoint of the permitted range for cool is 5 -models close to this value are likely to be insensitive to the value of . Figure 16. The mean value of mod , the lockup-modification factor, for successful models at each point in coll − space for the prototype simulations in which mod is bounded by 0 < mod ≤ 1. We note that the black strip at coll = 0 is the nonsensical combination coll = 0, > 0. As per Fig. 9 this was omitted it all other plots.
In Fig. 15 we see that for > 0.3 and coll < 13Gyr, the mean value of is ≈ 0.077 Gyr −1 , an order of magnitude below the commonly used value of ∼ 1 Gyr −1 . These models cool almost no gas into the cold gas phase, allowing the hot gas phase to become hyper-enriched relative to the cold gas phase. This shows that these successful models favoured extreme thermal fractionation as an explanation for the chemical history of the galaxy -whilst we do not doubt that a hot gas population is important for understanding GCE, the cooling rates indicated here strongly suggest that these models are highly unphysical.

Lockup Modification
The lockup modification factor, mod , encodes the rate at which synthesised material is removed from the cold gas reservoir by star formation, such that the lockup rate is ∝ mod sfr rather than the zeroth-order approximation ∝ sfr . mod differs from 1 since, though normal CCSN events may not synthesise new Eu material, they can recycle pollutant metals back to the ISM, thereby reducing the effective lockup rate of synthesised material.
We must be aware, however, of the extreme case mod ≈ 0 in which no synthesised material can be locked up. Star formation would continue to deplete the unenriched gas mass, such that the stars are preferentially forming from primordial material: the cold gas phase becomes chemically fractionated, and becomes ever more enriched without any additional synthesis events.
Collapsar models with coll < 4Gyr could exploit this property to 'hide' a reservoir of Eu which persists until late times, when by all physical reasoning, it should have been depleted by the lockup of the continual star formation. This was made evident in a prototype set of simulations (termed the Old-Style simulations), in which the lower bound on mod was set to be 0 for all simulations. Fig. 16 shows the behabiour of mod of successful models in the prototype simulation. We note that the prototype simulations differ from the final ones in several ways that make a direct comparison difficult, but we can see that the lower right section of the Fig. 16 is overwhelmingly dominated by models for which mod ∼ 0.001, the case where there is no chemical lockup. A similar pattern was observed for models U, W and V: an overwhelming bias towards extremely small values of mod in the region of high-and lowcoll .
Such models are evidently unphysical, and so in the final suite of simulations, we bounded mod from below by 0.3. This value was motivated by comparison with a reasonable IMF, noting that ∫ 1 0 ( )d > 0.35 for all commonly used IMFs (Salpeter 1955;Chabrier 2003). Since stars with = ⊙ have lifetimes ∼10Gyrs, at least this fraction of the gas must be locked-up for long timescales.
When bounded in this fashion, the bias in the values of mod vanished almost entirely, and so we surmise that we have closed off this unphysical route to achieving 'success'. We note that some values near the boundary of the exclusion zone of both W_Grad and V_Grad did show a small bias towards smaller values of mod , indicating that if we were to impose stricter and more physical constraints on the value of this parameter, these models would similarly be eliminated. However, in the spirit of our attempt to find the minimum possible set of constraints, we leave these potentially problematic models unchallenged, as the size of the exclusion zone is already sufficient to draw our conclusions.

Improperly Coupled Models
In the above discussion, we saw a number of traps which models can fall into: failing to properly lock up their materials, poor treatment of hot gas phases and associated cooling rates and torturing the SFR until it allows you to replicate your desired features. In utilising these 'traps', the models seemingly satisfy all chemical constraints which were placed upon them. However, on closer inspection these models were only able to reproduce the chemical properties due to an unphysical assumption elsewhere in the model.
The general theme of these assumptions was that they allowed the chemical reservoirs to become separated or fractionated in some way, such that the evolution of the SFR, cooling, chemical enrichment and subsequent lockup of the reservoirs was not functioning properly.
We note, for example, that Si19 use an SFR which is decoupled from the present gaseous or stellar mass in their model, which we suggest falls into a similar camp of improperly coupled models, and explains why their findings are in such strong tension to our own, despite seeming to meet all of the observable chemical criteria.
We also suggest that, due to the expansive sampling of parameter space, we would have been able to notice if some unusual combina-tion of physical properties were able to replicate both the chemical data, and not fall into one of the four 'traps' outlined above -in fact, we observed no such signal. Given that this is the case, it must be true that any successful GCE model which replicates a collapsar-dominated galaxy must strongly deviate from the physics encapsualted in the core equations of Appendix C.

CONCLUSIONS
In this paper, we have performed the first comprehensive chemical evolution study which examines the multi-dimensional parameterspace associated with the origin and evolution of galactic r-process material. In this extensive analysis, we could find no viable model with collapsars as the dominant source for today's r-process element budget.
In this work, we have introduced our newly developed Simple Analytical Chemical Evolution Model (SACEM). SACEM is an analytical framework, which incorporates the relevant physics (star formation histories, inflow, outflow of gas, yields to both a hot and cold gas phase, cooling of material from the hot phase, star-forming ISM, and different temporal and thermal properties for different sources of yields), but at the same time has run-times of fractions of a second, i.e. orders of magnitude faster than existing chemical evolution codes. Although SACEM relies heavily on some simplifying approximations (namely single-zone space with instantaneous mixing, empirical fixing of yields, and does not consider the lifetime of stellar populations), we found good agreement with the full chemodynamical models from Schönrich & Binney (2009a), which does not make such approximations.
Where there are divergences between the models, we found that SACEM was more generous to collapsar-dominated models than the SB09work, though we note that both SACEM and RaMiCES do not directly consider tertiary sources of r-process material (such as the i-process or magneto-rotational supernovae), the impact of dwarf galaxy accretion, and assume a simple form of the collapsar yields. Although these omissions might limit the generality of our results, we have justified their long-term impacts on the abundance patterns as negligible, or already encapsulated in part by features of our models. Hence, our results are robust against these approximations.
A central problem holding chemical evolution studies back has been a reliance on costly models in a high-dimensional parameter space, which has forced prior studies to operate with exploratory modelling of a small number of models. SACEM's performance allowed us to run > 10 11 models, mapping out the full parameter space of r-process chemical evolution with both collapsars and neutron star mergers, and allowing us to pursue an entirely different strategy: instead of trying to find models that match some observational constraints, we drew up a full set of "minimal consensus" observational constraints which models must replicate, and look for those models which fail to reach even these lax conditions: rough 'bounding boxes' that our chemical tracks in [X/Fe]-[Fe/H] must pass through (Fig. 1), chemical equilibrium imposed across the final 2Gyr of evolution, and a sustained rate of star formation. Our search through parameter space was bounded by imposed conditions on the allowed parameter values -the one which best represented the Milky Way (whilst still allowing for large degrees of variation) we termed the Viable set of constraints. We also explored the parameter space with unphysically lax constraints on parameter values (the Unconstrained and Mixed sets) -comparing these with the Viable constraint set reveals and analysis how classes of models with dominant r-process contributions discussed in the literature appeared to satisfy observational constraints. However, with the Viable constraints, we have found that: (i) no SACEM model could be found where collapsars contribute more than 30 per cent of the modern r-process budget, as long as collapsars were suppressed as in MacFadyen & Woosley (1999). Neutron Star Mergers were always required to be dominant (Fig. 10).
(ii) A significant collapsar contribution at early times was not eliminated: Many models in which NSM are responsible for > 99 per cent the Eu abundance at late times had > 50 per cent collapsar contributions at < 1Gyr (Fig. 4).
(iii) The RaMiCES code shows that the remaining parameter space allowed by SACEM still contains models that are in stark contradiction with the data. In particular the metallicity dependent cut-off can introduce a radial [Eu/Fe] increase in the galactic disc (Fig. 7), which starkly contrasts with observations: our limit of <30 per cent is likely still too high, and can be refined further.
We deliberately chose constraints which were overly generous, and our results should be seen as the maximum possible contribution of collapsars to the modern r-process budget. We leave further discussion regarding how far the constraints can be pushed and how far the allowed parameter space can be further shrunk to future studies, preferring to keep our argument simple: we have shown that even a minimal set of constraints permits no models with collapsars as dominant source of r-process elements, and thus leaving by exclusion (Sneden et al. 2008) only neutron star mergers as a dominant source.
The usage of a direct proxy might be called into question due to hints of differing trends between Eu, Gd and Dy (Guiglion et al. 2018), and the breakdown of the 'common fingerprint' for lighter rprocess elements (see Sn08). Whilst the effects of the assumptions of Local Thermodynamic Equilibrium have been studied in Europium (i.e. Zhao et al. 2016), such corrections have not been calculated for Gd and Dy. It is therefore plausible that the reported differences in the trends between Eu, Dy and Gd vanish upon full consideration of NLTE and 3D atmospheric modelling.
On balance, Eu serves as a convenient proxy for r-process elements.

APPENDIX B: UNDERSTANDING FIG. 1
We present a brief description of how the standard models of GCE predict and explain the behaviour exhibited in Fig. 1 Here ( ) is the Initial Mass Function (IMF) and X ( , ) is the net amount of element synthesised by a star with progenitor mass and metallicity . The integrals are bounded from below by min ( ), the minimum mass star which has reached the end of its lifetime at time . Under the approximation that the variation in min ( ) and ( ) does not alter the value of the integral over short timescales, the time dependence drops out, leaving a constant abundance ratio at early times: For a similar discussion, see Weinberg et al. (2019). As more stellar evolution is allowed to occur, SNIa events can kick in (initially being prohibited by longer progenitor lifetimes and subsequent inspiral or accretion phases). SNIa heavily favour the synthesis of iron-peak elements over elements such as magnesium (Iwamoto et al. 1999), so we will see a decrease in all [X/Fe] planes where X does not have significant SNIa production -this is the 'knee' seen at [Fe/H] ∼ -1. This simple outline (neglecting confounding factors such as the thermal phases of the ISM or metal loss from the galaxy) covers the main patterns seen in the chemical evolution of [Mg/Fe].
Naïvely, it is surprising that [Eu/Fe] in Fig. 1 behaves similarly to the canonical picture of [Mg/Fe], as we do not expect CCSN to contribute significantly to Europium synthesis (indeed, this is why europium was chosen as a tracer). This has often been used as an argument in favour of Collapsars being the dominant r-process source.

C1 Star Formation
Our analytical model, SACEM, derives its time-dependent star formation rate from a physically motivated model of the galaxy and the accretion and heating of three gas reservoirs ( , the cold gas, ℎ , the hot gas, and * the mass locked up in stars). We use a Kennicut-Schmidt style model for star formation: stars form at a rate given by * , = SFR , and we use a simple 'exponential death' model for stars returning their material back to the ISM, such that * , = * (note that this is solely for the purposes of the SFR and does not alter the chemical evolution. A more complex prescription could be used, but the final result would be equally replicated by altering sfr or , and therefore would add nothing tot he model except complexity). The returned material is split into the two gas reservoirs, with a fraction ℎ going to the hot reservoir and the remainder becoming immediately available for star formation in the cold reservoir.
The hot reservoir decays into the cold reservoir with a characteristic frequency cool , but we also include a mechanism for stellar feedback: when stars of mass form from the cold gas, an additional amount of cold gas is shifted into the hot reservoir. We initially assume the the galaxy is composed only of cold gas ( ( = 0) = 0 , ℎ (0) = * (0) = 0). Subsequent infall from the IGM is parameterised by exponential laws: Where the free parameters { } and { } = 1/ set the infalling mass and timescales respectively. Together, this produces the following coupled system of differential equations: This can be analytically solved for , and hence the star formation rate SFR ( ) = SFR ( ). Because of the linearised assumptions we have made, the solution is expressible in terms of a sum of exponential terms.

C2 Elemental Synthesis & Return of Metals
We are chiefly interested in the chemical compostion of the cold gas reservoir at any given time, as this determines today's observed stellar surface abundances (with minor modifications due to i.e. dredge up or gravitational settling).
If a nucleosynthesis pathway produces an amount , ( ) of element at time , then the amount of due to present in the cold gas is given by , the corresponding amount for the hot reservoir is given by ℎ . They are linked via: The final terms in these equations arise due to star formation, which removes a fractional amount of the element from the cold gas reservoir and either heats it up through stellar feedback, or locks it up in stars. This simplifies to: = (1 − ℎ, ) , ( ) + ℎ − (1 + ) sfr mod (C8) ℎ = ℎ, , ( ) − ℎ + sfr mod (C9) Here mod has been introduced as a 'lockup modification factor', such that the lockup rate is proportional to mod sfr , instead of just sfr . This modification is introduced to allow for the fact that is the rate at which new material is synthesised. Since stars are formed from polluted gas, as long as they do not destroy the material, they can release metals which they did not synthesise. If mod < 1, therefore, we reduce the rate at which material is being locked up by mimicking the recycling of previously synthesised material.

C3 Yield Functions
To produce the yield functions, , , we invoke a Delay Time Distribution (DTD). This function, Ψ ( ) gives the probability of a star undergoing stellar death a time after it was formed 5 . The mass-rate of events (i.e. the stellar mass loss rate through channel ) occuring at a time is therefore given by: Swapping the integration variable ′ = − , it follows that the yield from event is given by: Where ( ) is the Initial Mass Function (IMF), , ( , ) is the gross yield of from a star of mass and initial metallicty dying through process , and cg ( ) is the cold-gas metallicity at a time . With Eq. (C11) in hand, we are able to derive three equations for the cases of CCSN events, collapsars and delayed/inspiral events such as SNIa and NSMas follows:

C3.1 Core Collapse Supernovae
Eq. (C11) simplifies for the case of CCSN from high-mass progenitors. Such CCSN occur at the end of a lifetime , such that the DTD becomes a Dirac delta function: Ψ( ) = ( − ). In addition, for stars massive enough to go CCSN, this lifetime is short compared to the timescale over which galactic properties change, such that we may approximate ≈ 0: In the case where the yields are independent of the metallicty (which, by comparison with yield grids such as Chieffi & Limongi (2004); Limongi & Chieffi (2018), we see holds to a good approximaiton for both Fe and Mg), we may compute the integral over to find the characteristic yield function of via ,¯ , , giving the synthesis rate via CCSN as:

C3.3 Delayed Yields
Finally, we consider the yields of SNIa and Neutron Star Margers. These events do not occur uniquely at the end of a stellar lifetime. Instead, after the stellar lifetime has passed, there exists a period of probabilistic decay, whilst the system continues to evolve until finally the progenitors inspiral (for double degenerate SNIa and NSM events), or accrete enough matter from their companion (for single degenerate SNIa). In addition, there exists a non-trivial time delay before the first events can start occuring. Whilst the common DTD for SNIa is typically given as ∝ −1 , in order to continue our ability to easily analytically integrate them, we follow the work of SB09in using an exponential DTD for SNIa events -for a discussion of the validity of this approach, see Weinberg et al. (2017). Hence: Assuming metallicity independence, we find: The Iwamoto et al. (1999) W70 is a common metallicity-independent model for SNIa yields, though recent efforts such as have attempted to account for progenitor metallicity, Travaglio et al. (2005) for example shows that the Fe yields are to be altered by less than 6 per cent between 0.1Z ⊙ and ⊙ , such that we consider the metallicityindependent model a good approximation.

APPENDIX D: SIMULATION MODEL
RaMiCES models the the galaxy as a series of concentric rings. Each annulus (of ∼0.2kpc in width) contains a number of gas reservoirs -each of which has an independent chemical composition and is assumed to be chemically well mixed -in addition to containing stellar and stellar-remnant populations.
As per the prescription of SB09(and unlike the analytical model in 4.2) the composition of each annulus does not evolve independently: RaMiCES incorporates both radial migration of stars due to resonant scattering ('churning') and the oscillation of stars around their guiding centres due to epicyclic motion ('blurring'), allowing for stars (and hence the chemicals they produce upon death) to migrate away from their place of birth.
In addition to the migration of stars, galaxies require a steady inflow of fresh gas in order to sustain sufficient star formation rates and avoid depletion (Chiosi 1980), which must in turn drive radial flows of gas within the galaxy. We use the formalism of Bilitewski & Schönrich (2012) to account for the angular momentum balance. The material accreted from the IGM is not pristine, but otherwise its composition is poorly constrained. Schönrich & McMillan (2017) approximates the inflow composition using the abundance distribution of a ring in the mid-outer disc, but notes that this only materially affects the outer disc's metallicity.
The chemical yields from exploding stars are taken from a number of sources in order to cover the wide range of mass and metallicity required: we produce a compiled grid from the data of Marigo (2001); Chieffi & Limongi (2004); Maeder (1992) and data retrieved from the ORFEO database of Limongi & Chieffi (2008). The stellar lifetimes as a function of mass and metallicity are extracted from the BaStI database of Pietrinferni et al. (2004).

D1 r-Process Yields
The r-process synthesis contributions from CCSN and Collapsars are added into the usual CCSN yield network. We have allowed CCSN ( , ) = , a small, constant level of synthesis to arise from CCSN throughout history. This parameter is calibrated to give a 2-5 per cent contribution at simulation end. We then add a collapsar contribution derived from Eq. (3). Unlike SACEM, the Collapsar yields have the same thermal properties (i.e. distribution between hot and cold gas phase and cooling timescale) as standard CCSN gas, a consequence of the incorporation into the standard yield network.
NSM contributions are modelled as in SW19 -they are treated near-identically to SNIa events, with the exception that they have a shorter initial delay time (their projenitors are higher mass objects, and so have shorter lifetimes), and a lower hot-gas injection fraction -that is, a larger fraction of gas from NSM events is immediately available for star formation than CCSN or SNIa events, as justified in §4.2.5.
As with SACEM each of these pathways -NSM events, collapsars and small s-process contribution -has an undetermined prefactor, these are chosen either to reproduce the observational data, or to target some ideal collapsar/NSM/s-process fraction at simulation end (or some combination of both), though due to the numerical complexity of the simulation, the values have to be tuned by hand, rather than using simple analytical constraints. This paper has been typeset from a T E X/L A T E X file prepared by the author.