Runaway and Hypervelocity Stars from Compact Object Encounters in Globular Clusters

The dense environments in the cores of globular clusters (GCs) facilitate many strong dynamical encounters among stellar objects. These encounters have been shown to be capable of ejecting stars from the host GC, whereupon they become runaway stars, or hypervelocity stars (HVSs) if unbound to the galactic potential. We study high-speed stellar ejecta originating from GCs by using Monte Carlo N-body models, in particular focusing on binary–single encounters involving compact objects. We pair our model-discriminated populations with observational catalogs of Milky Way GCs (MWGCs) to compose a present-day Galactic population of stellar ejecta. We find that these kinds of encounters can accelerate stars to velocities in excess of 2000 km s−1, to speeds beyond the previously predicted limits for ejecta from star-only encounters and in the same regime of Galactic center ejections. However, the same ejections can only account for 1.5%–20% of the total population of stellar runaways, and only 0.0001%–1% of HVS, with similar relative rates found for runaway white dwarfs. We also provide credible regions for ejecta from 149 MWGCs, which we hope will be useful as supplementary evidence when pairing runaway stars with origin GCs.


INTRODUCTION
Hypervelocity stars (HVSs) are stars that have been accelerated beyond the local galactic escape speed, that is, to the point of becoming unbound from the galactic potential.Theorized by Hills 1988, the classical origin of HVSs is the dynamical disruption of a stellar binary by a super massive black hole (SMBH); this "Hills mechanism" is believed to be capable of accelerating stars to speeds up to 4000 km s −1 , and naturally makes HVSs a potential probe of the galactic center.Since the first HVS discovery (Brown et al. 2005), a handful of candidate objects have been identified in the Milky Way (MW) (e.g. Brown et al. 2014), including the object S5-HVS1, with a measured speed of ∼1700 km s −1 (Koposov et al. 2020).While S5-HVS1, along with several others, is well-explained by accepting an origin from the galactic center, for which there is overwhelming evidence for an SMBH (Ghez et al. 1998, Gravity Collaboration et al. 2018, Event Horizon Telescope Collaboration et al. 2022), recent studies have found examples of high-velocity stars that are not so easily read (e.g.Boubert et al. 2018;Erkal et al. 2019;Irrgang et al. 2021).Two consistent results from these lines of work are that many stars previously classified as HVSs are nonetheless bound to the galactic potential, and many of these "runaways" are more likely to have been ejected from the disk or a satellite as opposed to the galactic center.Altogether, the identification of systems of origin for HVSs and runaway stars is a necessary threshold problem to unlocking the full potential of observations of these objects.
Alternate acceleration mechanisms capable of producing high-velocity stars include the binary supernova and dynamical ejection scenarios (BSS (Blaauw 1961) and DES (Poveda et al. 1967), respectively).In the BSS, the more massive primary of a stellar binary undergoes a supernova, which subsequently accelerates the companion; this scenario has been predicted to accelerate stars to speeds of a few hundred km s −1 (Renzo et al. 2019;Igoshev et al. 2022), with exceptionally light companions potentially receiving velocities in excess of 1000 km s −1 (Tauris 2015).The DES concerns strong gravitational interactions among 3-4 stellar objects, and has been associated with higher speed limits: Leonard (1991) studied these encounters with numerical methods, and found that the upper limit for products of these encounters was approximately the escape speed from the surface of the most massive star involved.For a sun-like star, this escape speed is ∼620 km s −1 , while for a 60 M ⊙ late main-sequence star it is ∼1400 km s −1 .Dorigo Jones et al. (2020) found that for OB runaways in the Small Magellanic Cloud, the DES is expected to occur 2-3 times more frequently that the BSS.Globular clusters (GCs) are obvious candidate matrices for both of these events due to their high stellar densities, but the DES might be exceptionally amplified due to the presence of BH subsystems in the centers of GCs.Multi-band observations (Maccarone et al. 2007;Barnard et al. 2011;Chomiuk et al. 2013;Miller-Jones et al. 2015) employing a variety of techniques (Strader et al. 2012;Bahramian et al. 2016;Giesers et al. 2018) have been successful in finding tens of BH candidates in GCs, which are in themselves consistent with the presence of hundreds of BHs per cluster (Kremer et al. 2018;Giesers et al. 2019).BHs play important roles in the macroscopic evolution of the GC (e.g.Kremer et al. 2020) by dominating gravitational interactions in GC cores; this naturally places BHs in a privileged position when considering the DES in GCs, particularly when considering the maximum velocities attainable by this mechanism.
Several recent and ongoing studies have investigated these phenomena and the broader question of extra-tidal stars from GCs, in part encouraged by the continual improvement of astrometric measurements through efforts such as Gaia (Gaia Collaboration et al. 2022).Grondin et al. (2023) developed and used a particle spray code to model ejections from GC cores and demonstrate that kinematic cuts can be too aggressive when searching for ejecta, preferring the use of chemical abundances alone.Ferrone et al. (2023) presented a comprehensive catalog of extra-tidal features of MWGCs, produced by simulating tidal stripping of known MWGCs in the context of the MW potential.Of particular proximity to our work, Weatherford et al. (2022) used the same GC model catalog we employ to holistically examine stellar ejections from GCs and identify key mechanisms.The latter two of of these works are explicitly presented as the first in a series of papers seeking to compose a more complete picture of their respective objects.
This work investigates the capacity of GCs to produce HVSs through the DES.We focus on encounters between binary and single objects, as they are the most abundant kind of strong encounter, and specifically those that involve compact objects (COs, meaning BHs, a white dwarfs (WDs), and neutron stars (NSs)).In §2, we describe our method of generating binary-single compact object (BSCO) ejecta populations from realistic simulations of GCs.In §3, we examine these populations as functions of GC parameters, and the relationship they have to GC evolution.In §4, we combine our synthetic populations with observational MWGC catalogs to produce a MW-like population of ejecta, along with predicting rates and phase space distributions for these objects.We conclude that GCs are possible generators of HVSs in all velocity regimes thus observed, beyond the previously established limit for star-only encounters, albeit the rate at which these objects are produced from GCs is significantly lower than that for the galactic center.

