# Runaway and Hypervelocity Stars from Compact Object Encounters in Globular
Clusters

Tomás Cabrera McWilliams Center for Cosmology, Department of Physics, Carnegie
Mellon University, 5000 Forbes Avenue, Pittsburgh, PA 15213  Carl L. Rodriguez
McWilliams Center for Cosmology, Department of Physics, Carnegie Mellon
University, 5000 Forbes Avenue, Pittsburgh, PA 15213  Department of Physics
and Astronomy, University of North Carolina at Chapel Hill, 120 E. Cameron
Ave, Chapel Hill, NC, 27599, USA

###### Abstract

The dense environments in the cores of globular clusters (GCs) facilitate many
strong dynamical encounters among stellar objects. These encounters have been
shown capable of ejecting stars from the host GC, whereupon they become
runaway stars, or hypervelocity stars 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 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 ${\rm 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 hypervelocity stars, with similar relative
rates found for runaway white dwarfs. We also provide credible regions for
ejecta from 149 Milky Way GCs, which we hope will be useful as supplementary
evidence when pairing runaway stars with origin GCs.

††software: CMC (Joshi et al., 2000; Pattabiraman et al., 2013; Rodriguez et
al., 2022b), Fewbody(Fregeau et al., 2004), galpy (Bovy, 2015), astropy
(Astropy Collaboration et al., 2013, 2018, 2022), numpy (Harris et al., 2020),
pandas (Wes McKinney, 2010; pandas development team, 2020), matplotlib
(Hunter, 2007)

##  1 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 $\sim$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 ${\rm km~{}s^{-1}}$(Renzo et al., 2019;
Igoshev et al., 2022), with exceptionally light companions potentially
receiving velocities in excess of 1000 ${\rm 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
$\sim$620 ${\rm km~{}s^{-1}}$, while for a 60 $M_{\odot}$ late main-sequence
star it is $\sim$1400 ${\rm km~{}s^{-1}}$. Oey et al. (2023) recently 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.

##  2 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
binary-binary 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, 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_{\rm BS}$ for an encounter to occur within a single timestep
$\Delta T$ as

$P_{\rm BS}=n\Sigma w\Delta T~{},$ (1)

where $n$ is the local density of stars, $w$ is the relative velocity between
the neighboring star and binary, and $\Sigma$ is the cross section for
encounters to occur, given by

$\Sigma=\pi r_{p}^{2}\left(1+\frac{2GM}{r_{p}w^{2}}\right),$ (2)

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_{\rm BS}$, an encounter is performed with an impact
parameter, $b$, selected from a distribution proportional to $P(b)\sim 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
center-of-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 8th-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, and a tidal disruption criterion in the case of encounters
involving NSs and BHs, which multiplies the threshold separation by a factor
of $(m_{\rm 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 represent.

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 then 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_{\rm fin}$ as the hyperbolic excess velocity

$v_{\rm fin}=v\sqrt{\frac{U+K}{K}},$ (3)

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_{\rm esc}$ is provided by CMC with the other encounter parameters; any
object with $v_{\rm fin}\geq v_{\rm esc}$ is considered to escape the cluster.
The velocity of an ejected object once it has left the cluster is therefore

$v_{\rm out}=\sqrt{v_{\rm fin}^{2}-v_{\rm esc}^{2}}.$ (4)

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.

##  3 GC dynamics and BSCO ejections

###  3.1 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 timescaleBinney & 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 same time the
number of BHs drops below 10.

###  3.2 Ejections from individual clusters

The 148 CMC models are designated by cluster initial size $N$ (number of
stellar objects), initial virial radius $r_{\rm vir}$ (parsecs), distance from
the galactic center $r_{\rm 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

  * •

$N=8\times 10^{5}$

  * •

$r_{\rm vir}=0.5~{}{\rm pc}$

  * •

$r_{\rm gc}=8~{}{\rm kpc}$

  * •

$Z=0.01Z_{\odot}$.

Each of the other three sample models vary one of these parameters (either
$N\to 4\times 10^{5}$, $r_{\rm vir}\to 2~{}{\rm pc}$, or $Z\to Z_{\odot}$) to
promote understanding of the effects of each. The base $N$, $r_{\rm gc}$, and
$Z$ were chosen as they correspond to typical values for MWGCs; the notably
small $r_{\rm 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_{\rm out}$ for every MS star ejected
through a BSCO encounter in these models, plotted at the time of encounter
from initialization $t$. 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.

Figure 1:  Scatter plots of the cluster ejection velocity $v_{\rm out}$ 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.

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 tend to dominate the first generations of BSCO
encounters. This dominance usually continues until the first core collapse
after CO formation, where the heightened densities facilitate a rise in
encounters involving multiple COs. There is a stark difference in the ejection
velocities produced by these 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_{\rm vir}=2~{}{\rm 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. 1. The $N=4\times 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_{\rm vir}=2~{}{\rm 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 $\sim$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 same sample clusters, 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.

Figure 2:  The local escape velocities $v_{\rm esc}$ of all BSCO encounters
generated from the four sample CMC models (colored points), plotted alongside
the central escape velocity of the model (maroon). See Figure 1 for an
explanation of the color scheme.

Figure 3 shows histograms for cluster ejection velocities $v_{\rm out}$ and
masses $m$ of all stars ejected by BSCO encounters in all CMC Cluster Catalog
models, binned by model parameters $N$, $r_{\rm 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\times 10^{5}$, $r_{\rm 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_{\rm gc}\in\\{2,8,20\\}$
kpc) produce significantly more ejections than the other models ($\sim$2.9,
2.0, and $1.4\times 10^{4}$, respectively, versus the average number of
$1-2\times 10^{3}$ for the catalog).

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 $v_{\rm out}$ (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 $r_{\rm
vir}$ (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.

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\times 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\times 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 models; quantitatively, the 1% solar
models have about $1.2\times$ as many BHs at the time of first core collapse
as the solar models, and maximum BH masses about $2\times$ 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_{\odot}$. 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.

##  4 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.

###  4.1 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
(VandenBerg 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_{\rm 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 $\sim 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_{\rm 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_{\odot}$, $r_{c}$ is the Spitzer core
radius (Spitzer, 1987)

$r_{c}=\sqrt{\frac{3\sigma_{c}^{2}}{4\pi\rho_{c}}},$ (5)

(where $\sigma_{c}$ is the central velocity dispersion and $\rho_{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

$x_{\rm norm}=\frac{x-\bar{x}_{\rm\texttt{CMC}}}{\sigma_{\rm
x,\texttt{CMC}}},$ (6)

where $\bar{x}_{\rm CMC}$ and $\sigma_{\rm 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_{\rm 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_{\rm 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$.

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 $r_{\rm gc}$ v. [Fe/H] ($\log M$ v.
$r_{c}/r_{h}$) space.

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 Figure 4) 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 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.

###  4.2 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
galpy111http://github.com/jobovy/galpy (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 galpyMWPotential2014 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.

Figure 5:  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.

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 present-
day 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.

Figure 6:  Present-day positions (galactic longitude/latitude) and velocities
(projected onto the galactic longitude/latitude directions) for the runaway
objects from the sample GCs. The color scale is the same as in Figure 5,
depicting the eject 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.

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_{\rm gc}$
versus distance from the galactic plane $Z_{\rm 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 ${\rm 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 ${\rm km~{}s^{-1}}$.

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 $r_{\rm gc}$ (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 $r_{\rm gc}$ ($Z$)
positions of the MWGCs from our composite Baumgardt & Hilker (2018)+Harris
(2010) catalog, arbitrarily spread in the vertical dimension.  Figure 8:
Velocity distributions for our synthetic and some previous observational HVS
catalogs. $v_{\rm rf}$ is the heliocentric radial velocity of the star as
measured in the galactic rest frame.  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 by Brown (2015).

One result of note is that of the heliocentric radial velocity, $v_{\rm 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 $\gtrsim
700~{}{\rm 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
formation 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_{\rm 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). At 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.222Here, the cutoff of
this phase is fairly abrupt because we constrain our models to present-day
ages of 10-13.5 Gyr (§4.1). The physical rate will naturally be characterized
by the evolution of the GC formation rate in the past. Afterwards, 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.

##  5 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. BSCO encounters occurring in
realistic GCs are capable of accelerating stars to velocities in excess of
2000 ${\rm 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 by Leonard (1991); further study is warranted to discern 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 $\lesssim
500~{}{\rm 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.

##  Appendix 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\leq\phi\leq 2\pi$ 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 $\vec{v}_{\rm cm|GC}$.

If we assume that the direction of the post-encounter object velocity in the
center-of-mass rest frame $\vec{v}_{\rm f|cm}$ is isotropic with respect to
$\vec{v}_{\rm cm|GC}$, then the average speed after boosting back to the model
rest frame is

$\langle v_{\rm f|GC}\rangle=\int_{0}^{2\pi}\frac{d\phi}{2\pi}\sqrt{(v_{\rm
f|cm}+v_{\rm cm}\cos\phi)^{2}+(v_{\rm cm}\sin\phi)^{2}}.$ (A1)

Naturally, for $v_{\rm cm|GC}\ll v_{\rm f|cm}$ this average post-boost speed
approaches $v_{\rm f|cm}$, and in the limit $v_{\rm cm|GC}\gg v_{\rm f|cm}$ it
approaches $v_{\rm 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_{\rm cm|GC}=v_{\rm f|cm}$, where $\langle v_{\rm
f|GC}\rangle\sim 1.3v_{\rm f|cm}=1.3v_{\rm cm|GC}$.

The assumption of isotropy in calculating $\langle v_{\rm f|GC}\rangle$ 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 post-encounter 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.

##  Appendix B Ejecta credible regions

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.

We include with this publication phase credible regions in phase space for
ejecta populations from the MWGCs considered, which may be accessed at
https://zenodo.org/record/7599871. 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.fits file 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 in the GitHub repository;
click on the icon by the caption to be taken directly to the file).

We construct proper motion histograms and credible regions in a similar
manner. We use a domain of
$-30\leq\mu_{\alpha}\cos\delta~{}{\rm[km~{}s^{-1}]}\leq 30$,
$-45\leq\mu_{\delta}~{}{\rm[km~{}s^{-1}]}\leq 15$ divided into a 50 $\times$
50 grid. These bounds and resolution are included in the headers of the
pm_prob.fits files 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.fits files 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.py
script contains the generating code.

## References

  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33, doi: 10.1051/0004-6361/201322068
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123, doi: 10.3847/1538-3881/aabc4f
  * Astropy Collaboration et al. (2022) Astropy Collaboration, Price-Whelan, A. M., Lim, P. L., et al. 2022, ApJ, 935, 167, doi: 10.3847/1538-4357/ac7c74
  * Bahramian et al. (2016) Bahramian, A., Heinke, C. O., Tudor, V., et al. 2016, MNRAS, 000, 1
  * Barnard et al. (2011) Barnard, R., Garcia, M., Li, Z., Primini, F., & Murray, S. S. 2011, The Astrophysical Journal, 734, 79
  * Baumgardt & Hilker (2018) Baumgardt, H., & Hilker, M. 2018, MNRAS, 478, 1520, doi: 10.1093/mnras/sty1057
  * Binney & Tremaine (2008) Binney, J., & Tremaine, S. 2008, Galactic Dynamics: Second Edition (Princeton University Press)
  * Blaauw (1961) Blaauw, A. 1961, Bull. Astron. Inst. Netherlands, 15, 265
  * Boubert et al. (2018) Boubert, D., Guillochon, J., Hawkins, K., et al. 2018, MNRAS, 479, 2789, doi: 10.1093/mnras/sty1601
  * Bovy (2015) Bovy, J. 2015, ApJS, 216, 29, doi: 10.1088/0067-0049/216/2/29
  * Breen & Heggie (2012) Breen, P. G., & Heggie, D. C. 2012, MNRAS, 420, 309, doi: 10.1111/j.1365-2966.2011.20036.x
  * Breen & Heggie (2013) —. 2013, MNRAS, 432, 2779, doi: 10.1093/mnras/stt628
  * Breivik et al. (2020) Breivik, K., Coughlin, S., Zevin, M., et al. 2020, ApJ, 898, 71, doi: 10.3847/1538-4357/ab9d85
  * Brown (2015) Brown, W. R. 2015, ARA&A, 53, 15, doi: 10.1146/annurev-astro-082214-122230
  * Brown et al. (2014) Brown, W. R., Geller, M. J., & Kenyon, S. J. 2014, ApJ, 787, 89, doi: 10.1088/0004-637X/787/1/89
  * Brown et al. (2005) Brown, W. R., Geller, M. J., Kenyon, S. J., & Kurtz, M. J. 2005, ApJ, 622, L33, doi: 10.1086/429378
  * Brown et al. (2018) Brown, W. R., Lattanzi, M. G., Kenyon, S. J., & Geller, M. J. 2018, ApJ, 866, 39, doi: 10.3847/1538-4357/aadb8e
  * Chatterjee et al. (2013) Chatterjee, S., Umbreit, S., Fregeau, J. M., & Rasio, F. A. 2013, MNRAS, 429, 2881, doi: 10.1093/mnras/sts464
  * Chomiuk et al. (2013) Chomiuk, L., Strader, J., Maccarone, T. J., et al. 2013, The Astrophysical Journal, 777, 69, doi: 10.1088/0004-637X/777/1/69
  * Erkal et al. (2019) Erkal, D., Boubert, D., Gualandris, A., Evans, N. W., & Antonini, F. 2019, MNRAS, 483, 2007, doi: 10.1093/mnras/sty2674
  * Event Horizon Telescope Collaboration et al. (2022) Event Horizon Telescope Collaboration, Akiyama, K., Alberdi, A., et al. 2022, ApJ, 930, L12, doi: 10.3847/2041-8213/ac6674
  * Ferrone et al. (2023) Ferrone, S., Di Matteo, P., Mastrobuono-Battisti, A., et al. 2023, arXiv e-prints, arXiv:2301.05166, doi: 10.48550/arXiv.2301.05166
  * Fregeau et al. (2004) Fregeau, J. M., Cheung, P., Portegies Zwart, S. F., & Rasio, F. A. 2004, MNRAS, 352, 1, doi: 10.1111/j.1365-2966.2004.07914.x
  * Fregeau & Rasio (2007) Fregeau, J. M., & Rasio, F. A. 2007, ApJ, 658, 1047, doi: 10.1086/511809
  * Freitag & Benz (2001) Freitag, M., & Benz, W. 2001, A&A, 375, 711, doi: 10.1051/0004-6361:20010706
  * Freitag et al. (2006) Freitag, M., Rasio, F. A., & Baumgardt, H. 2006, MNRAS, 368, 121, doi: 10.1111/j.1365-2966.2006.10095.x
  * Gaia Collaboration et al. (2022) Gaia Collaboration, Vallenari, A., Brown, A. G. A., et al. 2022, arXiv e-prints, arXiv:2208.00211, doi: 10.48550/arXiv.2208.00211
  * Generozov & Perets (2021) Generozov, A., & Perets, H. B. 2021, arXiv e-prints, arXiv:2112.13864.  https://arxiv.org/abs/2112.13864
  * Ghez et al. (1998) Ghez, A. M., Klein, B. L., Morris, M., & Becklin, E. E. 1998, ApJ, 509, 678, doi: 10.1086/306528
  * Giesers et al. (2018) Giesers, B., Dreizler, S., Husser, T.-O., et al. 2018, Monthly Notices of the Royal Astronomical Society: Letters, 475, L1, doi: 10.1093/mnrasl/slx203
  * Giesers et al. (2019) Giesers, B., Kamann, S., Dreizler, S., et al. 2019, A&A, 632, A3, doi: 10.1051/0004-6361/201936203
  * Górski et al. (2005) Górski, K. M., Hivon, E., Banday, A. J., et al. 2005, ApJ, 622, 759, doi: 10.1086/427976
  * Gravity Collaboration et al. (2018) Gravity Collaboration, Abuter, R., Amorim, A., et al. 2018, A&A, 615, L15, doi: 10.1051/0004-6361/201833718
  * Grondin et al. (2023) Grondin, S. M., Webb, J. J., Leigh, N. W. C., Speagle, J. S., & Khalifeh, R. J. 2023, MNRAS, 518, 4249, doi: 10.1093/mnras/stac3367
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357, doi: 10.1038/s41586-020-2649-2
  * Harris (2010) Harris, W. E. 2010, arXiv e-prints, arXiv:1012.3224.  https://arxiv.org/abs/1012.3224
  * Hattori et al. (2018) Hattori, K., Valluri, M., Bell, E. F., & Roederer, I. U. 2018, ApJ, 866, 121, doi: 10.3847/1538-4357/aadee5
  * Heggie & Hut (2003) Heggie, D., & Hut, P. 2003, The Gravitational Million-Body Problem: A Multidisciplinary Approach to Star Cluster Dynamics
  * Hénon (1971a) Hénon, M. 1971a, Ap&SS, 13, 284, doi: 10.1007/BF00649159
  * Hénon (1971b) Hénon, M. H. 1971b, Ap&SS, 14, 151, doi: 10.1007/BF00649201
  * Hills (1988) Hills, J. G. 1988, Nature, 331, 687, doi: 10.1038/331687a0
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * Hurley et al. (2000) Hurley, J. R., Pols, O. R., & Tout, C. A. 2000, MNRAS, 315, 543, doi: 10.1046/j.1365-8711.2000.03426.x
  * Hurley et al. (2002) Hurley, J. R., Tout, C. A., & Pols, O. R. 2002, MNRAS, 329, 897, doi: 10.1046/j.1365-8711.2002.05038.x
  * Igoshev et al. (2022) Igoshev, A. P., Perets, H., & Hallakoun, N. 2022, MNRAS, doi: 10.1093/mnras/stac3488
  * Irrgang et al. (2021) Irrgang, A., Dimpel, M., Heber, U., & Raddi, R. 2021, A&A, 646, L4, doi: 10.1051/0004-6361/202040178
  * Joshi et al. (2000) Joshi, K. J., Rasio, F. A., & Portegies Zwart, S. 2000, ApJ, 540, 969, doi: 10.1086/309350
  * Koposov et al. (2020) Koposov, S. E., Boubert, D., Li, T. S., et al. 2020, MNRAS, 491, 2465, doi: 10.1093/mnras/stz3081
  * Kremer et al. (2019) Kremer, K., Lu, W., Rodriguez, C. L., Lachat, M., & Rasio, F. A. 2019, ApJ, 881, 75, doi: 10.3847/1538-4357/ab2e0c
  * Kremer et al. (2018) Kremer, K., Ye, C. S., Chatterjee, S., Rodriguez, C. L., & Rasio, F. A. 2018, ApJ, 855, L15, doi: 10.3847/2041-8213/aab26c
  * Kremer et al. (2020) Kremer, K., Ye, C. S., Chatterjee, S., Rodriguez, C. L., & Rasio, F. A. 2020, in Star Clusters: From the Milky Way to the Early Universe, ed. A. Bragaglia, M. Davies, A. Sills, & E. Vesperini, Vol. 351, 357–366, doi: 10.1017/S1743921319007269
  * Leonard (1991) Leonard, P. J. T. 1991, AJ, 101, 562, doi: 10.1086/115704
  * Lynden-Bell & Wood (1968) Lynden-Bell, D., & Wood, R. 1968, MNRAS, 138, 495, doi: 10.1093/mnras/138.4.495
  * Maccarone et al. (2007) Maccarone, T. J., Kundu, A., Zepf, S. E., & Rhode, K. L. 2007, Nature, 445, 183
  * Miller-Jones et al. (2015) Miller-Jones, J. C. A., Strader, J., Heinke, C. O., et al. 2015, Monthly Notices of the Royal Astronomical Society, 453, 3919, doi: 10.1093/mnras/stv1869
  * Morscher et al. (2013) Morscher, M., Umbreit, S., Farr, W. M., & Rasio, F. A. 2013, ApJ, 763, L15, doi: 10.1088/2041-8205/763/1/L15
  * Oey et al. (2023) Oey, M. S., Dorigo Jones, J., Phillips, G. D., et al. 2023, arXiv e-prints, arXiv:2301.11444, doi: 10.48550/arXiv.2301.11444
  * pandas development team (2020) pandas development team, T. 2020, pandas-dev/pandas: Pandas, latest, Zenodo, doi: 10.5281/zenodo.3509134
  * Pattabiraman et al. (2013) Pattabiraman, B., Umbreit, S., Liao, W.-k., et al. 2013, ApJS, 204, 15, doi: 10.1088/0067-0049/204/2/15
  * Poveda et al. (1967) Poveda, A., Ruiz, J., & Allen, C. 1967, Boletin de los Observatorios Tonantzintla y Tacubaya, 4, 86
  * Renzo et al. (2019) Renzo, M., Zapartas, E., de Mink, S. E., et al. 2019, A&A, 624, A66, doi: 10.1051/0004-6361/201833297
  * Rodriguez et al. (2022a) Rodriguez, C. L., Hafen, Z., Grudić, M. Y., et al. 2022a, arXiv e-prints, arXiv:2203.16547, doi: 10.48550/arXiv.2203.16547
  * Rodriguez et al. (2022b) Rodriguez, C. L., Weatherford, N. C., Coughlin, S. C., et al. 2022b, ApJS, 258, 22, doi: 10.3847/1538-4365/ac2edf
  * Rui et al. (2021) Rui, N. Z., Kremer, K., Weatherford, N. C., et al. 2021, ApJ, 912, 102, doi: 10.3847/1538-4357/abed49
  * Spitzer (1987) Spitzer, L. 1987, Dynamical evolution of globular clusters (Princeton University Press)
  * Strader et al. (2012) Strader, J., Chomiuk, L., Maccarone, T. J., Miller-Jones, J. C. A., & Seth, A. C. 2012, Nature, 490, 71, doi: 10.1038/nature11490
  * Tauris (2015) Tauris, T. M. 2015, MNRAS, 448, L6, doi: 10.1093/mnrasl/slu189
  * VandenBerg et al. (2013) VandenBerg, D. A., Brogaard, K., Leaman, R., & Casagrande, L. 2013, ApJ, 775, 134, doi: 10.1088/0004-637X/775/2/134
  * Weatherford et al. (2022) Weatherford, N. C., Kıroğlu, F., Fragione, G., et al. 2022, arXiv e-prints, arXiv:2211.16523, doi: 10.48550/arXiv.2211.16523
  * Wes McKinney (2010) Wes McKinney. 2010, in Proceedings of the 9th Python in Science Conference, ed. Stéfan van der Walt & Jarrod Millman, 56 – 61, doi: 10.25080/Majora-92bf1922-00a