Can Neutron-Star Mergers Explain the r-process Enrichment in Globular Clusters?
Michael Zevin Kyle Kremer Daniel M. Siegel Scott Coughlin Benny T. -H. Tsang Christopher P. L. Berry Vicky Kalogera




zevin@u.northwestern.edu
Department of Physics and Astronomy, Northwestern University, 2145 Sheridan Road, Evanston, IL 60208, USA
Center for Interdisciplinary Exploration and Research in Astrophysics (CIERA), 2145 Sheridan Road, Evanston, IL 60208, USA



Department of Physics and Astronomy, Northwestern University, 2145 Sheridan Road, Evanston, IL 60208, USA
Center for Interdisciplinary Exploration and Research in Astrophysics (CIERA), 2145 Sheridan Road, Evanston, IL 60208, USA



Perimeter Institute for Theoretical Physics, Waterloo, Ontario, Canada, N2L 2Y5
Department of Physics, University of Guelph, Guelph, Ontario, Canada, N1G 2W1



Center for Interdisciplinary Exploration and Research in Astrophysics (CIERA), 2145 Sheridan Road, Evanston, IL 60208, USA



Kavli Institute for Theoretical Physics, University of California, Santa Barbara, CA 93106, USA



Center for Interdisciplinary Exploration and Research in Astrophysics (CIERA), 2145 Sheridan Road, Evanston, IL 60208, USA



Department of Physics and Astronomy, Northwestern University, 2145 Sheridan Road, Evanston, IL 60208, USA
Center for Interdisciplinary Exploration and Research in Astrophysics (CIERA), 2145 Sheridan Road, Evanston, IL 60208, USA
CIFAR Fellow





Star-to-star dispersion of r-process elements has been observed in a significant number of old, metal-poor globular clusters (GCs). 
We investigate early-time neutron-star mergers as the mechanism for this enrichment. 
Through both numerical modeling and analytical arguments, we show that neutron-star mergers cannot be induced through dynamical interactions early in the history of the cluster, even when the most liberal assumptions about neutron-star segregation are assumed. 
Therefore, if neutron-star mergers are the primary mechanism for r-process dispersion in GCs, they likely result from the evolution of isolated, primordial binaries in the clusters. 
Through population modeling of double neutron-star progenitors, we find that most enrichment candidates are fast-merging systems that undergo a phase of mass transfer involving a naked He-star donor. 
Only models where a significant number of double neutron-star progenitors proceed through this evolutionary phase give rise to moderate fractions of GCs with enrichment; under various assumptions for the initial properties of GCs, a neutron-star merger with the potential for enrichment will occur in ∼ 15-60% (∼ 30-90%) of GCs if this phase of mass transfer proceeds stably (unstably). 
The strong anti-correlation between the pre-supernova orbital separation and post-supernova systemic velocity due to mass loss in the supernova leads to efficient ejection of most enrichment candidates from their host clusters. 
Thus, most enrichment events occur shortly after the double neutron stars are born. 
This requires star-forming gas that can absorb the r-process ejecta to be present in the globular cluster 30–50 Myr after the initial burst of star formation. 
If scenarios for redistributing gas in GCs cannot act on these timescales, the number of neutron-star merger enrichment candidates drops severely, and it is likely that another mechanism, such as r-process enrichment from collapsars, is at play. 









§ INTRODUCTION


Astrophysical mechanisms for synthesizing the heaviest elements in the universe are poorly understood, yet essential in explaining the nucleosynthetic abundances observable today. 
Roughly half the elements heavier than iron are formed through the rapid capture of neutrons in a dense, neutron-rich environment, known as r-process nucleosynthesis <cit.>. 
In these environments, the rate of neutron capture overcomes the rate of β-decay of radioactive nuclei, which converts the heavy nuclei into more stable isotopes with higher atomic numbers. 



Once the prevalent paradigm, regular core-collapse supernovae (CCSNe) are now strongly disfavored both theoretically <cit.> and observationally <cit.> as a production site for heavy r-process elements.
The neutron-rich ejecta from neutron-star mergers (NSMs) is theorized to fill these open gaps in the periodic table, polluting the universe with the heaviest naturally occurring elements <cit.>. Over the last two decades, numerical simulations of the merger and post-merger phase have revealed several mechanisms by which neutron-rich material is ejected from these systems, including dynamical ejecta of tidal and shock-heated nature <cit.>, neutrino-driven and magnetically driven winds from a (meta-)stable remnant <cit.>, and outflows from a post-merger accretion disk <cit.>. The fact that NSMs can indeed synthesize r-process elements was corroborated in 2017 August; the discovery of GW170817 and observations of the subsequent kilonova explosion <cit.> provided unequivocal evidence that NSMs are a site for heavy r-process nucleosynthesis <cit.>. 
If NSMs are the main channel for heavy r-process enrichment <cit.>, r-process abundances in various astrophysical environments can thus probe the veiled physical mechanisms of binary stellar evolution, the rates of compact object mergers, and dynamical processes. 

Globular clusters (GCs) are one environment where r-process-enhanced, metal-poor stars have been observed, with a significant number of GCs exhibiting star-to-star dispersion of r-process species such as Eu and La <cit.>. 
The homogeneity of iron-group abundances in these clusters, as well as the lack of correlation between the internal spread of r-process elements and light element dispersion, makes these observations difficult to explain in the standard CCSNe r-process scenario; such stellar explosions would likely introduce similarly high levels of Fe species into r-process-enhanced stars <cit.>. 
Furthermore, measurements of [Pb/Eu] ratios indicate that these clusters have negligible dispersion of s-process elements <cit.>. 
Therefore, the astrophysical mechanism for introducing r-process elements must be inefficient at s-process production, or the dispersion of s-process elements is washed out from other, more frequent events. 

A number of observational campaigns have investigated the ubiquity of r-process enhancement in GCs. 
From a sample of 11 low-metallicity GCs, <cit.> found that four GCs showed clear signs of large r-process dispersion, five GCs showed no dispersion, and two GCs were more ambiguous with smaller levels of dispersion. 
Out of the five GCs that had no clear r-process dispersion, two showed bimodal chemical abundance ratios, possibly suggesting formation through the merger of two separate GCs that may have washed out initial star-to-star dispersion <cit.>. 
Later work also found one of the enriched clusters from this sample to show no significant r-process dispersion <cit.>. 
More recent studies have expanded upon this sample and added more candidate GCs with and without evidence of r-process dispersion <cit.>. 
Together, these samples show that heavy-element dispersion is a common feature of GCs, although further work is needed to robustly quantify the fraction of GCs that are enriched. 
Since heavy-element dispersion is apparent in a significant number of GCs but not ubiquitous across the full population, the mechanism for introducing r-process dispersion into these environments must occur in many, but not all, GCs. 

These observations have a few immediate consequences. 
First, they indicate some kind of extended or secondary star formation episode early in the history of the GC, leading to stellar populations with distinct chemical abundances due to astrophysical processes that transpired over this time.
Second, there must be an interstellar medium in the young GC with a high enough density to reduce the energy and momentum from the relativistic r-process ejecta and prevent it from leaving the GC <cit.>. 
Third, to provide this dense interstellar medium, there is likely a mechanism that replenishes the natal GC environment with gas for forming a second generation of stars <cit.>, as early-time SNe and feedback from of massive-star winds will efficiently remove the natal gas embedded in the GC. 
Finally, since GCs are believed to form the majority of their stars within the first ∼ 10 Myr after formation <cit.>, the enrichment mechanism had to proceed on a relatively rapid timescale. 

If NSMs are assumed to be responsible for the r-process enrichment in the second generation of stars formed in GCs, there are multiple stringent constraints on the properties of double neutron-star (DNS) systems at birth, most notably inspiral times and post-SN systemic velocities. 
Given the shallow gravitational potential expected of GC progenitors of ∼ 10 – 100 km s^-1 and small physical sizes of only a few parsecs, if the newly formed DNSs attain appreciable post-SN systemic velocities, they are typically ejected from the cluster and must rapidly merge before fully evacuating the cluster environment. 
If newly formed DNSs receive small enough barycentric kicks to remain in the cluster environment, they still must merge before the cluster is void of the gas that will form the second generation of (enriched) stars. 