METHODS: SAMPLING & INTEGRATING BSCO ENCOUNTERS
The CMC Cluster Catalog was generated using CMC, a Hénon-style N -body code for collisional stellar dynamics.Developed over two decades (Joshi et al. 2000;Pattabiraman et al. 2013;Rodriguez et al. 2022b), CMC relies upon the technique originally developed by Hénon (1971a,b), where the cumulative effect of two-body encounters is modeled as an "effective" encounter between neighboring particles (in a radially sorted, spherically symmetric potential).Because these neighboring particles are individual stars (or binaries), the Hénon method allows detailed stellar and strong dynamical encounters to be considered as well.To that end CMC includes prescriptions for three-body binary formation from single BHs (Morscher et al. 2013), binary-single and binarybinary gravitational encounters using the Fewbody package (Fregeau et al. 2004;Fregeau & Rasio 2007), and galactic tidal fields (Chatterjee et al. 2013).The version of CMC used to create the CMC Cluster Catalog (which used identical physics to the public version described in Rodriguez et al. 2022b) treats single and binary stellar evolution for stars with the COSMIC code for population synthesis (Breivik et al. 2020).COSMIC is based upon the original Binary Stellar Evolution (BSE) code (Hurley et al. 2000(Hurley et al. , 2002)), but with updated prescriptions for compact-object formation and massive star evolution; see Breivik et al. (2020) and Rodriguez et al. (2022b) for details.
When one of the two neighboring particles in a cluster is a binary, CMC calculates whether to perform a strong dynamical encounter by calculating the probably P BS for an encounter to occur within a single timestep ∆T as where n is the local density of stars, w is the relative velocity between the neighboring star and binary, and Σ is the cross section for encounters to occur, given by where M is the total mass of the system and r p is the radius within which a strong encounter is assumed to occur (equal to twice the binary semi-major axis by default).
During each timestep, CMC determines whether to perform a strong three-body encounter between a neighboring star and binary by computing the probability from 1 and comparing it to a random variable drawn from [0,1], X.If X < P BS , an encounter is performed with an impact parameter, b, selected from a distribution proportional to P (b) ∼ b db out to a maximum integrated area set by 2. At this point the encounter is handed off to Fewbody (Fregeau et al. 2004), a small-N -body gravitational scattering code.Fewbody randomly selects all remaining parameters, such as the phase of the binary and the orientation of the angular momentum and Runge-Lenz vectors; of course, this means that any encounter produced in a single CMC integration is only a single realization of all the possible encounters that could have occurred in the cluster at that time.For binary-single encounters, a hyperbolic encounter between the binary centerof-mass and the single object is initiated with infinite separation, and is analytically advanced until the tidal perturbation of the binary reaches a set threshold.An 8 th -order Runge-Kutta Prince-Dormand integrator then evolves the encounter, and the participating objects are classified as singles or into binaries at regular timestep intervals.If the reduced encounter among the present single(s) and possible binary has a positive energy and the single(s)/binary are moving away from each other, then integration is terminated (termination also occurs if the encounter becomes analytic e.g. a merger results in a Keplerian interaction).
We implement additional features into Fewbody pertinent to the study of compact object dynamics.We add additional parameters to track stellar object type, distinguishing between stars and compact objects.We use these object types to decide whether to use Fewbody's default sticky-sphere collision criterion (which triggers a collision when the separation of the objects' centers is less than the sum of their radii) in the case of stellar interactions, or a tidal disruption criterion in the case of encounters involving NSs and BHs, which multiplies the threshold separation by a factor of (m CO /m * ) 1/3 (as was done in Kremer et al. 2019).To the utility of this collision metric, we conservatively assign BHs radii of 5 times their Schwarzschild radii to group the more extreme encounters with mergers, both of which require more careful treatment to accurately compute.
To establish realistic binary-single encounter populations for clusters of various parameters, we take the 148 models from the CMC Cluster Catalog and extract the initial conditions of all strong binary-single encounters that involve at least one luminous object (star, WD, or NS; to limit the focus to encounters that can produce observable ejecta) and one CO.Over the entire catalog, our encounter sample makes up about half of all strong encounters that occur in the models, with each model contributing a few ten thousand BSCO encounters on average.We realize each encounter in isolation with our modified Fewbody computing 10 realizations of each encounter while re-drawing the randomly selected parameters (e.g.binary phase) to obtain a better statistical representation of the binary-single encounter population.
The resulting objects that leave the model to become runaway stars or HVSs are identified as follows.We calculate the final velocity of an object after an encounter v fin as the hyperbolic excess velocity where v is the velocity of the object leaving the encounter, and U and K are the Keplerian potential and kinetic energies of the top-level binary-single system (all of these quantities are evaluated at the termination of integration, which as noted above is a time where the potential may still have some significant magnitude).The local escape velocity of the star cluster v esc is provided by CMC with the other encounter parameters; any object with v fin ≥ v esc is considered to escape the cluster.The velocity of an ejected object once it has left the cluster is therefore The initial conditions used for these encounters are calculated in the center-of-mass rest frame of the encounter, i.e. they do not contain information about the center-of-mass velocity of the encounter in the frame of the GC model.We do not attempt to correct for this, which leads to a minor underestimation of final velocities for encountering objects and subsequently ejection rates (we justify this in Appendix A).

Core collapse in GCs
Core collapses are well-documented features of GC evolution (e.g.Lynden-Bell & Wood 1968;Freitag & Benz 2001;Freitag et al. 2006;Binney & Tremaine 2008;Kremer et al. 2020), and the heightened densities during such phenomena makes them relevant to the present study.This process occurs due to the energy transport by two-body relaxation from the central regions of the cluster outwards, causing a core to develop and contract while the outer regions expand (Heggie & Hut 2003).The presence of massive objects (e.g.BHs) can significantly accelerate this process, as dynamical friction slows the most massive objects in the system, causing them to segregate into the core on a much more rapid timescale (Binney & Tremaine 2008).Following Breen & Heggie (2013), we refer to this initial BH-lead collapse as the "first" core collapse.While this process can be effected substantially by the metallicity of the star cluster, with more metal-rich clusters experiencing greater mass loss due to stellar winds and producing lower-mass BHs with correspondingly longer dynamical friction timescales (fewer of which are retained due to natal kicks) (Rodriguez et al. 2022a), core collapse typically occurs within a few 100 Myr for most GCs.
First core collapse is stopped by the formation and hardening of binaries in the core (which at this point is largely composed of BHs), and the liberated energy supports the universal expansion of the GC for a time on the order of Gyr (Breen & Heggie 2012).The BH population in the core is gradually depleted as strong encounters accelerate their participants beyond the local escape speed, until no more than a handful of BHs remain in the cluster (Breen & Heggie 2013;Kremer et al. 2020).The loss of the dynamical heat source causes the core to contract again, as the system adjusts to being supported by WDs, NSs and stars in the same mass regime.Unlike the first collapse, this second, permanent core collapse is observable in that the surface brightness profile of the GC becomes cuspy to the limit of the cluster center (and is what observers would classically call a "core-collapsed cluster").We follow Kremer et al. (2018) in pairing these demographic and observable features, marking the time of second core collapse in our models as the time when the number of BHs drops below 10.