Dynamical interactions may also play a role in expediting NSMs in GCs <cit.>, though these will need to have a significant impact early in the cluster lifetime. 
Though the NSM scenario for r-process enhancement has been explored in other environments such as ultra-faint dwarf galaxies <cit.>, GCs provide a unique and complementary probe for investigating NSMs as the primary site for r-process enrichment due to their vastly different physical sizes and masses, shorter star formation timescales, and the possible role of dynamical encounters. 

In this paper, we examine multiple scenarios for rapid r-process enrichment in GCs using both of semi-analytic arguments and numerical modeling.
In Section <ref>, we investigate NSMs produced by dynamical interactions, finding that this channel does not contribute to NSMs at early times. 
Section <ref> examines the impact that compact binaries formed from isolated primordial stellar pairs could have on enrichment. 
We show that primordial binaries can only lead to an appreciable number of r-process enhanced GCs if systems proceed through a phase of mass transfer (MT) involving a naked He-star donor. 

In Section <ref>, we discuss the effect that MT and SNe prescriptions have on the enrichment candidates in our models, review the implication of our results on the properties of GC progenitors and star formation timescales, compare our results to other environments that have observed r-process enhanced stars, and suggest alternative scenarios that may lead to r-process enhancement in GCs. 
Finally, we highlight our main conclusions in Section <ref>. 
In all population modeling, we assume a Kroupa <cit.> initial mass function (IMF) ranging from 0.08 M_⊙≤ M ≤ 150 M_⊙, unless otherwise specified. 





§ DYNAMICAL ASSEMBLY OF DNS SYSTEMS


GCs host complex dynamics between stars and compact objects. 
Through dynamical friction, the most massive objects in a cluster sink to its core, where they readily take part in strong gravitational interactions. 
These interactions can induce rapid mergers of compact objects <cit.>. 
However, via energy equipartition, as the most massive objects migrate to the cluster core, lighter objects move further away from the cluster center. 
As a consequence, black holes (BHs) shape the dynamical evolution of the lower-mass (M ≲ 1 M_⊙) main-sequence stars and, because these lighter stars make up the bulk of the total cluster mass, the cluster as a whole <cit.>. 
At early times when a large BH population is present, BHs therefore inhibit the dynamical segregation of the NSs that could potentially merge and enrich the cluster with r-process elements. 
Nevertheless, it has been argued that dense stellar environments such as nuclear star clusters in the early universe can yield an appreciable rate of dynamically assembled NSMs and possibly explain the r-process enhancement of metal-poor stars <cit.>. 

Over time, BHs will be ejected from the cluster through strong encounters with other BHs in the core <cit.>, allowing lighter objects in the cluster to follow suit and migrate toward the core through a similar segregation process, ultimately leading to cluster core-collapse <cit.>. 
The depletion of BHs in the core also permits the segregation of NSs to the cluster core which leads to an increased formation rate of millisecond pulsars, DNSs, and possibly NSMs <cit.>. 

Mass segregation timescales can be approximated by the time necessary for a GC to settle into equilibrium. 
The half-mass relaxation time is given by

    t_ relax = 0.138 M_ c^1/2R_ h^3/2/⟨ m ⟩ G^1/2lnΛ,

where M_ c is the total cluster mass, ⟨ m ⟩ is the average stellar mass, R_ h is the half-mass radius, and lnΛ is the Coulomb logarithm, where Λ≃ 0.4  N for GCs with N being the total number of stars <cit.>. 
For simple dynamical friction in a two-component model, the mass segregation timescale is

    t_ ms^i ∼⟨ m ⟩/m_i t_ relax,

where m_i is the mass of the segregating population <cit.>. 
For a typical cluster with M_ c≃ 4 × 10^5 M_⊙, ⟨ m ⟩≃ 1 M_⊙, R_ h≃ 1 pc, and N ≃ 8 × 10^5, and the half-mass relaxation time is t_ relax≃ 100 Myr <cit.>. 
Therefore, 30 M_⊙ BHs segregate in t_ ms^ BH≃ 3 × 10^6 yr, whereas 1 M_⊙ NSs segregate in t_ ms^ NS≃ 10^8 yr. 

To investigate this channel, we model GCs using the Hénon-style Monte Carlo code  <cit.>. 
This code monitors and evolves the global properties of the GC via two-body relaxation, while accounting for binary stellar evolution <cit.> and small-N gravitational encounters using the  package <cit.>, which now includes post-Newtonion effects in N-body integrations <cit.>. 
We record all strong encounters between objects, tracking the times at which new compact binary pairs are formed and calculating the gravitational-wave (GW) inspiral times of binaries synthesized due to these dynamical encounters. 

We simulate two clusters: one with typical GC properties and one with a number of liberal assumptions in an attempt to eliminate BHs from the environment so NSs will segregate to the cluster cores and interact as rapidly as possible. 
In this liberal model, we assume BHs have no mass fallback following the SN at their formation and therefore receive the full natal kicks typical of NSs, and truncate the IMF at 20 M_⊙ (which corresponds to a remnant mass of ∼ 5-8 M_⊙ at low metallicities; ). 
This acts to reduce the number of BHs that are created from the initial stellar population, and efficiently eject those that do happen to form. 
Furthermore, we reduce all NS natal kicks, drawing their magnitude from a Maxwellian with a dispersion of 20 km s^-1 instead of a dispersion of 265 km s^-1 <cit.>. 
Figure <ref> shows the time that compact binaries are dynamically formed in these models compared to their inspiral times. 
The sum of these two quantities approximates the time in the history of the cluster at which these compact binaries would merge, and in the case of systems with a NS component, possibly enrich the cluster with r-process material. 





As expected, in the standard model binary BHs are the dominant dynamically induced systems and inhibit the formation of other types of compact binaries for most of the GC lifetime. 
Only a small number of BHs merge within 100 Myr of the formation of the GC, and none merge within 10 Myr. 
Primordial mass segregation in star clusters may help expedite dynamically induced BH mergers <cit.>, though the interaction rate of other types of compact objects will be stifled until the higher-mass BH population is processed and ejected. 

In our liberal model that artificially removes most BHs and amplifies the number of retained NSs, we still do not find NSMs until a few Gyr after cluster formation — far too late to enrich a second generation of stars with r-process material. 
Though BHs are efficiently removed from the cluster in this model, the 1–2 M_⊙ NSs are still less massive than the surviving main-sequence stars. 
Therefore, mass segregation is still inefficient for ∼ 3 Gyr until these stars finish their main-sequence evolution. 
Though dynamically formed neutron star-black hole (NSBH) systems may also play a role in enriching our liberal GC model, the NSBH rate is lower than the DNS rate and they share similar inspiral time distributions; hence, NSBH enrichment should be subdominant to DNS enrichment. 
We therefore conclude that dynamical interactions in GCs do not contribute to early-time NSMs, and it is unlikely that dynamics plays a significant role enhancing second-generation GC stars with r-process material. 








§ MERGERS FROM PRIMORDIAL CLUSTER BINARIES


We next investigate NSMs from primordial binaries in GCs following their isolated binary stellar evolution.[In this context, primordial binaries refer to stars that are born as binary pairs in the GC and have yet to segregate and undergo any dynamical interactions.]  
To pollute the cluster with r-process elements shortly after formation, DNSs must acquire short inspiral times through tight orbital separations and/or significant orbital eccentricities at birth. 
The standard picture by which compact binaries achieve hardened orbits during isolated binary evolution is via (at least) one common envelope (CE) phase, where one of the stars experiences Roche-lobe overflow (RLO) and unstable MT, enveloping both bodies in the outer layers of the donor star and exerting a drag on the orbiting binary system. 