Ejections from individual clusters
The 148 CMC models are designated by cluster initial size N (number of stellar objects), initial virial radius r vir (parsecs), distance from the galactic center r gc (kiloparsecs, used in CMC to calculate tidal effects from the galactic potential), and metallicity Z (used to prescribe star evolution), with values chosen to span much of the MWGC parameter space; see Kremer et al. (2020) for more details.To get a clear idea of how cluster evolution is linked to BSCO ejections, we examine four sample models; our base CMC Cluster Catalog model has initial conditions Each of the other three sample models vary one of these parameters (either N → 4 × 10 5 , r vir → 2 pc, or Z → Z ⊙ ) to promote understanding of the effects of each.The base N , r gc , and Z were chosen as they correspond to typical values for MWGCs; the notably small r vir is chosen for the samples because the dynamic range of cluster evolution and stellar ejections is greatest among CMC models of this size, making it easier to identify the same features present in larger models.
Figure 1 shows the ejection velocity v out for every MS star ejected through a BSCO encounter in these models, plotted at the time of encounter from initialization t; generally, these ejecta account for ∼1-6% of all stellar escapers from a model over its lifetime.What is immediately apparent from these single cluster models is a strong correlation between BSCO encounters and the evolution of the cluster core: whenever a cluster undergoes a core collapse, the heightened densities lead to predictably high encounter rates.
An interplay between stellar evolution and cluster dynamics is revealed in the kinds of BSCO encounters that occur at certain times.The most massive stars in a cluster are both the first to form binaries and the first to evolve into COs.Hence, encounters between a mixed binary (composed of a MS star and a CO) and a single MS star have the first opportunity to contribute significantly to the overall BSCO encounter population; in the densest models, this primacy leads to a majority of BSCO encounters being of this type over the 14 Gyr evolution.Multiple-CO encounters first become major players at first core collapse after CO formation, and can take over as the majority in the case of models with weak or nonexistent pre-core collapse BSCO encounter phases.There is a stark difference in the ejection velocities produced by these multiple-CO encounters in comparison to the previous mixed binary-single MS star encounters, and the fastest ejections over the cluster's lifetime are produced in the first tens of Myr after the associated first core collapse.This regime lasts for a few to tens of Gyr, during which the COs in the core (predominantly BHs) preferentially form CO-CO binaries, as evidenced in the drop of mixed binary-single CO encounters during this phase.
For cluster models that evolve to BH-depletion second core collapse (in our sample, all but the r vir = 2 pc model), the corresponding increase in BSCO encounters is dominated by stellar binary-single CO encounters.We also find that the majority of these encounters do not involve BHs, but rather WDs (or in some cases a NSs).These behaviors are consistent with mass segregation, where it is only after the higher-mass BHs are ejected that these lower-mass COs can migrate into the core.
This general evolution varies with model parameters, as can be observed by comparing different panels of Fig.The points are color-coded by the kind of encounter they originated from: encounters between a binary star and a CO are in red, encounters between a mixed binary (1 star and 1 CO) and a CO are in blue, encounters between a mixed binary and a star are in yellow, and encounters between a CO binary and a star are in purple.The core density (in code units) is plotted in above the scatter plot.
1.The N = 4 × 10 5 model lacks a core energetic enough to sustain a larger spatial profile, and the resulting contraction of the core leads to an ejection of 5% of the cluster mass via BSCO encounters in the first 100 Myr.This model does reach a BH-ejection core collapse by the end of the integration time, but this collapse is less pronounced than that for the base model.The r vir = 2 pc model evolves to first core collapse later than the previous two, and noticeably lacks the pre-core collapse abundance of mixed binary-single star encounters seen above.When first core collapse does occur, the distribution of encounters is similar to the respective distributions for the base model.
The difference between the solar metallicity model and the others is striking.This model does not evolve to first core collapse until ∼200 Myr (by far the latest out of the four), and leads to typical densities no more than 10% those of the other models; this occurs due to the lower BH masses and abundances associated with the greater mass loss (c.f.§3.1 above, and §5.2 in Rodriguez et al. (2022a)).Also unlike the other models, here the number of ejections does not peak until the second, BH-depletion core collapse near 7 Gyr.
The dominance of the core in overall BSCO encounter production is even more clear when considering the radial localization of the encounters.CMC does not record the radius where each strong encounter occurred, but it does record the local escape velocity.For comparison, the escape velocity from the center of the model can be computed from the central and tidal boundary potentials, which are recorded throughout the evolution of the cluster.Figure 2 plots these data for the first sample cluster, where the distribution of encounter escape velocities is strikingly correlated to the central escape velocity, indicating the degree to which the core dominates these dynamics.Throughout the evolution of the model, most BSCO ejections occur at or near the central escape velocity (the same plots for the other sample models also display this feature).There is a small fraction of points that lie above the core velocity line, which we attribute to the difference in methods employed to calculate escape velocities in the ejection and core contexts.
Figure 3 shows histograms for cluster ejection velocities v out and masses m of all stars ejected by BSCO encounters in all CMC Cluster Catalog models, binned by model parameters N , r vir , and Z.The number of CMC models corresponding to each parameter value varies, and so each histogram is averaged over the respective models, in addition to being divided by the factor of 10 in encounter multiplication.It is worth noting the prominence of the N = 4×10 5 , r vir = 0.5 pc, Z = 0.0002 models in these averaged histograms: they are the cause of the peaks above the other distributions, as the three (one for each value of r gc ∈ {2, 8, 20} kpc) produce significantly more ejections than the other models (∼2.9, 2.0, and 1.4 × 10 4 , respectively, versus the average number of 1 − 2 × 10 3 for the catalog).
Acknowledging these especially fecund models, predicable trends are visible in these histograms: increasing mass and decreasing size are both associated with higher ejection velocities.However, while the number of ejections increases as the models become more compact, increasing the number of particles/mass of the model leads to a slight decrease in ejecta.The N = 4 × 10 5 models consistently produce the most ejections in comparison with otherwise identical models of different population sizes, suggesting that at this mass there are not enough massive objects to prevent a severe initial core collapse through binary burning, but there are enough to facilitate a high number of BSCO encounters (the number of strong encounters in general is also maximized in the N = 4 × 10 5 models).
The most noticeable distinction among models of different metallicities is the number of ejections, which decreases with increasing metallicity.This is understandable, as in practice the lower metallicity models have greater BH populations than the higher metallicity mod- In the the mass histograms, the data are further divided by whether the ejection occurred before or after the "second," BH-depletion core collapse of the cluster, if one occurred within the integration time.
els; quantitatively, the 1% solar models have about 1.2× as many BHs at the time of first core collapse as the solar models, and maximum BH masses about 2× those of the same.This weakens BSCO dynamics before and during first core collapse, as was first made visible in Figure 1.Because this effect takes place in the earlier stages of model evolution -namely, when mass segregation has had less time to separate lighter objects from the strong dynamics of heavier objects -this has the most significant effect of reducing the number of ejecta with m < M ⊙ .Notably, the different metallicity models have quite similar ejecta mass distributions after second core collapse, as the effect of metallicity on the masses of the remaining WDs and NSs is much less pronounced than on the now-ejected BHs.

A MW-LIKE POPULATION OF BSCO EJECTA
In the interest of predicting realistic statistics and rates for a MW-like GC population, we seek to assemble a synthetic MW-like population of GC-runaways from BSCO encounters.The two steps involved here are 1) selecting representative CMC models for galactic GCs, and 2) integrating the post-ejection orbits in the context of the MW to produce a present-day picture.

Pairing CMC models to MWGCs
We predominantly use the observational catalog of Baumgardt & Hilker (2018) to obtain parameters for MWGCs.This catalog lacks metallicity measurements for the objects; therefore we supplement with that of Harris (2010).12 of the GCs in the former catalog do not have metallicity measurements in the latter, leaving us with the 149 MWGCs we use in this analysis.
Assigning a representative model to each MWGC is nontrivial.While many GCs are expected to be older than 12 Gyr, some could be as young as 9 Gyr (Vanden-Berg et al. 2013); to reflect this, we find 11 timesteps as evenly spaced as possible between 10 and 13.5 Gyr for each CMC model, yielding a set of model snapshots spanning the different CMC initial conditions and sampling the models at different evolutionary stages.
These snapshots are plotted with our composite catalog of MWGCs in Figure 4.The upper plot shows the models/GCs in r gc -[Fe/H] space; note the discretization of the CMC models, as the respective parameters are held constant over integration, and so each blue point represents the ∼ 15 CMC models that are initiated with those particular values.When choosing a representative snapshot for each GC, we effectively follow Rui et al. (2021) by first finding the blue point closest to the GC in this space and constrain our search to the respective snapshots.
Having constrained the snapshots by r gc and metallicity, we then choose one of an appropriate size and evolutionary state by comparing the snapshots with the model in normalized log M -r c /r h space, where M is the mass of the cluster in M ⊙ , r c is the Spitzer core radius (Spitzer 1987) (where σ c is the central velocity dispersion and ρ c is the central density) and r h is the 3D half-mass radius of the snapshot (the latter two parameters are used because both are immediately accessible in the Baumgardt & Hilker (2018) catalog and standard CMC output).The models and GCs are plotted in the non-normalized space of these parameters in the bottom plot of Figure 4.The transformation to the normalized space is simply  practice, we find that the models that are chosen for several MWGCs are average models for the CMC Cluster Catalog, and so we expect that while characteristics of more extreme GCs are not necessarily well-represented, they are replaced by typical models nonetheless.