The common envelope causes the binary to spiral inwards and harden its orbit <cit.>. 
For compact binaries that merge within a Hubble time, a CE phase will typically need to occur after one of the stars in the binary has already evolved into a compact object, so that the orbit can harden sufficiently without the bodies coming into contact and merging. 
Though the system circularizes following this evolutionary stage, eccentricity imparted into the system via the subsequent SN that forms the second NS will further reduce the GW inspiral time. 




 §.§ Enrichment Efficiency


In addition to altering the orbital properties of the compact binary, the second SN can significantly kick the system through a combination of asymmetries in the explosion mechanism (the natal kick) and from symmetric mass loss of the exploding star <cit.>. 
The post-SN center-of-mass velocity is referred to as the systemic velocity. 
Given a final mass M_ f = m_1+m_2 and initial mass M_ i = m_1+m_2+Δ M_ SN where m_1, m_2, and Δ M_ SN are the first-born NS mass, the second-born NS mass, and mass lost during the SN that creates the second NS, we can write the magnitude of the post-SN systemic velocity as

    v_ sys^2 = 1/M_ i^2[
        (m_2 v_ k)^2 
        + m_1 Δ M_ SN v_ r/M_ f^2(m_1 Δ M_ SN v_ r - 2 m_2 M_ f v_ k ∥)
        ],

where v_ k is the magnitude of the natal kick, v_ r is the relative orbital velocity between the two objects prior to the SN, and v_ k ∥ is the component of the natal kick aligned with the instantaneous orbital velocity of the exploding star <cit.>. 

We consider DNS systems viable r-process polluters if the system merges within some enrichment radius of the cluster, R_ enrich, and within some span of time when new stars are still being born after the initial burst of star formation, Δτ_ SF. 
This can happen in two ways: 
(i) the post-SN systemic velocity is greater than the cluster escape speed, but the DNS merges before leaving the cluster environment; (ii) the post-SN systemic velocity is less than the cluster escape speed, and the DNS delay time is less than Δτ_ SF.
To derive an approximate escape velocity, we assume the mass distribution in our fiducial clusters follow a Plummer profile <cit.>: 

    ρ_p(r) = (3 M_ c/4 π R_p^3) ( 1 + r^2/R_p^2)^-5/2,

where M_ c is the mass of the cluster and R_p is the Plummer radius. 
Given a Plummer radius, we assume DNS systems form at the half-mass-radius, which is r_h ≈ 1.3 R_p for a Plummer sphere. 
The escape velocity is thus 

    v_ esc(r_h) = √(2 |Φ_p (r_h)|)≈ 72 (M_ c/10^6 M_⊙)^1/2(R_ p/1  pc)^-1/2 km s^-1,

where Φ_p is the gravitational potential of the Plummer model. 
Though this is a simplistic description for the true potential of a young GC, we find our results to be robust to changes in the form of the potential since the post-SN systemic velocity of most enrichment candidates greatly exceeds the GC escape velocity. 

Assuming negligible deceleration from the gravitational potential, if v_ sys > v_ esc, the time it takes for the DNS to go beyond the enrichment radius of the cluster is 

    τ_ eject≈R_ enrich/v_ sys.
 
We typically assume that the enrichment radius is R_ enrich = R_ vir≈ 1.7 R_p <cit.>, though later we will relax this assumption. 

Systems that remain bound to the cluster after being kicked and merge within Δτ_ SF are also viable r-process polluters. 
The delay time is defined as t_ delay = t_ DNS + t_ insp, where t_ DNS is the time from zero-age main sequence (ZAMS) of the DNS progenitors to DNS formation and t_ insp is the DNS inspiral time due to GW emission <cit.>. 
We therefore categorize DNS systems that remain bound as viable polluters if t_ delay < Δτ_ SF. 
Thus, the DNS enrichment efficiency is 

    ϵ_ DNS = 1/N_ DNS∑_i = 1^N_ DNSΘ(Δτ_ SF^i - t_ delay^i) ×𝒯,

where N_ DNS is the total number of DNS systems, Θ is the Heaviside step function, and 𝒯 = Θ(τ_ eject^i - t_ insp^i) if v_ sys > v_ esc, and 𝒯 = 1 if v_ sys≤ v_ esc. 





 §.§ NS Natal Kicks


Proper motions of isolated pulsars in the Milky Way indicate that many NSs receive large natal kicks at birth, on the order of a few hundred kilometers per second (, though see also ). 
As described above, post-SN systemic velocities of DNS systems are affected by natal kicks, and can lead to systemic velocities that are ≲ 50% larger than the pre-SN orbital velocity <cit.>. 

Since the escape speeds of GCs are typically a few tens of kilometers per second, DNSs in relatively tight pre-SN orbits that are not disrupted from standard CCSNe natal kicks will usually lead to post-SN systemic velocities that unbind the systems from their host cluster. 

Evidence for some NSs receiving lower natal kicks has been determined by examining the proper motions of both the isolated NS population <cit.> and Galactic DNSs <cit.>. 
These lower kick magnitudes of a few tens of kilometers per second are predicted for DNSs that explode due to electron capture in a strongly degenerate ONeMg core, known as electron-capture SNe <cit.>.
ECSNe are typically assumed to occur when a star has a He core mass of around 2 M_⊙ at the base of the AGB branch <cit.>, though the extent and placement of this range is debated and can be strongly affected by binary interactions <cit.>. 
This pathway will make it more likely for a DNS to remain bound to the cluster (as well as survive the SN), since lower natal kicks and smaller amounts of mass loss lead to smaller post-SN systemic velocities. 
However, the eccentricities imparted into these systems from the SN will likewise be smaller, increasing the typical GW inspiral time. 

In our population modeling, we apply a bimodal distribution for natal kick magnitudes. 
Standard CCSNe have natal kicks drawn from a Maxwellian distribution with scale parameter σ_ high = 265 km s^-1 <cit.>, whereas for stars with He core masses at the base of the AGB branch in the range 1.4 M_⊙≤ m_ core≤ 2.5 M_⊙, which are predicted to undergo ECSNe <cit.>, we draw natal kicks from a Maxwellian distribution with σ_ low = 20 km s^-1. 
Another possible mechanism for stifling natal kicks of DNSs is by significantly stripping the atmosphere of the progenitor star via binary interactions prior to the SN, a scenario known as ultra-stripped SNe <cit.>, which we also take into account in one of our population models detailed below. 
We also apply natal kicks from the σ_ low distribution for systems that form through accretion-induced collapse of an ONeMg white dwarf <cit.> or a merger-induced collapse <cit.>, though these channels (especially the latter) do not typically lead to DNS formation. 






 §.§ Case BB MT


One aspect of binary evolution that affects DNS populations is the onset, stability, and outcome of MT that results from a post-He main-sequence star overflowing its Roche lobe, known as Case BB MT <cit.>. 
Case BB MT is believed to occur during the He-burning analog of the Hertzprung gap (HG), when core He burning has ceased and shell burning causes the star to expand. 

Originally, it was assumed that the RLO of low-mass He stars that entered the HG was unstable and led to a successful CE phase <cit.>, causing the spiral-in of the already-formed NS and the core of the He-star donor. 
However, the ability of binary systems to proceed through this phase is highly uncertain. 
For example, dynamical instability during RLO may lead to merger for donor stars without a core-envelope structure or a clear entropy jump at the core-envelope transition <cit.>. 
<cit.> demonstrated with detailed stellar evolution simulations that most unstable MT during this evolutionary phase leads to delayed dynamical instability and merger, and for DNS progenitors the MT typically proceeds stably during He-shell burning. 
More recent work by <cit.> has confirmed this, and found that the stability may be even more prevalent in DNS progenitors. 
This phase of binary evolution has important implications on fast-merging DNSs: Case BB MT is predicted to lead to hardened DNS systems with orbital periods of ≲ 10^-2 days at formation in the most extreme cases, corresponding to GW inspiral times of 𝒪(10^3) yr, <cit.>. 






 §.§ Population Models