Integrating runaways to the present day
Having chosen representative snapshots for all MWGCs, we now use our synthetic ejecta populations to compose a MW-like population of BSCO ejecta from GCs.All orbit integrations in the galactic frame are done with galpy1 (Bovy 2015).Setting the snapshot to the present-day galactic orbital parameters from Baumgardt & Hilker (2018), we back-integrate the orbit of each GC in the galpy MWPotential2014 potential to the initialization time of the respective model.We then locate each ejected object at the appropriate place in the GC orbit, using the time of BSCO encounter (naturally ignoring any ejections that occurred after the time of the chosen snapshot).A random ejection direction is chosen in the GC rest frame for each object, and the velocities are then transformed to the galactic rest frame.Finally, the orbits of all ejecta are then integrated to the present day in the same galactic potential.
Figure 5 shows the integrated orbits for a selection of GCs, and the points at which objects are ejected from the cluster.As was seen earlier, ejections occur much more frequently in the early stages of the cluster, and here the slower velocities at larger distances from the galactic center leads to a higher percentage of objects being ejected far away from this focus.
The present-day synthetic populations of runaways for the same sample CMC/GC pairings are shown in Figure 6.While it is clear how the GC orbits influence the distribution of ejected objects, the ejecta wander to a broad enough spread that any one object loses some of the information of its cluster of origin by the present day.The proper motion distribution is less affected in this way (especially for more circular GC orbits), which is coherent with previous works that use these velocities to study the origins of such objects.What is important to note here is that the current proper motion of the cluster is not necessarily the best locus to use when comparing stellar proper motions: a well-informed back-integration of the orbit can reveal the average proper motion of the cluster, which may be distinguishable from its presentday proper motion.We include histograms, and 50% and 90% credible regions in sky area and proper motion space as a product of this work, for use as evidence when assigning runaway stars to GC origins; see Appendix B for a description of these products.
From a galactocentric perspective, the population of runaways is fairly isotropic in position, as can be seen in Figure 7.The distributions of distance from the galactic origin, and the local distributions of galactocentric velocity are similar when distributing over radius r gc versus distance from the galactic plane Z gc .Most runaways end up at a distance on the order of 10 kpc from the galactic center, and with a velocity on the order of 100-300 km s −1 in the galactocentric rest frame.Few runaways make it past the 100 kpc mark, but those that do naturally retain the highest velocities, reaching upwards of 1000 km s −1 .
One result of note is that of the heliocentric radial velocity, v rf , distribution of the synthetic population.Generozov & Perets (2021) in part studied the galactic center origin of HVSs, and found that there is an apparent tension between the observed and predicted runaway populations; specifically, the predicted number of stars with velocities ≳ 700 km s −1 was much higher than the observed rate from the HVS sample of Brown et al. (2018).The same work found that a decreasing star for-mation rate over last millions of years would reduce the recent (and thus more likely to be observed) HVS production rate enough to overcome the discrepancy.GCs naturally follow this pattern (c.f.§3), as the BSCO runaways studied here are produced at much higher rates when the models are young, and few HVSs are produced in the later stages of the cluster.The resulting v rf distribution is much closer to the observed distribution than that for the galactic center origin case, albeit skewed towards slightly smaller velocities (Figure 8).This latter difference grows when comparing to the HVS sample of Hattori et al. (2018), who focused on metal-poor stars, which are more akin to the objects that populate GCs.
Finally, Figure 9 shows the HVS/runaway rates for our synthetic population, along with observational rate estimates from Brown (2015).We note that 50(90)% of ejections occur before our synthetic population is 3.21(8.59)Gyr old, and that naturally the earliest ejections are most dependent on any time variability in the MW potential.At these earlier times, BSCO ejecta could have made as much as 20% (1%) of all runaways (HVSs), during the times when GCs underwent their first core collapses.2Afterwards, the ejection rates decay to the present day, where they are closer to 1.5% (0.0001%) of the observational values.It is important to note nonetheless that GCs are able to produce HVSs with speeds high enough to make them comparable to those accelerated by mechanisms in the galactic center.

DISCUSSION & CONCLUSIONS
In this work we have studied strong 3-body encounters in GCs as means of producing stellar runaways.We composed a synthetic MW-like population of ejecta by matching observed GCs to realistic N -body models of these systems and embedding the models in the orbits of their real counterparts.In particular we considered binary-single encounters involving at least one compact object; this selection includes about half of all strong encounters in the catalog of models.
BSCO ejections were found to be closely linked to the evolution of the cluster core, where the closest encounters among the densest stellar objects occur.GCs lose mass and expand as they evolve; accordingly, the majority of and the fastest BSCO ejecta were produced in the early stages of the models.High-metallicity models had overall weaker ejection mechanisms (in frequency and maximum velocity), due to smaller stellar masses.The back-integrated orbits are shown as the gray trajectories, and the blue "x" is the position/velocity of the GC as measured by Baumgardt & Hilker (2018).The set of synthetic ejecta shown here is the result of downsampling the total set by a factor of ten, to account for the repeated-realizations method described in §2.BSCO encounters occurring in realistic GCs are capable of accelerating stars to velocities in excess of 2000 km s −1 , which complicates the identification of ejection mechanisms for HVSs when their origins are not easily recognizable.We also note that these velocities appear to pass the speed limit on star-only encounters found The evolution of the ejection rate for different components of our synthetic population.The blue (red) dashed line shows the estimated from observations for runaways from the galactic disc (HVSs from the galactic center), as reported by Brown (2015).The threshold times by which 50% and 90% of the runaway ejections have occured are indicated along the bottom of the figure .by Leonard (1991); further study is warranted to dis-cern the impact of compact objects on small-N -body dynamics.
While ejected objects evolve to be largely indistinguishable from other MW stars in terms of position, they were found capable of retaining some information about the motion of their GC of origin, particularly in the case of GCs with near-circular orbits.The overall population of ejecta was usually concentrated around the average proper motion of the GC throughout its orbit.It is important to recognize that the present-day proper motion of a GC may not reflect this average proper motion, and that in general a better kinematic picture is accessible through back-integration of the orbit.
In the galactic context, the velocity distribution of the synthetic ejecta was found to be similar to that of HVS observations for velocities ≲ 500 km s −1 ; specifically, our population was skewed towards these velocities with respect to observations.With galactic-center origin studies finding distributions skewed towards higher velocities in the same respect (Brown et al. 2018), it is possible that a mixture of the two could be used to more accurately model the real population of these objects.Such a calculation must be done in light of the relative rates of the two mechanisms: our study concludes that while the GC BSCO runaway rate might have been a few 10% of the overall rate in the first few Gyr of the MW, in the present day it is likely no more a few 1% of the same.This work was supported by NSF Grant AST-2009916 and NASA ATP Grant 80NSSC22K0722.CR also acknowledges support from a Charles E. Kaufman Foundation New Investigator Research Grant, an Alfred P. Sloan Research Fellowship, and a David and Lucile Packard Foundation Fellowship.This work was prepared in part with the reproducibility software showyourwork (Luger et al. 2021), which leverages continuous integration to transparently connect the paper to the publicly-available dataset and relevant code.The GitHub repository for this project can be found at https://github.com/tomas-cabrera/hvss-bsco; the icons next to paper content link to the particular scripts used to generate that content.The dataset is stored at https://doi.org/10.5281/zenodo.7599870,which also includes copies of the simulation and figure scripts from the GitHub for insurance of future access.

A. VELOCITIES IN ENCOUNTER REST FRAMES VERSUS THE GC REST FRAME
We justify here our claim that in neglecting the transformation from the encounter rest frame to the GC model rest frame we obtain slower -and subsequently fewer -ejections.
Hénon's Monte Carlo method assumes spherical symmetry, and subsequently the orbit of a particle is characterized by its energy, total angular momentum, and radial position; radial and tangential velocity magnitudes are drawn from the orbit by weighting by the amount of time the object spends at each point in its orbit (see §2.5 of Rodriguez et al. (2022b)).The sign of the radial velocity is randomly chosen as positive or negative with equal weighting, and when setting up a strong interaction an angle 0 ≤ ϕ ≤ 2π between the tangential velocities of the two objects is chosen from a uniform distribution.These two features ensure isotropy of either velocity with respect to the other, and of the center-of-mass velocity of the encounter with respect to the GC model rest frame ⃗ v cm|GC .
If we assume that the direction of the post-encounter object velocity in the center-of-mass rest frame ⃗ v f|cm is isotropic with respect to ⃗ v cm|GC , then the average speed after boosting back to the model rest frame is Naturally, for v cm|GC ≪ v f|cm this average post-boost speed approaches v f|cm , and in the limit v cm|GC ≫ v f|cm it approaches v cm|GC , which recall must be less than the initial speed of the other object in the encounter: in either of these cases, the speed augmentation caused by the boost does not favor faster or slower post-encounter velocities.The maximum average "acceleration" resulting from the boost occurs at the limit v cm|GC = v f|cm , where ⟨v f|GC ⟩ ∼ 1.3v f|cm = 1.3v cm|GC .The assumption of isotropy in calculating ⟨v f|GC ⟩ is appropriate for resonant encounters wherein the dynamics are chaotic.For flyby encounters where the objects travel on roughly hyperbolic trajectories, there is a preference for postencounter velocities in the same direction as the initial velocity, if the impact parameter distribution is sufficiently expansive and weighted by the square of the parameter and the final velocities are marginalized over the angle that orients the plane of the 2-body encounter between the star and the binary center-of-mass.The isotropic case is therefore a conservative limit on this flyby case, and predicts a greater boost-propagated speedup than if the calculation was done in full detail.
In summary, neglecting to return to the model rest frame after the Fewbody step for strong encounters does not favor faster or slower post-encounter speeds when the final speed of an object is much greater than or much less than the center-of-mass speed of the encounter in the model rest frame.In the case that the final speed of an object is comparable to the same center-of-mass speed, the average post-transformation boost of 30% is unlikely to significantly increase the number of ejections, as the two speeds are generally less than the maximum initial speed among objects entering the encounter.

B. EJECTA CREDIBLE REGIONS
We include with this publication phase credible regions in phase space for ejecta populations from the MWGCs considered, which may be accessed at https://doi.org/10.5281/zenodo.7599870.We intend these constructions to be used as secondary evidence for cluster membership of runaway stars, and to encourage further studies into understanding the respective phase space distributions.The 2D position and proper motion distributions are separated into two different files, the conventions of which are described below.
The distribution of ejecta on the sky is quantified by discretizing the sphere with an order 4 nested HEALPix map (Górski et al. 2005), binning the sky into 3072 equal-area tiles.We choose this resolution to enable identification of interesting credible regions while minimizing the apparent effect of isolated points whose exact location is dependent on the RNG seeds used.We create the respective histogram by counting the number of synthetic ejecta that are found in each tile in the present day, and normalize by the total number of ejecta for the GC.The resulting probability histogram is stored in the first column (PROBS) of the hp probs.fitsfile for each GC; the NEJECT field in the header of the same file contains the number of ejections for the GC.We calculate our credible regions by cumulatively adding the highest probability bins until the target percentage is reached.The last two columns (CR50 and CR90) are boolean masks of the same convection as the HEALPix probability histogram corresponding to the credible regions (50% and 90%, respectively), where bins with entries of 1 are included in the region.Figure 10 shows the histogram and 50% and 90% credible regions for NGC 104, as an example (the script used to generate this plot is check credible regions.py,which may be found in the GitHub/Zenodo repositories).
We construct proper motion histograms and credible regions in a similar manner.We use a domain of −30 ≤ µ α cos δ [km s −1 ] ≤ 30, −45 ≤ µ δ [km s −1 ] ≤ 15 divided into a 50 × 50 grid.These bounds and resolution are included in the headers of the pm prob.fitsfiles for each GC, where calling np.linspace(PM[D/RCD]MIN, PM[D/RCD]MAX, PMNUM) will return the bin edges used for the appropriate dimension.The header also includes a COVERAGE field containing the fraction of ejecta that lie in the specified domain; for all GCs this fraction is at least 0.98, and in most cases is greater than 0.999.The total number of ejected objects (including those outside of the domain) is stored in the NEJECT field, as for the HEALPix histograms.The proper motion histogram and respective masks for the 50% and 90% credible regions are stored in the pm prob.fitsfiles as separate HDUs; these items for the same example GC as the HEALPix plot is show in Figure 11, and the same check credible regions.pyscript contains the generating code.