To determine ϵ_ DNS, we simulate multiple populations of merging DNSs using the  population synthesis code <cit.>.[https://cosmic-popsynth.github.io/cosmic-popsynth.github.io] 
 is a modified version of  <cit.>, which relies on polynomial fitting formulae for single-star evolution <cit.> and includes physical prescriptions for binary evolutionary processes such as tidal evolution, MT, CEs, and GW decay <cit.>. 
 is updated to include state-of-the-art prescriptions for mass loss in O and B stars <cit.>, metallicity dependence in the evolution of Wolf-Rayet stars <cit.>, new prescriptions for fallback and post-SN remnant masses <cit.>, variable prescriptions for the CE λ parameter <cit.>, as well as prescriptions for ECSNe <cit.>, USSNe <cit.>, and (pulsational) pair instability SNe <cit.>. 
A modification to  pertinent to this study is a correction to an inconsistency regarding the masses and orbital separations used when an SN occurs directly after a CE phase; more details can be found in Appendix <ref>. 
In addition,  determines when particular populations of compact binaries have been adequately sampled by repeatedly checking for convergence in the distributions of various binary properties <cit.>. 

The four models we explore in this study cover a range of uncertainties in binary evolution by varying certain aspects of Case BB MT, survival through CEs, and natal kick prescriptions. 
In all models, initial binary properties are determined according to the prescriptions in <cit.>, and we assume that all stars form in a single burst of star formation with a metallicity of Z_⊙/20.[Though DNS merger rates are strongly impacted by metallicity, properties of merging DNSs vary only slightly with metallicity <cit.>, and thus our choice of metallicity should have little effect on our results.]
In addition to providing astrophysically motivated distributions of post-SN orbital properties and systemic velocities,  tracks the total sampled mass necessary to generate our population of DNSs. 
The converged simulations resulted in 𝒪(10^4) DNS systems per population, corresponding to a total sampled mass of ≈ 10^9 M_⊙. 
Details of each population model are as follows. 












    

    


    
    
  * Standard model, which allows for evolved He stars to proceed through a successful CE.
    NS masses are calculated using the Delayed prescription from <cit.>.[Models run using the Rapid prescription from <cit.> only led to percent-level differences in our results.]  
    For naked He stars, we adopt maximum mass ratios for stable MT from <cit.>: q_ crit = 1.7 for main-sequence naked He stars (denoted as stellar type k_⋆=7 in ) and q_ crit = 3.5 for HG or giant branch naked He stars (k_⋆=8, 9). 
    Other stellar types use the default  q_ crit values. 
    
    
    
    

    
  * Same as Model A, except that unstable MT from donor stars without a well-developed core-envelope structure always lead to a merger. 
    This is also assumed to be the case for donors on the HG, since they lack a clear entropy jump at the core-envelope boundary <cit.>. 
    Analogous to the pessimistic CE model in <cit.>, this drastically limits the number of short-period systems since it eliminates potential DNS progenitors that would have otherwise gone through a Case BB CE phase. 

    
  * Same as Model A, except orbital properties and mass loss following MT with a He-star donor are calculated according to the prescriptions in <cit.>, who found that this phase of MT to typically proceeds stably and does not initiate a CE phase, and provided fitting formulae that map pre-RLO masses and separations to their values following stable Case BB RLO (i.e. the values immediately prior to the second SN).[Though further phases of stellar evolution occur before the SN, they proceed on rapid timescales and any stable or unstable MT will not strongly affect the orbital properties of the pre-SN system. ]
    Though these fitting formulae are derived from systems at solar metallicity with pre-CE He-star masses between 2.5 ≤ M_ He/M_⊙≤ 3.5 and companion masses of m_ NS = 1.35 M_⊙, we use these fits for all He-star–NS systems that undergo Case BB MT. 
    Extrapolating to larger pre-RLO He-star masses may lead to larger uncertainties in post-RLO parameters, and exact results could potentially differ at the lower metallicities considered here. 

    
  * Same as Model C, except all systems that undergo Case BB RLO are assumed to become ultra-stripped and receive the same natal kicks as ECSNe. 


The time t_ DNS between ZAMS and DNS formation for our population ranges from ≈ 11 to 68 Myr, with a median value of t_ DNS = 41 Myr. 
When calculating t_ delay for an individual system, we use the t_ DNS found from our modeling of that system and integrate its post-SN semi-major axis and eccentricity according to <cit.> to determine t_ insp. 



 §.§ Enrichment probability from primordial DNSs





Figure <ref> shows the orbital properties immediately following the second SN for all DNSs that merge within a Hubble time in our population models, as well as the corresponding post-SN systemic velocities and inspiral times. 
We assume that the first generation of stars is born as a single burst at t=0 and that the second generation of (enriched) stars must be born within Δτ_ SF = 100 Myr, though we investigate variations in the latter assumption in Section <ref>. 
Black lines separate DNS systems that are viable enrichment candidates (lower left) from those that either merge outside the cluster or after the star formation has ceased (upper right) for three assumed cluster sizes. 
Therefore, the enrichment fraction ϵ_ DNS is given by the fraction of systems that lie to the lower left of these lines. 
For the remainder of the population, systems that undergo a Case BB MT phase are colored in purple whereas those that proceed through only a single CE after the first NS formed (typically during the hydrogen giant branch of the secondary) are in gray. 
Systems that undergo Case BB MT provide the largest contribution to ϵ_ DNS, otherwise viable enrichment candidates can only result from systems that are kicked into an extremely eccentric orbit at birth (e.g., the small number of red points in Model B). 
However, as shown in Eq. (<ref>), the post-SN systemic velocity scales with the orbital velocity of the binary prior to SN. 
Therefore, tighter pre-SN binaries will have larger Blaauw kicks; this anti-correlation between systemic velocity and inspiral time for the Case BB systems in evident in Figure <ref>. 

To calculate the fraction of GCs we expect to exhibit r-process dispersion, we must first determine whether a single NSM can distribute enough r-process material to enrich the stars in a GC. 
The total Eu mass of a GC can be roughly estimated by M_Eu≈ X_Eu,⊙ 10^[Eu/H] M_gas, where X_Eu,⊙ is the solar Eu mass fraction and M_gas is the mass of intra-cluster gas in which the r-process was mixed into. For typical values, we have

    M_Eu≈ (4.2 × 10^-5) × 10^[Eu/H](M_gas/10^5 M_⊙) M_⊙,

where we have used the solar abundances from <cit.>. 
This amounts to M_Eu∼ 4 × 10^-6 M_⊙ assuming typical GC abundances from <cit.>. 
Assuming that an NSM can eject ≈ 0.05 M_⊙ of r-process material, consistent with GW170817 <cit.>, and assuming a solar abundance pattern starting at mass number A=69 <cit.>, this translates into a Eu mass of M_Eu≈ 5.2× 10^-5 M_⊙. 
We therefore conclude that a single NSM may be sufficient to enrich the second generation of stars in a typical GC.

Each model in Figure <ref> is synthesized using a sample of binaries far larger than the stellar mass of a young GC. 
To quantify the probability that an enrichment event will occur in a GC with a stellar mass at birth of M_ c, we must determine the rate at which an NSM enrichment event would occur in a GC with a given stellar mass budget. 

The sampling procedure used in  from <cit.> determines the probability that a system with a given primary mass and orbital period will have a companion.[The prescriptions from <cit.> are based on observations of binaries in the Galactic field. 
Though there is significant uncertainty in the initial binary fraction of high-mass stars in GCs <cit.>, observations of young stellar clusters still find that most massive stars are in binary pairs <cit.>, and we expect slight variations in binary fractions to have a percent-level impact on our results. ]
Though  only evolves binary systems, the total mass sampled includes the mass of both single and binary star systems. 


Therefore, the expected number of enriching NSMs in a GC is 

    Λ_enrich = ϵ_ DNS N_ DNSM_ c/M_ samp,

where M_ samp is the total mass sampled in our population models. 
Given the mean number Λ_enrich, the probability that a GC has at least one enriching NSM follows a Poisson process as

    P_enrich = 1 - exp(-Λ_enrich).

The expected enrichment probabilities for our four population models using multiple assumptions for R_ enrich and M_ c are in Table <ref>.[We do not consider NSBH mergers in our enrichment probability calculations. 
The upper limit on the NSBH merger rate density measured by LIGO–Virgo is already lower than the median DNS merger rate density <cit.>. 
The mass fraction of binary stars that results in NSBH mergers should therefore be lower than that for DNSs, and thus their enrichment contribution in GCs should be subdominant. ]











For a given population model, enrichment fractions generally increase with increasing cluster mass and cluster size. 
Increasing R_ vir from 1 to 10 pc increases the enrichment probability in most models by a factor of ≈ 2-3. 
For low to moderate values of P_ enrich, enrichment fractions are typically ≈ 2 times larger for M_ c = 10^6 M_⊙ compared to M_ c = 5 × 10^5 M_⊙: while the change in M_ c affects the GC escape velocity (and thus ϵ_ DNS), the primary contribution to Λ_enrich comes instead from the increased number of potential DNS progenitors, which scales linearly with M_ c.
However, as one pushes to extreme assumptions about the initial GC mass, this linear scaling breaks down since many more DNSs remain bound to the natal GC. 
To this end, we also examined the enrichment fractions recovered when assuming an extreme assumption about the initial mass: M_ c = 10^7 M_⊙. 
In this case, regardless of the cluster size, Models A, C, and D all have expected numbers of enriching NSMs much greater than unity (Λ_ enrich≫ 1), and therefore enrichment probabilities of P_ enrich ≈ 100%. 
However, for this assumed cluster mass, Model B still only reaches enrichment probabilities of P_ enrich ≈ 50%. 


The P_enrich results in Table <ref> assume that an NSM must occur within the virial radius of the young GC in order to enrich the still-forming second generation of stars. 
However, it may be possible to enrich the stellar populations from outside the cluster environment, though this will lead to a geometrical R^-2 reduction of r-process material available for pollution, assuming the ejecta is isotropically dispersed. 
We examine P_enrich as a function of enrichment radius in Figure <ref> for Model D. 
For a GC with an initial mass of 5 × 10^5 M_⊙, enrichment probabilities hit ≈ 50% at R_ enrich = 20 R_ vir (2 R_ vir) if R_ vir = 1 pc (10 pc). 
For most assumptions about initial GC properties, Model B does not exceed enrichment probabilities of even 20% unless R_ enrich ≳ 100 R_ vir. 
This indicates that the pessimistic CE model is incompatible with substantial enrichment fractions, regardless of the cluster initial conditions or assumed enrichment radius. 

















§ DISCUSSION



In the following sections, we highlight the implications of GC enrichment fractions on the physical assumptions of binary stellar evolution, discuss how initial properties and star formation timescales of GCs can be constrained by such observations, touch on other possible mechanisms for r-process enrichment in GCs, and comment on implications for other environments that have been observed to be r-process enhanced. 




 §.§ r-process Enrichment from Primordial NSMs in GCs


Model A, which allows for successful Case BB CEs, has enrichment probabilities ranging from ∼ 30% to 90%. 
Though the majority of DNS systems have post-SN systemic velocities greater than the escape velocity of the fiducial GCs, the tight orbital separations achieved from successful Case BB CEs still allow for an appreciable number of DNS mergers to occur within the virial radii of the GC. 
This model is consistent with most, if not all, GCs with virial radii ≳ 3 pc and/or initial masses ≳ 10^6 M_⊙ showing star-to-star r-process dispersion if star formation does indeed persist for ≈ 100 Myr. 
Since star-to-star r-process dispersion is not a ubiquitous phenomenon in GCs, this model may overpredict NSM enrichment candidates for most assumptions about initial GC properties. 

In Model B, which implements a pessimistic CE scenario not allowing for CEs involving an evolved He-star donor <cit.>, we find few enrichment candidates. 
Regardless of initial GC size and compactness, P_enrich is consistently ≲ 0.06; the only enrichment candidates are those that proceed through a single CE and are kicked into highly eccentric orbits from the second SN. 
This shows that without some form of Case BB MT, NSMs are unable to explain r-process enhancement in more than a few percent of GCs. 

Following <cit.>, in Models C and D we assume that Case BB MT typically proceeds stably (Case BB MT does not lead to a CE phase). 
However, as the mass of the donor star is always greater than that of the first-born NS when considering DNS progenitors, stable MT will still cause the binary orbit to shrink. 
Since the orbital hardening is less drastic than if a CE occurred, we find lower Λ_ enrich compared to Model A by a factor of ∼ 2 for Model C and ∼ 3 for Model D. 
The reason for lower enrichment levels in Model D compared to Model C is the stifled natal kicks for systems that go through Case BB MT. 
This leads to less scatter in the post-SN orbital properties, a tighter trend in v_ sys - t_ insp space, and a slightly lower enrichment probability. 
Both of these models are consistent with many, but not all, GCs having an NSM within its first 100 Myr. 
This demonstrates that in the NSM scenario with extended star formation, GC enrichment can be explained without invoking a Case BB CE phase. 



















 §.§ Sensitivity to CE Prescriptions for HG He stars


As shown in Section <ref>, the vast majority of viable enrichment candidates proceed through a phase of Case BB MT. 
The degree of orbital hardening from this stage of binary evolution strongly affects the post-SN systemic velocities and inspiral times, and thereby has an important impact on the number of NSM enrichment candidates in a GC. 

Since pre-SN orbits of systems that went through Case BB MT are slightly wider for Models C and D than for Model A, mass loss during the SNe in Model A typically leads to larger changes in barycentric velocity. 
In Model A the median post-SN systemic velocity for systems that undergo Case BB MT is ∼ 600 km s^-1, whereas Models C and D have median systemic velocities of ∼ 400 km s^-1 and ∼ 200 km s^-1, respectively. 
Though these systems are less likely to be ejected, as seen in Figure <ref>, the majority still exceed the escape velocity of a typical GC. 

Case BB systems in Models C and D have slightly longer inspiral times, which leads to slightly lower enrichment fractions than when Case BB MT proceeds unstably. 
Systems evolving through Case BB CE (Model A) generally lead to DNSs with smaller separations and hence shorter inspiral times than those evolving through stable Case BB MT (Models C, D; in these models 41% and 27% of systems have t_ insp < 10^5 yr, respectively, whereas this fraction reaches 73% in Model A). 
This is particularly evident in Model D, where Case BB donors are assumed to be ultra-stripped and have stifled natal kicks. 

In the case of Model B, we impose the pessimistic CE scenario, causing all stars that undergo unstable MT as an evolved He-star to merge. 
This removes Case BB systems and thereby depletes the population of DNSs with short inspiral times, eliminating most viable enrichment candidates. 
With only large-separation DNS progenitors surviving, in this model only systems that get kicked into eccentricities close to unity merge quickly enough to potentially enrich the GC at early times; only ≈ 0.2% of DNS systems merging within a Hubble time in this model have inspiral times less than 100 Myr and systemic velocities less than 100 km s^-1. 

Overall, significant hardening of the pre-SN binary has a counteractive effect on the enrichment efficiency. 
Though hardened pre-SN orbits typically lead to shorter inspiral times, they also amplify post-SN systemic velocities; even if the natal kicks are small, high post-SN systemic velocities will be achieved through Blauuw kicks <cit.>. 
This causes hardened systems to be more easily ejected from the GC. 
For example, DNS progenitors with inspiral times of ≲ 100 Myr will typically require an envelope mass prior to the second SN of ≲ 0.1 M_⊙ in order to remain bound to the GC. 
However, the hardest DNS progenitors, with pre-SN orbital separations of ≲ 0.2 R_⊙, will typically be ejected from the GC regardless of the amount the secondary was stripped (see Appendix <ref>). 
Larger pre-SN orbital separations will reduce the effect of the Blaauw kick, but will also increase inspiral times. 
In the context of GCs, enrichment candidates are typically in the first of these regimes: ejected from the GC but residing in the tail of the inspiral time distribution such that they merge before leaving the GC entirely. 
We explore implications of Blauuw kicks on rapid-inspiral DNS progenitors in Appendix <ref>. 













 §.§ Initial Cluster Properties



The GCs with observed r-process dispersion have varying properties; both core-collapsed and non-core-collapsed clusters show dispersion, and present-day masses of these clusters range from ∼ 2 × 10^5 to 5 × 10^5 M_⊙ <cit.>. 
Since most of the Milky Way's GCs formed ∼ 10 Gyr ago, initial GC properties are highly uncertain. 
Observations of super star clusters, which are believed to be the progenitors of GCs, help to elucidate the initial conditions of these old stellar systems. 
These observations indicate that the stellar masses of GCs found in the Milky Way today are a factor of a few times less massive than they were at formation <cit.>. 





More massive GCs at formation lead to more enrichment candidate NSMs due to deeper gravitational potentials (higher escape velocities) and, more importantly, more stars sampled from the IMF that can potentially form DNS binaries. 
In our analysis, we investigated two representative natal GC masses of 5 × 10^5 and 10^6 M_⊙, working under the assumption that present-day GCs were a few times more massive at birth <cit.>. 
However, with knowledge of the true population model representing DNSs in young GCs, 

the inference can be reversed to place constraints on initial conditions of GCs necessary to match the number of enriched GCs.
Figure <ref> demonstrates this by highlighting the regions of initial cluster mass–initial cluster size space that lead to the enrichment fractions between 25% and 75%, broadly consistent with the fraction of GCs with star-to-star r-process dispersion in the limited observational samples to date <cit.>. 

Assuming the enrichment event must occur within the virial radius of the cluster and Case BB MT proceeds stably (Models C and D), this maps to typical GCs with initial stellar masses of ∼ 8 × 10^5 (5 × 10^5) M_⊙ given a physical size of ∼ 3 pc (10 pc) at formation. 
Though this is a few times more massive than the present-day masses of some of the GCs with observed r-process dispersion, simulations of massive star cluster formation do find stellar masses of up to ≈ 1–2 × 10^6 M_⊙ (, see also ). 
However, constraints on initial cluster properties should be taken with caution, since theoretical uncertainties in the evolution of DNS progenitors as well as the limited number of observational results concerning the r-process elements in GCs currently prohibit any definitive statements. 
If future observations of DNSs in both the Milky Way and via GWs can further constrain uncertain aspects of high-mass binary stellar evolution, the enrichment fraction observed in GCs could help probe these uncertain properties of GC progenitors. 




 §.§ Extended Star Formation in Young GCs








The timespan over which stars form from the natal gas in GCs is uncertain, though typically believed to be ≲ 10 Myr <cit.>. 
Observations of young GCs in the local universe find that they are already gas-free after ∼ 2-4 Myr <cit.>, and have no evidence for extended star-formation episodes <cit.>. 
The primary drivers for removing natal gas from young GCs are believed to be radiation pressure, photoionization, ram pressure, and feedback from the SNe of high-mass stars, which proceed through the entirety of their stellar evolution in only a few Myr. 
Despite this, mechanisms have been proposed for redistributing gas within young GCs, which needs to be sufficiently dense to halt the relativistic r-process ejecta from NSMs and to form new, r-processed-enhanced stars. 
For example, <cit.> suggest that this may be accomplished through the conglomeration of AGB ejecta from 6 to 9 M_⊙ stars, which can form a high-density, compact gaseous region in the center of a young GC. 
This gas would begin to pollute the cluster ∼ 100 Myr after the initial burst of star formation, past the time that gas chemically polluted by SNe of high-mass stars has been removed. 

We take a conservative star formation timescale of 100 Myr for most of the values presented in this paper, which is corroborated by certain observations of young massive star clusters <cit.>. 
Though the freefall timescale of a 10^6 M_⊙ molecular cloud about 10 pc in radius is only ≈ 1 Myr, star clusters in general do not form from a single isolated cloud but rather through the hierarchical merging of smaller clusters. 
If clusters form hierarchically in turbulent molecular clouds with freefall times of tens of millions of years, it may be possible for new star formation to occur ≈ 100 Myr into the cluster evolution as gas clouds continue to stream down the potentials of the newly formed clusters. 

Our results are largely insensitive to increases in Δτ_SF; as seen in Figure <ref>, few systems fall directly above the horizontal black line at t_ delay = 100 Myr, so increasing Δτ_ SF has little impact on ϵ_ DNS. 
For example, assuming a progenitor GC mass of 10^6 M_⊙ and virial radius of 3 pc, increasing Δτ_ SF from 100 Myr to 1 Gyr increases the enrichment probabilities by only 0.5%, 2.3%, and 3.2% in Models A, C, and D, respectively. 
Therefore, cluster formation scenarios that invoke more extended star formation episodes will not lead a to significant amplification of enrichment candidates. 

In contrast, results can change drastically if we decrease Δτ_ SF. 
In Figure <ref>, we show the enrichment fraction as a function of Δτ_ SF for various assumptions about the initial conditions of natal GCs. 
Regardless of GC properties, P_enrich plateaus at ≈  40 Myr, which is the median formation time of DNS systems in our models. 
For example, a value of Δτ_ SF = 30 Myr produces only ≈ 4–7% of the enrichment candidates as Δτ_ SF = 50 Myr across our population model. 
This is because the evolutionary timescales of DNS progenitors typically exceed Δτ_ SF for Δτ_ SF≲ 30 Myr. 

The steep buildup of P_enrich at 30-50 Myr has important implications for when the redistribution of gas in the GC must have occurred by. 
Since most DNS mergers that contribute to enrichment are those that are efficiently ejected from the GC, NSMs that are viable enrichment candidates must merge shortly after DNS formation. 
Therefore, for NSMs to be a viable scenario for r-process enhancement, star formation in GCs must be ongoing or rejuvenated at the time that DNSs form, ≈ 40 Myr after the initial burst of star formation. 
This criterion may be hard to accomplish using, for example, lower-mass AGB stars as the mechanism for re-polluting the GC with gas ≈ 100 Myr after formation, and may indicate that another mechanism for r-process enhancement in young GCs is at play. 








 §.§ Alternative Scenarios for Enrichment


As demonstrated in the above sections, the NSM scenario for r-process enrichment can explain moderate to high fractions of GCs with enrichment if star formation for forming the second generation of stars is ongoing at the time of DNS formation, ≈ 30–50 Myr after the initial bout of star formation. 
Though this may be possible, either through extended star formation or
a redistribution of gas from stellar processes after the first wave
of massive stars heat and eject the natal gas, other scenarios that can contaminate the cluster with r-process material at early times may also be viable. 

Various other astrophysical scenarios have been proposed for contributing to the observed abundances of r-process elements <cit.>.
One such scenario is via the formation of a dense accretion disk during the collapse of a massive, rapidly rotating star — the collapsar scenario <cit.>. This scenario has also been suggested to be the cause of long γ-ray bursts (GRBs; e.g., ). 



Collapsars arguably better satisfy existing constraints on r-process enrichment and may overcome problems or alleviate tensions that enrichment scenarios solely based on NSMs face (see  for a discussion of some of the issues). 
This relates to issues of prompt enrichment at low metallicities (e.g., ), and also to the high-metallicity stars of the Milky Way disk <cit.>. 
Though potential concerns have been raised regarding the co-production of iron in the context of extremely metal-poor stars <cit.>, this may not be applicable to r-process enrichment in GCs; given the total mass and metallicity of GCs that show internal r-process dispersion, many SNe likely occurred in young GCs.




Unlike enrichment from NSMs, the collapsar scenario is not susceptible to ejection from the GC due to large post-SN systemic velocities. 
Furthermore, as collapsars result from massive star evolution, they could occur within a few million years from the initial burst of star formation, when natal gas is still present in the cluster environment and the first wave of stars are still forming. 
However, for this scenario to be viable, the rate density of collapsars per unit of stellar mass must be in line with the total stellar mass available to the young GC. 

By assuming collapsars are the cause of long GRBs, the number of collapsars per unit stellar mass formed can be empirically estimated from the observed local long GRB rate. 
Assuming the rate of collapsars (events per unit comoving volume per unit time) tracks the star formation rate ρ_SF (stellar mass per unit comoving volume per unit time) with negligible delay (∼ few Myr), R_coll(t) ∝ρ_SF(t), where t denotes cosmic time. Thus, one can write the total number of collapsars per unit stellar mass formed as

    n_coll = ∫ R_coll(t) dt/∫ρ_SF(t) dt = R_coll(z=0)/ρ_SF(z=0).

The collapsar rate is given by R_coll(z=0)= R_LGRB(z=0)/f_b, where R_ LGRB(z = 0) ≈ 1.3_-0.7^+0.6 Gpc^-3 yr^-1 is the rate of local (z=0) long GRBs <cit.> and f_b≈ 5× 10^-3 is the GRB beaming fraction <cit.>. 
This GRB rate does not include the separate class of low-luminosity GRBs <cit.>, which likely produce little to no r-process material at all <cit.>. 
Employing the cosmic star formation history as reported by <cit.>, we obtain from Eq. (<ref>):

    n_coll = 2.6^+1.2_-1.4× 10^-5(f_b/5× 10^-3)^-1 M_⊙^-1.

This estimate also roughly applies to collapsars in our Milky Way. 

The estimate in Eq. (<ref>) is unchanged when considering that collapsars and GRBs may only occur up to a certain metallicity threshold. 
Host galaxy studies of long GRBs show that they preferentially occur below a certain stellar metallicity and thus may have shut off in recent Galactic history <cit.>. 
The metallicity threshold for collapsars is slightly sub-solar <cit.>, and is thus not relevant at the metallicities of old, metal-poor GCs considered here. 
Assuming that the estimate in Eq. (<ref>) also roughly applies to star formation in GCs, this would indicate that, given typical initial GC masses of 5×10^5 M_⊙–10^6 M_⊙, most young GCs could have of order unity collapsar events early in their cosmic histories. 
As collapsars tend to produce more r-process material per event than NSMs <cit.>, a single collapsar event per GC may be enough to explain the observed internal r-process dispersion of GCs (see Sec. <ref>). 
This simple estimate thus indicates that collapsars may indeed be a viable r-process enrichment scenario for GCs. 
Future work, however, is required to provide a better quantitative estimate. 
One open question is whether a collapsar event early in the history of the GC can lead to sufficient inhomogeneity of its r-process ejecta into a new generation of stars to explain the observed dispersion in r-process elements. While this seems conceivable at first sight, future work is needed to address this question quantitatively.

Finally, magnetohydrodynamic (MHD) SNe <cit.> might provide another source of r-process enrichment. 
If MHD SNe indeed occur at a rate of 0.1 - 0.3% of CCSNe as assumed by <cit.>, we obtain in analogy to Eq. (<ref>) n_MHDSN≈ (7× 10^-6 - 2× 10^-5) M_⊙^-1, again order unity enrichment events per GC. 
Here, we have used the observed local CCSN rate of 7.05^+1.43_-1.25× 10^-5 Mpc^-3yr^-1 <cit.>. 
However, when considering the three-dimensional stability of the magnetized jets that give rise to the fast-expanding neutron-rich ejecta in these systems <cit.>, such events are challenged to eject significant amounts of heavy r-process nuclei <cit.>. 
Additionally, if MHD SNe did produce significant amounts of heavy r-process elements, the high opacity of the lanthanide material would be mixed with the ^56Ni of the supernova ejecta in a way that would likely be incompatible with present observations of CCSNe <cit.>.






 §.§ r-process Enrichment in Other Environments


Enhancement of r-process material in metal-poor stars has been observed in various other environments, including UFDGs <cit.> and the Galactic halo <cit.>. 
Similar to GCs, the enhancement in UFDGs is difficult to explain using the NSM scenario due to their shallow gravitational potential, and depends strongly on SN kicks and distance traveled prior to merger <cit.>. 
Even with stifled natal kick velocities for DNS systems, the extremely shallow gravitational potentials of UFDGs lead to NSMs that are typically at large offsets from their hosts <cit.>. 

Using the population synthesis models from <cit.>, which use the  code <cit.>, <cit.> found that the only way NSMs can explain the enrichment in the UFDG Tucana III is if their progenitors proceed through unstable Case BB MT and form extremely hardened DNSs that merge before their post-SN systemic velocity causes them to evacuate the UFDG. 
However, the modeling in <cit.> still found it difficult to explain the enrichment in Reticulum II from this scenario. 
From our population models, we find less tension in the NSM enrichment scenario for UFDGs; for Tucana III-like and Reticulum II-like initial stellar masses and r_ vir≈ 4.5 kpc, we find that ≈ 32% and ≈ 10% of Tucana III-like and Reticulum II-like UFDGs would have an NSM enrichment event, respectively. 
This is broadly consistent with the observed number of r-processed enhanced UFDGs, which constitutes ≈ 20% of the population <cit.>. 
These numbers are nearly identical whether or not unstable (Model A) or stable (Models C and D) Case BB MT is assumed; the reason that unstable Case BB models produce about double the enrichment candidates than stable Case BB models in GCs is due to the small physical sizes of GCs (and therefore the necessity for extremely short inspiral times to merge before ejection) compared to UFDGs.
This indicates that if the NSM scenario is the main contributor to r-process enhancement in UFDGs the onset of Case BB CEs is not necessary, which is in better agreement with studies that find this phase of MT to typically proceed stably <cit.>. 






§ CONCLUSIONS


In this study, we examine the efficiency of NSMs at enriching young GCs with r-process material. 
The enrichment mechanism must occur relatively early in the history of the GC in order to pollute the second generation of stars. 
We focus on NSMs that are induced from dynamical encounters, and those that are born and evolve as isolated binary pairs in the GC. 
Our main conclusions are as follows.


    
  * Dynamically hardened DNSs merge far too late to contribute to r-process enrichment, even when the most liberal assumptions about NS segregation are assumed. 
    
    
  * For the primordial binary population to contribute enough NSM enrichment candidates, DNSs must be allowed to proceed through a phase of Case BB MT. 
    
    
  * The stability of Case BB MT affects the number of enrichment candidates in GCs by about a factor of 2; if Case BB MT typically proceeds stably we find NSM enrichment events in ∼ 15-60% of GCs, and if it typically proceeds unstably, ∼ 30-90% of GCs. 
    This is comparable to the fraction of GCs with star-to-star r-process dispersion in current observational samples <cit.>.  
    
    
    
  * Significant stripping (such that the envelope mass prior to the second SN is ≲ 0.1 M_⊙) of the donor star during Case BB MT is required for DNS enrichment candidates to remain bound to the GC.
    However, the tightest DNS progenitors, with pre-SN orbital separations of ≲ 0.2 R_⊙, will typically be ejected from the GC regardless of the amount the donor was stripped. 
    
    
  * Since most enrichment candidates are efficiently ejected from the GC and quickly evacuate the cluster environment, for NSMs to be viable enrichment candidates there must be dense, star-forming gas in the GC at the formation time of DNSs, ∼ 30-50 Myr after the initial burst of star formation. 



Though the current sample of r-process enriched GCs cannot provide a definitive enrichment fraction, future observations may be able to provide a representative sample with well understood biases. 
As we show in this paper, this can act to constrain the late mass-transfer phase of high-mass binaries — which is essential in the classical evolutionary scenario for the formation of merging DNSs — as well as the uncertain initial properties of GC progenitors. 
However, if mechanisms for quickly replenishing young GCs with star-forming gas cannot be established, another mechanism for r-process enhancement in GCs is needs to be at play. 
The predicted rate density of collapsars is in line with the stellar mass available to a young GC, though a more detailed investigation of this scenario's ability to produce the observed chemical abundances observed in r-process enriched, metal-poor GCs is required to asses this channel's viability. 






M.Z. thanks Jeff Andrews, Pablo Marchant, Mario Spera, Katie Breivik, and Chris Belczynski for useful discussions regarding binary evolution and DNS populations, Steinn Sigurdsson and Fred Rasio for insights into GC initial conditions and multiple populations, Enrico Ramirez-Ruiz for many helpful talks about r-process enrichment, Nate Bastian for useful references regarding star formation in young GCs, Claire Ye for discussions on NS populations in GCs, Nick Stone for conversations about other potential r-process scenarios, Seth Jacobson for constructive interpretations of r-process dispersion, and Chris Pankow for general all-round helpfulness. 
We also thank the anonymous referee for thoughtful and constructive comments on the content and structure of the manuscript. 
M.Z. greatly appreciates financial support from the IDEAS Fellowship, a research traineeship program supported by the National Science Foundation under grant DGE-1450006. 
C.P.L.B. is supported by the CIERA Board of Visitors Research Professorship.
V.K. is supported by a CIFAR G+EU Fellowship and Northwestern University. 
Research at Perimeter Institute is supported in part by the Government of Canada through the Department of Innovation, Science and Economic Development Canada and by the Province of Ontario through the Ministry of Economic Development, Job Creation and Trade.
The majority of our analysis was performed using the computational resources of the Quest high performance computing facility at Northwestern University, which is jointly supported by the Office of the Provost, the Office for Research, and Northwestern University Information Technology. 
This work was performed in part at Aspen Center for Physics, which is supported by National Science Foundation grant PHY-1607611. 





 <cit.>,  <cit.>,  <cit.>,  <cit.>,  <cit.>,  <cit.>.


aasjournal








§ TREATMENT OF SN KICKS FOLLOWING A CE PHASE




Besides the standard updates previously implemented in the  version of , one notable change in our population models is how SN kicks are implemented during CE evolution. 
By default,  treats this in a self-inconsistent way that leads to artificially inflated post-SN systemic velocities. 

Following a CE phase, the final separation between the core of the donor star and its compact companion is calculated using standard energetics arguments as in  <cit.>. 
The orbital hardening during a CE phase is calculated by the energy necessary to eject the envelope: 

    E_ bind, i  = α_ CE (E_ orb, f - E_ orb, i),


    -G/λM_ c(M_⋆-M_ c)/R = - α_ CEG M_ c m_ NS/2( 1/a_ f - 1/a_ i).

This can be solved for the final orbital separation as 

    a_ f = α_ CE λ R m_ NS a_ i/2 (M_⋆-M_ c) a_ i + α_ CE λ R m_ NS,

where α_ CE is the CE efficiency parameter, λ is the envelope binding energy parameter that adjusts the amount of energy needed from the orbit to eject the envelope, R is the envelope radius, and M_⋆ is the initial mass of the donor star, and M_ c is the core mass of the donor star. 

If a star evolves into an NS or BH immediately after the CE, the function for implementing the SN kick is called. 
However, in default , at this point the pre-SN mass of the exploding star is reset to its value prior to the CE (it does not account for the envelope mass that is ejected during the CE spiral-in to bring the pre-SN binary to its tightened orbital configuration), leading to a large amount of mass loss during the SN at a tight pre-SN orbital separation. 
This issue is especially drastic when dealing with CEs involving an evolved naked He star, as the cores of the stars can reach extremely hardened orbits with separations of ≈ 0.1 R_⊙. 
The contribution of the mass-loss Blaauw kick on the post-SN systemic velocity is 

    v_ sys, Blaauw≈ 220 
        (Δ M_ SN/M_⊙)
        (a_ pre/R_⊙)^-1/2[2.8 + (Δ M_ SN/M_⊙)]^-1/2km s^-1,

where Δ M_ SN is the mass lost in the SN, a_ pre is the pre-SN semi-major axis, and NS masses are assumed to be 1.4 M_⊙. 

The increased mass loss therefore leads to systematically higher post-SN systemic velocities. 
For example, a DNS progenitor at an orbital separation of 0.1 R_⊙ that loses 1.5 M_⊙ in the SN will achieve a post-SN systemic velocity solely from the Blaauw kick of ≈ 500 km s^-1. 
Therefore, even without SN natal kicks, newly born DNSs can be boosted to post-SN systemic velocities of ≳ 1000 km s^-1 solely due to mass loss in the SN, as seen in Figure <ref>. 
Though the successful Case BB CE can lead to hardened binaries that merge in ≲ 1 Myr, the amplified systemic velocities for these DNS systems cause most to escape the cluster environment before gravitational radiation leads to inspiral and merger. 
Thus, even with a successful Case BB CE, enrichment probabilities are only a few percent. 
In addition to systemic velocities, this inconsistency can strongly affect post-SN orbital properties, SN survival, and merger rates from populations that are modeled using a modified  framework. 








§ EFFECT OF SN MASS LOSS ON ENRICHMENT CANDIDATES





Small natal kicks for some NSs at birth are necessary to explain the retention fraction of NSs in GCs <cit.>. 
These retained systems are predicted to form from ECSNe or the accretion-induced collapse of a white dwarf <cit.>. 
However, retention of NSs formed in binary systems also depends on the Blaauw kick, which is anti-correlated with the square root of the binary orbital separation at the time of the SN. 

As shown in Figure <ref>, most DNSs born as a binary pair have post-SN systemic velocities that would exceed the escape velocity of young GCs. 
This is particularly apparent for systems that undergo Case BB MT, as even small amounts of mass loss in their extremely hardened state prior to SN will result in post-SN systemic velocities of ≳ 100 km s^-1. 
This does not conflict with the number of pulsars observed in GCs today since isolated NSs do not receive a Blaauw kick at formation — the natal kick alone controls the post-SN velocity and therefore whether or not it remains bound to the cluster. 

Hardened DNS systems that underwent Case BB MT make up the bulk of the systems that contribute to GC enrichment fractions. 
In their tight orbital configurations, the natal kick is typically small compared to the pre-SN orbital velocity, and the post-SN systemic velocity is mostly controlled by the mass lost in the SN rather than the natal kick that the exploding star received <cit.>. 
Lower amounts of mass loss in the SN could therefore lead to more of these systems being retained and amplifying the number of enrichment candidates (i.e. moving the purple points in the left column of Figure <ref> to the left). 

In Figure <ref> we show the impact of mass loss in the SN and pre-SN separation on the post-SN systemic velocity. 
For simplicity, we assume no natal kicks, a circular pre-SN orbit, and NSs masses of m_ NS=1.4 M_⊙, such that the post-SN systemic velocity reduces to an exact expression that is solely dependent on the orbital separation and mass loss: 

    v_ sys|_v_ k=0 = Δ M_ SN/2√(G/a_ pre (2m_ NS+Δ M_ SN)).

The entire hatched region marks the pre-SN separations and mass loss that would form a DNS that is bound to a fiducial GC and merges within 100 Myr. 
Stable Case BB MT can deplete the envelopes of the donor stars such that the mass lost in the subsequent SNe can be as low ∼ 0.1 M_⊙ and harden pre-SN orbits to a_ pre ≲ 1 R_⊙ <cit.>, safely within the hatched region. 
However, regardless of the level at which the envelope of the donor star is stripped, mass lost from neutrinos during the collapse of the stellar core results in a non-negligible contribution to the systemic velocity of the post-SN binary <cit.>, particularly in binaries that have been significantly hardened. 
The dotted–dashed line in Figure <ref> further constrains the parameter space by assuming 10% of the NS mass is lost in neutrino emission during the SN; under this assumption bound systems that are viable enrichment candidates only reside in the narrow cross-hatched region. 

Fortuitous orientations of the natal kick can increase the parameter space that will lead to DNSs with low enough systemic velocities such that they remain bound to the cluster. 
However, assuming isotropic natal kicks, this will only happen for a small number of DNS systems. 
Given how fine-tuned these conditions must be, we expect the GC enrichment contribution from bound, fast-merging DNSs to be small regardless of the physical prescriptions underlying Case BB stripping, with the bulk of the enrichment candidates coming from DNSs that are unbound from the GC and merge before evacuating the cluster environment. 