Figure 1 .
Figure 1.Scatter plots of the cluster ejection velocity vout versus encounter time t for every escaping object from the integrated encounters for the four sample CMC models (see the beginning of §3.2 for details); the histograms show the distribution of velocities.The points are color-coded by the kind of encounter they originated from: encounters between a binary star and a CO are in red, encounters between a mixed binary (1 star and 1 CO) and a CO are in blue, encounters between a mixed binary and a star are in yellow, and encounters between a CO binary and a star are in purple.The core density (in code units) is plotted in above the scatter plot.

Figure 2 .
Figure 2. The local escape velocities vesc at the locations of all BSCO encounters generated from the first sample CMC model (colored points), plotted alongside the central escape velocity of the model (dark blue).The top plot shows the data in units of km s −1 , while the lower plot shows the same data normalized to the continuum of the core escape velocity (dark blue line in the upper plot).See Figure 1 for an explanation of the color scheme.

Figure 3 .
Figure 3. Histograms for all MS stars ejected from the CMC Cluster Catalog models as a result of BSCO encounters.The top (bottom) row displays the distribution of ejection velocities from the models vout (masses m of the ejected objects).Each each separates the data by different CMC model parameters: either size N (number of objects), initial virial radius rvir (parsecs), or metallicity Z.Each histogram is averaged over all models computed with the respective value of model parameter.In the the mass histograms, the data are further divided by whether the ejection occurred before or after the "second," BH-depletion core collapse of the cluster, if one occurred within the integration time.
where xCMC and σ x,CMC are respectively the mean and standard deviation of the parameter x over the complete set of CMC snapshots.The snapshot from the r gc -[Fe/H] subset that is closest to the MWGC in the normalized log M -r c /r h space is chosen to represent the GC; as an example, 47 Tuc (r gc = 7.52 kpc, Z = 0.0019; log M = 5.95, r c /r h , m = 0.125 (Baumgardt & Hilker 2018; Harris 2010)) is represented by a snapshot from model N1.6e6 rv2 rg8 Z0.002, which at the 10.8 Gyr time of the snapshot has log M = 5.91 and r c /r h = 0.126.While the CMC Cluster Catalog has fairly representative models for most of the MWGCs, there are a number of larger-cored GCs (right side of the bottom plot of Figure4) that are relatively distant from the nearest CMC model; in practice, this means that a single CMC model can be used to represent several MWGCs.In Figure 4.The clusters and models used in this work.Galactic GCs are represented by black triangles, and CMC model snapshots by blue dots.The top (bottom) plot shows the systems in rgc v. [Fe/H] (log M v. rc/r h ) space.

Figure 5 .Figure 6 .
Figure5.Plots showing the back-integrated orbits for some sample MWGCs (gray curves), and the points in the orbit where an object is ejected from the GC (scatter points).The color scale communicates the mass of the ejected star.The concentration of ejections in the first few orbits is clear, and an increased density of ejections when clusters are farther away from the galactic center is visible as well.Viewing the figure electronically makes it easier to find the few high-mass ejecta amid the abundance of lower mass objects.

]Figure 7 .Figure 8 .
Figure 7. Histograms (top) and velocity quantiles (bottom) for the synthetic stellar ejecta.The left (right) plots show the profile over radial distance from the galactic center rgc (distance from the galactic plane Z).The quantiles in the lower plots are calculated from the present-day velocities of our population.The red dots show the MWGCs from our composite Baumgardt & Hilker (2018)+Harris (2010) catalog.
Figure 9.The evolution of the ejection rate for different components of our synthetic population.The blue (red) dashed line shows the estimated from observations for runaways from the galactic disc (HVSs from the galactic center), as reported byBrown (2015).The threshold times by which 50% and 90% of the runaway ejections have occured are indicated along the bottom of the figure.

Figure 10 .
Figure 10.The HEALPix histogram for a sample GC (NGC 104), with the 50% and 90% credible regions in the right subplot.Galactic coordinates are used, and the histogram weights are log scaled.
Figure 11.The proper motion histogram and credible regions for the same sample GC.