Stellar Escape from Globular Clusters. II. Clusters May Eat Their Own Tails

We apply for the first time orbit-averaged Monte Carlo star cluster simulations to study tidal tail and stellar stream formation from globular clusters (GCs), assuming a circular orbit in a time-independent spherical Galactic potential. Treating energetically unbound bodies—potential escapers (PEs)—as collisionless enables this fast but spherically symmetric method to capture asymmetric extratidal phenomena with exquisite detail. Reproducing stream features such as epicyclic overdensities, we show how returning tidal tails can form after the stream fully circumnavigates the Galaxy, enhancing the stream's velocity dispersion by several kilometers per second in our ideal case. While a truly clumpy, asymmetric, and evolving Galactic potential would greatly diffuse such tails, they warrant scrutiny as potentially excellent constraints on the Galaxy’s history and substructure. Reexamining the escape timescale Δt of PEs, we find new behavior related to chaotic scattering in the three-body problem; the Δt distribution features sharp plateaus corresponding to distinct locally smooth patches of the chaotic saddle separating the phase-space basins of escape. We study for the first time Δt in an evolving cluster, finding that Δt∼(EJ−0.1,EJ−0.4) for PEs with (low, high) Jacobi energy E J, flatter than for a static cluster ( EJ−2 ). Accounting for cluster mass loss and internal evolution lowers the median Δt from ∼10 Gyr to ≲100 Myr. We finally outline potential improvements to escape in the Monte Carlo method intended to enable the first large grids of tidal tail/stellar stream models from full GC simulations and detailed comparison to stream observations.


INTRODUCTION
Recent astronomical surveys, especially the exquisite kinematics from Gaia (Gaia Collaboration et al. 2016), have revealed extensive substructure in the Milky Way (MW), including nearly 100 stellar streams-for a recent review and catalog, see Helmi (2020) and Mateu (2023), respectively.Many of these strands of stars sharing similar orbits are debris from either dwarf galaxies recently accreted and tidally disrupted by the MW or dissolved/dissolving globular clusters (GCs; e.g., Bonaca et al. 2021).Aiding this linkage, tidal tails directly emanate from ≈15 MWGCs (e.g., Piatti & Carballo-Bello 2020), some extending into full streams.
As kinematically cold structures sensitive to subtle perturbations on long timescales, stellar streams are outstanding probes of the MW's mass and gravitational potential, dark matter (DM) halo, internal substructure, and assembly history (e.g., Koposov et al. 2010;Bonaca et al. 2014;Bonaca & Hogg 2018;Küpper et al. 2015;Bovy et al. 2016).Stream morphology and kinematics can provide especially tantalizing hints to the nature of particulate DM (e.g., Banik et al. 2021); for example, cold DM should produce cuspier (more centrally dense) DM subhalos than alternatives such as warm DM (Peebles 1982), which suppresses primordial density fluctuations on small scales, or ultra-light scalar 'fuzzy' DM (Hui et al. 2017), which forms less small-scale substructure.So the frequency and morphology of gaps, and kinematic heating, induced in streams by passing DM subhalos depends on DM's identity (e.g., Ibata et al. 2002;Johnston et al. 2002;Carlberg 2009).
Yet, capitalizing on the prospect of DM inference with stellar streams has proven challenging since many phenomena simultaneously and degenerately affect stream morphology and kinematics.Stream gaps, overdensities, and fanning naturally arise via the epicyclic trajectories of slow GC escapers even in static MW potentials without any clumpy substructure (e.g., Capuzzo Dolcetta et al. 2005;Küpper et al. 2008), and also from perturbations by many baryonic substructures, including giant molecular clouds (Amorisco et al. 2016), the MW spiral arms or disk (Banik & Bovy 2019;Carlberg & Agler 2023), infalling MW satellites (Garavito-Camargo et al. 2019), and other GCs (Doke & Hattori 2022).Further complicating matters, MW substructure such as a rotating bar (Hattori et al. 2016;Pearson et al. 2017) or even modified gravity (Thomas et al. 2018;Kroupa et al. 2022) can induce stream asymmetry and observational artifacts can cause false gaps/overdensities (Ibata et al. 2020).Internal GC dynamics is impactful, too; in particular, dynamical heating from central black hole (BH) binaries can dramatically accelerate the evaporation rate at the end of a GC's life (e.g., Banerjee & Kroupa 2011;Whitehead et al. 2013;Contenta et al. 2015;Chatterjee et al. 2017;Giersz et al. 2019;Weatherford et al. 2021), thereby increasing the tidal tail density (Gieles et al. 2021;Gieles & Gnedin 2023;Roberts et al. 2024).Disentangling the individual contributions of such diverse phenomena in observations (e.g., to measure the properties of DM subhalos specifically) will require comprehensive modeling of stream formation and disruption.
The importance of stellar streams to Galactic archeology strongly motivates further theoretical study of escape from GCs.This work is the second of a series on this topic that began with a thorough exploration of escape mechanisms (Weatherford et al. 2023, hereafter W23).W23 emphasized more strongly high-speed ejection, relevant to production of runaway or hypervelocity stars in the MW halo.Focusing more on low-speed escape, we now apply for the first time orbit-averaged GC models-using the Hénon (1971a,b) Monte Carlo method-to study stellar stream formation.
At first glance, orbit-averaged GC models are not an obvious choice for this application.The impacts of the GC potential and collisional dynamics steeply diminish beyond the tidal boundary, so it is common when simulating streams to neglect the GC's internal dynamics in favor of collisionless 'streakline' or 'particle spray' techniques (e.g., Varghese et al. 2011;Küpper et al. 2012;Lane et al. 2012;Bonaca et al. 2014;Gibbons et al. 2014;Fardal et al. 2015;Shipp et al. 2018;Grondin et al. 2023).Stars are simply sprayed at a prescribed rate and velocity distribution from a point representing the GC (or the Lagrange points on its tidal boundary) orbiting in an external Galactic potential.Key benefits of fully modeling the GC are then to accurately determine the stellar escape rates and velocities, and internal stream details such as stellar masses, luminosities, and binary properties.Direct summation N-body codes (e.g., nbody6; Aarseth 1999) are the "gold standard" for such purposes due to their exquisite accuracy and accommodation of GC asymmetry (essential for tidal physics) and stages of the GC's life spent out of virial equilibrium or with small N (essential for late-stage dissolution).Yet direct summation is computationally intensive; despite steady computing advances the method still has not been used to model a GC with a density and initial N typical of MWGCs over a full Hubble time (e.g., Wang et al. 2016;Arca Sedda et al. 2023).So direct summation is illsuited to generating large model grids of extratidal structures from dense GCs, necessary to optimally match specific MW streams or explore the many aforementioned factors that can significantly influence stream properties.
The Monte Carlo method is a much faster alternative to direct summation, capable of simulating MWGCs of typical mass and density using 10 3 -10 4 -times fewer CPU-hours.Of the two main algorithmic lineages (excellently summarized by Vasiliev 2015), the orbit-averaged version descended from Hénon (1971a,b) is the only one still in regular usespecifically as the codes MOCCA and our CMC (respectively overviewed by Giersz et al. 2013 andRodriguez et al. 2022).Instead of direct orbit integration, Hénon's method applies a statistical treatment of stellar dynamics by modeling the cumulative effect of many distant two-body encounters as a single effective scattering between neighboring bodies, radially sorted into pairs.This occurs each timestep on the relaxation timescale t r rather than the dynamical timescale, leveraging that distant two-body encounters dominate evolution in each body's energy E and angular momentum J.Only especially strong encounters are handled via a small-N direct integrator, greatly hastening the computation.At the end of the timestep, each body's position and velocity (r, v r , v t ) are randomly sampled from the time-averaged orbit consistent with the body's new E and J.
Hénon's method compares well to both direct summation and observations of typical GCs in most regimes (e.g., Giersz et al. 2013;Rodriguez et al. 2016;Kremer et al. 2020), but imposes three key assumptions: (1) the radial sorting and orbit averaging impose spherical symmetry (but not velocity isotropy), (2) the statistical treatment imposes a large N and virial equilibrium, preventing study of the final stages of GC dissolution, and (3) the orbit averaging neglects global evolution occurring on timescales ≪ t r , such as sharply evolving tides (tidal shocks).All three assumptions are relevant to escape, but most critically the assumed spherical symmetry means that unmodified the method cannot resolve tidal tails.
We test a simple solution to this challenge by removing bodies from CMC once they obtain sufficient energy to escape (usually in the GC's core) and evolving their trajectories from there in the full asymmetric tidal field.We explore how well this approach-implemented as a post-processing step applicable to any CMC simulation-reproduces known features of tidal tails and stellar streams in the case of a circular GC orbit in an unevolving, spherical MW potential.While a vast simplification, this well-studied regime-standard for Hénon's method-is ideal for a first test of our approach.Crucially, it also enables us to investigate several nuances that have yet to be explored in even this simple case, including chaotic behavior in the escape (survival) timescale of energetically unbound bodies within GCs, modification of this timescale by evolution of the GC potential, and impact of return trajectories.Note our modification does not generalize the internal dynamics of the Monte Carlo method to non-spherical potentials (e.g., Vasiliev 2015) or its escape prescriptions to evolving external tides (e.g., Sollima & Mastrobuono Battisti 2014;Rodriguez et al. 2023).Such enhancements are complementary to our approach and the latter especially are similar to our intended future upgrades to CMC (see Section 6).
As in W23, we shall often refer to escape before and after core collapse, the GC's observable change from a flat (noncore-collapsed; NCC'd) to a steep (core-collapsed; CC'd) central surface brightness.This occurs upon ejection of the BH population, which weakens binary burning-hardening of binaries via encounters with passing bodies (e.g., Heggie 1975;Hills 1975).The potential energy released by the binaries heats the bodies involved, supporting the GC's core against collapse.BH binaries are strong heat sources due to their mass, but GCs born especially dense quickly harden and eject them; binary burning then relies on (less massive) white dwarfs, reducing heating and allowing the core to observably collapse (e.g., Chatterjee et al. 2013;Kremer et al. 2019Kremer et al. , 2020Kremer et al. , 2021;;Rui et al. 2021a).We recap how this affects GC escape mechanisms in Section 5.1, but see also W23.
The paper is organized as follows.We review relevant theory in Section 2 and describe our GC simulations in Section 3. We explain how we integrate (in post-processing) the trajectories of unbound bodies removed from the simulations in Section 4. In Section 5, we first examine the escape energies and escape timescale from typical NCC'd and CC'd GCs.We then examine how GC mass loss affects the escape timescale and how return trajectories impact stream morphology and kinematics.In Section 6, we discuss observational implications, additional complexities and caveats, and potential upgrades to escape in CMC.We conclude with a summary of our findings and planned future work in Section 7.

THEORY
We examine escape in the simplified case of an evolving spherical GC potential ϕ c circularly orbiting within an unevolving spherical Galactic potential ϕ g .We assume that once a cluster member becomes a potential escaper (PE) by gaining enough energy to eventually escape it ceases to be influenced by scattering with other stellar bodies.This collisionless approximation allows the PE mass distribution to follow the non-spherical tidal field despite the spherical ϕ c in CMC.The PE trajectories in this case are an example of the circular restricted three-body problem (CR3BP, e.g., Szebehely 1967;Marchal 1990;Hénon 1997;Valtonen & Karttunen 2006;Koon et al. 2011).In this Section, we review the elements of the CR3BP most essential to our new results on the escape timescale and return trajectories, leaving details of our algorithmic implementation to Section 3.While the classic CR3BP features unevolving ϕ c and ϕ g , much of its conceptual framework remains useful in our case, and we comment where relevant on the impact of ϕ c 's time-dependence.

Problem Setup
We study trajectories of stellar bodies in the gravitational fields of their birth GC-of slowly evaporating mass m(t)and the MW, each on instantaneously circular orbits about their mutual center-of-mass at a constant angular speed ω = v g /R gc , where v g is the GC's circular speed at Galactocentric distance R gc .To satisfy Kepler's third law at t = 0, we define the MW's mass enclosed within R gc to be a constant . Technically, a constant M, ω, v g , and R gc is mildly inconsistent with a decreasing m(t) that would transfer about half of its lost mass to the volume within R gc , but these assumptions are very accurate since m/M ∼ 10 −6 for a MWGC and since the circular speed profile in the MW halo is very flat.Defining coordinates in the center-of-mass frame XY Z in units of R gc , the MW and GC centers instantaneously trace circles in the XY -plane with respective dimensionless radii [µ(t), 1 − µ(t)], where µ(t) ≡ m(t)/[M + m(t)]. 1e also define a clustercentric frame xyz, in units of R gc , that co-rotates with the GC as it orbits the center-of-mass.In this frame, the GC's center is the origin (r c = 0) with velocity v c = [1−µ(t)]v g ŷ while the MW's center has constant position r g = −x.The effective potential at position r in this frame is where the last term is centrifugal.ϕ eff sets the tidal boundary's shape and the minimum energy a body in the GC must have to escape.If ϕ eff were static and the cluster dynamics collisionless, then the motion of each body with speed v in the xyz frame would conserve a quantity known as the Jacobi energy (in our case time-dependent): Here A(t) = max [ϕ eff (z = 0,t)], achieved at the L4/L5 Euler-Lagrange points, tips of both in-plane equilateral triangles whose base vertices are r c and r g .The convention of subtracting A (e.g., Spitzer 1987;Fukushige & Heggie 2000) has no impact on trajectories but affects normalized energy definitions like Ẽ in Equation (3), so is essential when later comparing to the latter study, hereafter referred to as FH00.Since v2 ≥ 0, Equation (2) implies that bodies with energy E J cannot enter regions where ϕ eff (r) > E J + A. Figure 1 shows the xy-plane's intersection with these 'forbidden realms' (gray) for several E J (in its normalized form Ẽ; see below).As E J increases, the zero-velocity surface bounding these realms (ϕ eff = E J + A) ranges from separately enclosing the GC and MW centers (upper left/middle panels), to allowing passage through one or two openings directly toward and away from the latter (upper right-lower center) to disappearing entirely from the xy-plane (lower right).
Slowly shrinking due to global mass loss, the GC's tidal boundary lies on the highest-E J zero-velocity surface that (at the time) fully encloses the GC (upper center).This boundary terminates nearest to the Galactic center at a saddle point of ϕ eff (r) known as the L1 Euler-Lagrange point and furthest just shy of another saddle, L2. 2 Their respective locations, r 1 (t) and r 2 (t), are numerically solvable from the definition of a critical point, ∇ϕ eff (r,t) = 0, but for MWGCs (µ ∼ 10 −6 ), a 2nd-order expansion about µ → 0 is very accurate-to about one part in 4µ −1/3 (Szebehely 1967).In this limit |r 2 − r c | → |r 1 − r c | and ϕ eff (r 2 ) → ϕ eff (r 1 ), so the tidal boundary is symmetric, terminating exactly at both saddles.Yet for generality we define the tidal radius as the maximum clustercentric distance to the tidal boundary, r t ≡ |r 1 − r c |R gc .For Keplerian ϕ c and ϕ g , the expansion under µ → 0 yields r t /R gc ≈ (µ/3) 1/3 , becoming r t /R gc ≈ (µ/2) 1/3 for a logarithmic ϕ g closer to that of the MW halo (e.g., Spitzer 1987).
Bodies within the tidal boundary can only pass beyond it if they have E J (t) > E crit (t), where E crit (t) ≡ ϕ eff (r 1 ,t) − A(t).We shall refer to this escape criterion as the raw criterion  (Li) when v = 0.These points and the centers of the Galaxy (G) and cluster (C) are labeled in the first panel but appear in all.The forbidden realm (gray) is bounded by the zero-velocity surface and cannot be entered by bodies with the given Ẽ or lower.The blue curves in the last four panels exemplify trajectories at the given Ẽ initiated with (x, y, z) = (µ, 0, 0), velocity v0 ∥ ŷ, and integrated for 20 full orbits of the cluster in the Galaxy.Upper left: Ẽ < 0, so bodies cannot transit between the cluster interior, Galactic interior (within Rgc), or Galactic exterior (beyond Rgc).Upper center: Ẽ = 0, so the zero-velocity surface can expand no further before opening a neck between the cluster and Galactic interiors at L1.The closed portion of the surface passing through L1 is the tidal boundary.Upper right: Ẽ1 < Ẽ < Ẽ2 opens the neck at L1.The example trajectory-nearly periodic in the rotating frame-transits back-and-forth through the cluster/Galactic interiors and illustrates an escaper taking many tc to escape before eventually returning to the cluster.Lower left: Ẽ2 < Ẽ < Ẽ3 opens a second neck at L2, connecting all three domains.The example trajectory illustrates an escaper that immediately finds the neck to the Galactic interior and briefly passes again through the cluster after nearly 20 full Galactocentric orbits.Lower center: Ẽ3 < Ẽ < Ẽ4 opens a third neck at L3.The example trajectory transits between all three domains, with the escaper temporarily returning to the cluster after only four orbits about the Galaxy (to the cluster's three in that time).Lower right: Ẽ = Ẽ4 = 1, at which point the entire xy-plane is available to escapers, though the zero-velocity surface still exists out-of-plane.The example trajectory alternates between orbiting in the Galactic interior/exterior, passing back through the cluster during each transit.Though at higher-Ẽ and not nearly as regular, this trajectory is loosely analogous to those of some Jovian comets.
to distinguish it from more complex modeling alternatives discussed later.Following FH00 and to emphasize the trajectories' extreme sensitivity to E J as E J → E + crit , we henceforth express this criterion in normalized form as This definition has the benefit that Ẽ = Ẽ4 = 1 when ϕ eff = A.
Here and elsewhere, the Ẽ subscript i indicates it is the value of Ẽ at the Euler-Lagrange point Li for v = 0.The space accessible to PEs (white space in Figure 1) varies with Ẽ and reduces to three domains: the cluster, and the Galactic interior and exterior (beyond the tidal boundary at Galactocentric distances <R gc and >R gc , respectively).
For Ẽ ≤ 0, transit between domains is energetically disallowed, but once Ẽ > Ẽ1 = 0, a neck in the zero-velocity surface opens at L1 to allow transit between the cluster/Galactic interior.A similar neck at L2 allows transit between the cluster/Galactic exterior once Ẽ > Ẽ2 .This is the most relevant geometry to tidal tail formation from MWGCs, for which Ẽ1 ≈ Ẽ2 .Only once Ẽ > Ẽ3 (corresponding to L3, the third and final saddle point of ϕ eff ) is direct transit between the Galactic interior/exterior possible.The forbidden realm disappears entirely from the xy-plane once Ẽ ≥ Ẽ4 = 1, though it still exists at z ̸ = 0, receding away from the xy-plane as Ẽ grows.So escape at high Ẽ is nearly unconstrained while escapers at low Ẽ must pass near L1/L2.

Tidal Tail and Stellar Stream Formation
The equations of motion in the CR3BP in the xyz frame are The second term in each equality, the Coriolis acceleration, drives much of the trajectories' behavior, including production of tidal tails from low-Ẽ PEs.Since these must escape near L1/L2, their velocities in the necks are biased to be parallel to ∓x.The Coriolis effect then bends these trajectories into tails leading/trailing the GC, respectively, and induces epicycles in each trajectory's projection onto the xyplane (Figure 1).The epicycles' characteristic size depends on µ and Ẽ, but for (µ, Ẽ) ≪ 1, the Coriolis effect keeps escapers near to the GC's orbital path, forming elongated stellar streams along it.Since each epicycle has a turning point minimizing the velocity parallel to the stream, there are periodic over-densities in the streams near these points, spaced ∼10r t apart for Ẽ ≪ 1 (e.g., Capuzzo Dolcetta et al. 2005;Küpper et al. 2008Küpper et al. , 2010Küpper et al. , 2012;;Just et al. 2009).The Coriolis effect also stabilizes trajectories retrograde to the GC orbit, slowing escape of bodies on retrograde orbits and allowing those with Ẽ > 0 within and even beyond the tidal boundary to remain near the GC for an arbitrarily long time (Hénon 1970;FH00;Read et al. 2006;Ernst et al. 2008).

Dynamical Systems Theory and Return Trajectories
Stellar escape from GCs can be given extensive mathematical formalism from a standpoint of dynamical systems theory (e.g., FH00; Ernst et al. 2008;Tanikawa & Fukushige 2010;de Assis & Terra 2014;Zotos 2015a,b;Zotos & Jung 2017).For each neck, twin, infinitely winding/branching tubes in phase space known as invariant manifolds enclose all possible transit trajectories into or out of the GC.Trajectories outside the invariant manifolds cannot transit between the CR3BP's three domains, implying that some portion of even bodies with Ẽ > 0 would never escape without the aid of perturbations from other bodies or evolution in ϕ eff (e.g., via GC mass loss, an evolving ϕ g , an eccentric GC orbit, or passage near MW substructure).These effects induce additional phase space diffusion, enabling bodies otherwise stuck forever on non-transiting orbits outside the invariant manifolds to move to transiting orbits within them.Reframed energetically, the internal gravitational scattering and evolution of ϕ eff cause diffusion in Ẽ.Since the shapes of the manifolds change with Ẽ, this effectively smears them out to encompass more phase space, promoting escape.
Trajectory families known as heteroclinic orbits asymptotically connect the saddle points to each other.Through the union of such orbits, one can design heteroclinic chains that, e.g., periodically alternate between orbiting in the interior and exterior Galactic domains, crossing (and temporarily orbiting within) the GC during each transit.The intentional design of such itineraries in the solar system context is essential to many space missions (e.g., Koon et al. 2011) but these trajectories exist naturally, too.For example, some comets, including Oterma and Gehrels 3 (e.g., Belbruno & Marsden 1997;Koon et al. 2000), periodically transition from heliocentric orbits outside to inside the orbit of Jupiter and vice versa.During each transition, the comets become temporarily bound to Jupiter.The existence of such orbits (with the Sun-Jupiter system analagous to the MW-GC system) means that GCs can re-capture past members that previously escaped.In GC models, the recapture rate is typically assumed negligible due to a time-dependent and/or asymmetric ϕ g and diffusion of the return trajectories by MW substructure.Yet the rate may be significant in the case of a nearly circular GC orbit in an unevolving ϕ g .Indeed, we demonstrate that robust 'returning tidal tails' resulting from such trajectories can form in this case in Section 5.4.Though idealized, this case is a reasonable starting point since the impact of returning bodies on stellar streams has not been studied before.The idealization also enables us to determine upper bounds on the impact of returners on stream observables (e.g., surface density and velocity dispersion) under more realistic conditions.

Escape Timescale
Weak, diffusive two-body relaxation has long been known to dominate energy transport in GCs (e.g., Ambartsumian 1937;Spitzer 1940;Chandrasekhar 1942).Despite some early caveats in the idealized context of isolated GCs (e.g., Hénon 1960Hénon , 1969)), relaxation therefore naturally dominates escape from GCs in a tidal field (e.g., Spitzer & Shapiro 1972;W23).So while various strong encounters propel some bodies to Ẽ ≫ 1, ejecting them on the crossing timescale t c , most bodies first acquire sufficient energy to escape by gradually random-walking in energy to Ẽ > 0, remaining at Ẽ ≪ 1 for potentially many relaxation times.
Since the necks about L1/L2 are narrow for Ẽ ≪ 1, relaxation-driven PEs may take many t c to 'find' and escape through these necks, a fact highlighted in the seminal work by FH00 (see also Tanikawa & Fukushige 2010).Via direct N-body modeling with analytical support, they showed that the typical escape timescale ∆t esc -between first satisfying Ẽ > 0 and crossing beyond r t -is of order an entire Hubble time for MWGCs in the idealized case of a static ϕ c .They also found that for Ẽ ≪ 1 and static ϕ c , then ∆t esc ∝ Ẽ−2 .Yet ∆t esc has attracted little further study, likely due to its limited relevance to direct N-body models (though see Baumgardt 2001), which typically only remove bodies once they pass beyond several r t anyway.∆t esc is important for models that use energy-based escape criteria, however, like most Fokker-Planck or Monte Carlo codes, including CMC and MOCCA.For reasons we discuss in the next Section, such codes usually remove bodies immediately once they satsify Ẽ > 0 or a similar energy-based criterion, effectively assuming ∆t esc = 0 (e.g., Spitzer & Shull 1975;Lee & Goodman 1995;Giersz et al. 2008;Chatterjee et al. 2010).Especially since the MOCCA code now implements delayed escape based on FH00's idealized findings (Giersz et al. 2013), ∆t esc 's extreme length merits revisiting.In particular, ∆t esc in a realistically evolving ϕ c -a highly practical scenario that should hasten escape due to GC mass loss-has yet to be explored at all, motivating our later attention to this case.

GLOBULAR CLUSTER SIMULATIONS
We employ the latest public version of the Cluster Monte Carlo code (CMC; Rodriguez et al. 2022) to simulate GCs.CMC includes numerous physical processes essential to GC evolution, including internal stellar evolution with COSMIC (Breivik et al. 2020b)-an updated version of SSE/BSE (Hurley et al. 2000(Hurley et al. , 2002))-two-body relaxation (Joshi et al. 2000;Pattabiraman et al. 2013), Galactic tides (Joshi et al. 2001;Chatterjee et al. 2010), strong binary encounters and physical collisions (Fregeau et al. 2003;Fregeau & Rasio 2007), and three-body binary formation (Morscher et al. 2013(Morscher et al. , 2015)).CMC simulates strong binary encounters via the small N-body direct integrator fewbody, which accounts for post-Newtonian dynamics (Fregeau et al. 2004;Antognini et al. 2014;Amaro-Seoane & Chen 2016;Rodriguez et al. 2016Rodriguez et al. , 2018a,b),b).CMC also allows for two-body binary formation through gravitational wave dissipation and tidal capture (Kremer et al. 2021;Ye et al. 2022), but for simplicity we only allow the former in this study.We study escape in four cases (Table 1): an archetypal noncore-collapsed (NCC'd) and core-collapsed (CC'd) MWGC, each under two distinct escape criteria.The corresponding GCs from the CMC Cluster Catalog (Kremer et al. 2020) are numbered 2 and 8 in W23.Other than the escape criteria (see below) these simulations differ only in the initial virial radius r v = 0.5 pc (for the CC'd GC) and 2 pc (for the NCC'd GC).All other initial parameters are identical, such as the initial total number of particles (singles plus binaries) N i = 8 × 10 5 , Galactocentric distance R gc = 8 kpc, and metallicity Z/Z ⊙ = 0.1.As in W23, stellar masses (primary mass, in the case of a binary) draw from the standard Kroupa (2001) stellar initial mass function from 0.08-150 M ⊙ .We also assume that all neutron stars born in core-collapse (electron-capture) supernovae receive natal kicks drawn from a Maxwellian with dispersion σ = 265 km s −1 (20 km s −1 ).Natal kicks for BHs share the core-collapse kick distribution of the neutron stars, but reduced in magnitude by the fraction of the stellar envelope's mass that falls back onto the BH via the prescription from Hobbs et al. (2005)-see also Fryer et al. (2012).For any further details or justifications regarding the models other than the two updates from W23 mentioned below, see W23 and references therein.
For completeness, we note two minor differences from the earlier versions of the models used in W23.Both changes are purely to bolster consistency with concurrent and future works using CMC.First, instead of a uniform initial binary fraction f b = 5%, we now set f b = 5% for stars born with mass M < 15 M ⊙ (in accord with present-day MWGC observations, e.g., Milone et al. 2012) and f b = f b,high = 50% for those born with M ≥ 15 M ⊙ (more in line with observations in young star clusters, e.g., Sana et al. 2009Sana et al. , 2012;;Moe & Di Stefano 2017).Though still simplistic, this two-component initial binary fraction is more consistent with models of binary formation in star clusters embedded in molecular clouds (e.g., Cournoyer-Cloutier et al. 2021;Guszejnov et al. 2023), in which f b increases steeply with stellar mass.The update has minimal impact on the GC evaporation rate or distribution of escape mechanisms and energies (evident by comparing later figures to those in W23).We also now utilize the Fryer et al. (2012) 'delayed' supernova prescription-CMC's new default-to compute compact remnant masses since this treatment produces a robust remnant population in the BH lower-mass gap (2-5 M ⊙ ), in line with observations from the third Gravitational Wave Transient Catalog (Abbott et al. 2023).This was not the case under the Fryer et al. (2012) 'rapid' supernova prescription used in the CMC Cluster Catalog models used in W23.This change has negligible dynamical impact but expands the population of bodies labeled as BHs and so explains an increase in BH ejections seen in later figures relative to those in W23.

Escape Criteria in CMC
While CMC allows for an arbitrary time-varying tidal field (specified by a tidal tensor; see Rodriguez et al. 2023), we explore CMC's default tidal scenario-the time-dependent CR3BP from Section 2 with a logarithmic Galactic potential and ωR gc = v g = 220 km s −1 , typical of MWGCs (e.g., Spitzer 1987;Binney & Tremaine 2008).CMC's tidal radius is then Note the time-dependence due to GC mass-loss via escape and stellar evolution.The GC's spherical potential ϕ CMC c (r,t) is similarly time-dependent, computed at each CMC timestep from the positions of all bodies.Unless otherwise noted, this time-dependence carries over to our integration of escape trajectories in post-processing.Its main impact is to slowly raise (make less negative) ϕ c and thereby Ẽ, promoting escape.
Since the raw escape criterion Ẽ > 0 involves ϕ eff , CMC's spherical version of this criterion is approximate, becoming where and As explained in W23's Appendix A.2, a term ≈ − (r/r t ) 2 /B (where B = 12 for logarithmic ϕ g or 9 for Keplerian ϕ g ) is omitted from the coefficient of Equation ( 8) since most bodies first satisfy the escape criterion with r/r t ≪ 1.The raw energy criterion for escape has never strictly been used in CMC.Instead, CMC's current default is to remove bodies in a two-step process at every timestep: first any bodies with clustercentric apocenter distance r a > r CMC t , and then any remaining bodies satisfying a modified version of the raw criterion we call the α criterion (Giersz et al. 2008): where In the Coulomb logarithm, ln Λ = ln (γN i ), we use γ = 0.01, appropriate for GCs with realistic stellar initial mass functions (e.g., Freitag et al. 2006;Rodriguez et al. 2018cRodriguez et al. , 2022)).The first step's purpose is simply to ease comparison to GC models using apocenter-based criteria (e.g., CMC before 2010 or most direct N-body codes).3Physically, the two steps reduce to solely the α criterion, which all bodies with r a > r CMC t must satisfy for α > 1 (N i ≳ 5 × 10 3 , below which Hénon's method is unreliable anyway; Aarseth et al. 2008).
The α criterion's purpose is to make escape slightly harder to account for scattering of PEs back to Ẽ < 0 (e.g., Chandrasekhar 1942;King 1959;Baumgardt 2001).This would occur if the PEs were to continue participating in collisional dynamics on their way out of the GC, since weak encounters induce a random walk in Ẽ.Over the escape timescale (typically long; FH00), this could either increase or decrease Ẽ-in the latter case, potentially delaying escape.The effect's strength scales inversely with N i ; specifically, Baumgardt (2001) found that averaged over the half-mass time of a GC, the fraction of bodies within r t that are PEs scales roughly as ∼[ln(Λ)/N i ] 1/4 .Using this factor as a measure of the strength of back-scattering, Giersz et al. (2008) tuned the second coefficient in Equation ( 10) so that Hénon-type codes like MOCCA and CMC best match the evaporation rate from direct N-body codes.Since the only intent was to improve accuracy for internal GC evolution, the impact on tidal tails (e.g., their velocity dispersion and corresponding width) was not considered and remains untested.
It is unclear if the α criterion should lead to more accurate tidal tails than the raw criterion.By roughly accounting for the impact of ongoing collisional dynamics on PEs, the former may better reproduce the escape timescale for PEs, and hence the evaporation rate, as it slightly raises the Ẽ at which removal from CMC (promotion to 'PE status') occurs.This should hasten escape of PEs, qualitatively matching the expectation if CMC were to allow continued relaxation for PEs rather than removing them immediately.But by raising Ẽ, the α criterion may artificially increase the velocity dispersion (width) in the tidal tails.We therefore study and compare escape under both criteria.Carefully doing so now may aid future efforts to improve escape physics in CMC by providing a rough sense of how much accounting for the impact of collisional dynamics on PEs changes the morphology and kinematics of tidal tails.

ESCAPE TRAJECTORY INTEGRATION
We now describe how we evolve PEs removed from CMC with the galactic dynamics code Gala (Price-Whelan 2017).All of these steps are currently executed in post-processing of CMC output using Python scripts; the functionality is not yet embedded into CMC's base code, though see Section 6.4 for a discussion of such potential enhancements in the future.

Coordinate Systems and Trajectory Initialization
Since CMC imposes spherical symmetry, the output phase space coordinates of removed bodies are (r rmv , v r,rmv , v t,rmv ), the radial position, and radial and tangential velocities, respectively-all in the rotating clustercentric frame xyz at removal time t rmv .(To tidy notation, we forgo subscript 'rmv' in the rest of this subsection.)Yet trajectories in the true non-spherical ϕ eff require full phase space coordinates, so we isotropically project the PEs' positions/velocities at removal into full 6-dimensional phase space.We first orient each PE in the xyz frame with position r = (0, 0, r) and velocity v = (v t , 0, v r ).To distribute v t isotropically relative to v r , we rotate v about ẑ by angle ψ = U(0, 2π), where U(a, b) indicates a random sample from the uniform distribution between a and b.To isotropically distribute r and v, we then rotate each about ŷ by angle θ = arccos[U(−1, 1)] and again about ẑ by angle ϕ = U(0, 2π), yielding positions [x, y, z] = r[sin θ cos ϕ, sin θ sin ϕ, cos θ] and velocities Whenever we convert to the non-rotating XY Z coordinates we also need to specify a phase, so we define the cluster to be located at (X,Y, Z) = (1 − µ, 0, 0) at time t = 0.The coordinate transformation from xyz to XY Z at any time is then

Trajectory Integration
We use the galactic dynamics code Gala (Price-Whelan 2017) to integrate the PE trajectories in the true ϕ eff from the time and location of removal from CMC (usually deep within the GC; W23).The trajectory integration takes place in the rotating frame and uses Gala's wrapper for the DOP853 SciPy integrator, an implementation of the Dormand-Prince method-within the Runge-Kutta family-of order 8(5,3).This implementation internally yields positional errors less than one part in 10 10 per output interval, which we set to be every Myr.When solving the trajectories using the truly time-dependent ϕ c , we do so up to our simulations' final age 14 Gyr (i.e., an integration time of 14 Gyr − t rmv ) or until the first time in the above output that the PE's Galactocentric distance exceeds 2R gc , whichever comes first.This limits the computational burden while enabling resolution of full stellar streams and return trajectories.When solving the trajectories using a constant ϕ c (t rmv ) to compare directly to FH00, we use a full 14 Gyr integration time, regardless of t rmv .

Ensuring Consistency in Ẽ Between CMC and Gala
. To smoothly transition PEs removed from CMC into an orbit integration in Gala, we must ensure consistency between the sphericalized Ẽ estimated by CMC-Equation (6)-and Ẽ numerically computed in Gala from ϕ c , ϕ g , and the PE's position and velocity at t rmv .In Gala, we use the logarithmic Galactic potential ϕ g (R) = v 2 g ln(R) with v g = 220 km s −1 .This is identical to the ϕ g assumed by CMC when computing r CMC t and thereby E CMC crit -Equations ( 5) and ( 8), respectively.It thus ensures that ϕ g 's contribution to Ẽ via E crit is consistent between CMC and Gala.Consistency in the contribution to Ẽ from the cluster potential ϕ c is harder to ensure.While a recent Gala upgrade allows users to input cylindrical spline potentials built from CMC-like lists of masses and positions, this functionality is still being optimized; at present, the computation is orders of magnitude faster for analytic potentials.
An analytic fit to ϕ CMC c must be chosen carefully.While ϕ g 's functional form has a greater impact on the tidal boundary,4 ϕ c more strongly affects Ẽ since PEs typically first satisfy the escape criterion in the GC's core (W23).So even slight inconsistencies in ϕ c between CMC and Gala greatly affect Ẽ (Figure 2), and thereby the escape timescale.For example, a Keplerian fit ϕ fit c parameterized by the GC mass m(t) is too steep in the core relative to the true ϕ CMC c .This causes Gala to underestimate Ẽ(t rmv ), enough to prevent all but the most energetic PEs from crossing beyond r t (t) within a Hubble time.Meanwhile, a Plummer (1911) ϕ fit c set by m(t) and the GC half-mass radius causes the opposite issue; this choice is too shallow in the core and escape occurs far too rapidly, on the crossing timescale.Finally, though our simulations sample initial positions/velocities from a King (1966) profile, which also fits well ϕ c after core collapse, we find it is not an ideal match prior to collapse (most ages; explained shortly).In all cases, these discrepancies are far more visually apparent in the logarithmic enclosed mass profile as it magnifies changes in ∇ϕ c at small r.
We find a three-component Plummer potential-Equation ( 13)-provides a generally excellent fit to ϕ CMC c (and the GC's enclosed mass profile) across all t and r.We determined this via trial and error upon observing that once BHs segregate to the core, but before their eventual loss precipitates core collapse, the enclosed mass profile in CMC deviates strongly from a typical Plummer or King profile.Specifically, the mass density profile becomes bimodal due to the BHs forming a denser sub-cluster in the deep core, naturally leading to a two-component profile (one for BHs, one for stars).That a three-component fit also performs slightly better in Figure 2 is due more to the additional degrees of freedom.We use SciPy's curve-fitting functionality to fit such a potential to each simulation snapshot, linearly interpolating the fitted parameters both spatially and in time to yield ϕ fit c (r,t) for the Gala integration: The six fitted parameters m i (t) and b i (t) are the characteristic masses and Plummer scale lengths, respectively, for each piece of ϕ fit c , interpolated to time t.This definition guarantees the GC's enclosed mass tends to m(t) = 3 i=1 m i (t) as r → ∞.While a vast improvement over the above alternatives, the interpolated three-component Plummer fit still allows small inconsistencies between ϕ CMC c and ϕ fit c , enough to affect the Ẽ distribution in Figure 2.So we add a final step to correct for this.The equivalent of Equations ( 6)-( 8) after the fit/interpolation are: and To correct the inconsistency in Ẽ from fitting/interpolating ϕ c , we then slightly adjust v so that Ẽfit = ẼCMC rmv ≡ ẼCMC (t rmv ).The corrected speed that achieves this is The resulting Ẽ is Ẽcorr In the few cases (≲1 in 10 4 PEs) this results in imaginary v corr , we instead use v corr = v rmv .The (10th, 50th, 90th) percentiles of v corr /v rmv are ≈(0.99,1.00, 1.01).While small, this correction does impact Ẽ in Figure 2; in the top panel, the uncorrected Ẽfit (dashed light blue) deviates significantly from ẼCMC (solid gray).The correction works as intended, bringing Ẽcorr (dashed dark blue) in line with ẼCMC .For reference, we also show Ẽfit when instead setting ϕ fit c to be Keplerian, Plummer, or two-component Plummer (dotted red, yellow, and cyan lines, respectively).
Technically the above correction only accounts for inconsistency in Ẽ between CMC and Gala directly attributable to inaccuracy in fitting/interpolating ϕ CMC c .There are two other minor discrepancies.First, the raw criterion in CMC is only a spherical approximation, whereas the true Ẽ-defined in  6) and ( 9).The α criterion results in significantly higher Ẽ.
Gala as in Equation (3)-depends on the PE's full positional coordinates.So projecting the PEs into full six-dimensional phase space for initialization in Gala changes Ẽ from what CMC measured.The inconsistency has a small impact, causing ∼5% of CMC's 'PEs' to actually have Ẽ < 0 when measured under the true raw criterion in Gala (solid black).These mostly have low |z(t rmv )| and high |x(t rmv )|-see, e.g., W23's Equation (A4).To avoid bias, we exclude them from our analysis of the escape timescale in Section 5. Secondly, r CMC t and E CMC crit are both first-order approximations in µ (albeit excellent ones) to the actual location and ϕ eff of L1, so both criteria in CMC-Equations (6) and ( 9)-are only approximate whereas Equation (3) in Gala is exact to within machine precision.Other uncertainties in GC physics dwarf this tiny inconsistency, so, like all cluster-modeling codes we are aware of, we do not account for it in CMC.

RESULTS
Figure 3 shows the evolution of our archetypal NCC'd and CC'd GCs (solid and dashed curves, respectively) under both the raw (black) and α (blue) escape criteria.For reference, we also include the nearly-identical models from W23 (red), which used the α criterion.The top panels show the retention fractions of the initial number of particles N(t)/N i (left) and cluster mass M(t)/M i (right).The CC'd GCs (r v = 0.5 pc) evaporate faster than the NCC'd GCs (r v = 2 pc) due to their higher density and correspondingly shorter dynamical and relaxation timescales.As intended, the α criterion lowers the evaporation rate by raising the escape threshold-to Ẽ ≳ 0.1 (Figure 2).Yet the impact on the evaporation rate is smallless than the typical stochastic variation between separate statistical realizations of CMC models.This is unsurprising given that Giersz et al. (2008) saw relatively modest changes for N i = 10 4 , 80 times smaller than our N i .Since α scales inversely with N i in Equation ( 9) the two criteria are much closer in our case.Together, Figures 2 and 3 suggest that changing from the raw to α criterion primarily affects Ẽ, and thereby the escape timescale and escaper velocities, rather than the evaporation rate.So tuning the escape criterion to most accurately capture the latter (due to its greater relevance to internal GC evolution) may not be ideal for Galactic archeology applications, where the former are also important.
The lower left panel of Figure 3 shows a rolling average of the theoretical core radius r c , expressed as a ratio to the half-mass radius r h .The steep drop in r c /r h between 8 and 13 Gyr for the CC'd GCs demonstrate how their cores indeed collapse within a Hubble time, accompanying the transition from a centrally flat to a centrally steep surface brightness (an observationally CC'd state) upon loss of most BHs (e.g., Kremer et al. 2020Kremer et al. , 2021;;Rui et al. 2021b).The lower right panel shows the evolution of r 99 /r t , the normalized radius enclosing 99% of the GC's mass.Each GC significantly underfills its tidal boundary at birth but tends toward a tidally-filling state where r 99 /r t ≲ 0.626-the minimum clustercentric distance to the tidal boundary for a logarithmic Galactic potential (in the z-direction; Claydon et al. 2017).First, two-body relaxation (yellow) dominates overall, producing escapers with Ẽ ≪ 1.About half of escapers from relaxation originate within the typical core radius at removal r c (t rmv ), indicated by the vertical line and shaded interval in each scatter plot.This reflects that even bodies with Ẽ just below 0 in the GC halo typically first cross to Ẽ > 0 only after first plunging back through the core, where the higher density greatly enhances relaxation's efficiency (e.g., Spitzer & Shapiro 1972).Strong fewbody encounters dominate escape at high Ẽ and from the deep core (though relaxation still dominates in the core overall).Strong encounters are especially prolific in the CC'd GC (dashed curves), which also features several times more escapers from strong binary-single (light blue) and binary-binary (dark blue) interactions, as well as two-body relaxation.These reflect the increased density and correspondingly faster dynamics.Meanwhile, the faster loss of BHs in the CC'd clusters quenches three-body binary formation (from three singles; red) due its steep mass dependence (see W23).This mechanism dominates high-Ẽ escape prior to observable core collapse, corresponding to the present in most MWGCs (Trager et al. 1995).

Escaper Energy Distribution
The expression of energy as a fractional difference from E crit is an important qualitative difference from the similar figure in W23, since it emphasizes behavior at Ẽ ≪ 1.This reveals that two-body relaxation in CMC applies stronger kicks at higher density (smaller r).This arises because the average squared velocity kick applied to each body per (spatially uniform) timestep in CMC's relaxation algorithm is proportional to the local density-see Equation ( 9) of Rodriguez et al. (2022).This discretization introduces some uncertainty to Ẽ from relaxation since it is truly a continuous diffusive process.Ẽ from three-body binary forma-tion also decreases with increasing r because CMC limits the semi-major axis of the newly formed binary to a minimum a min ≡ Gm 1 m 2 /(50⟨m⟩σ 2 ), where m 1 and m 2 are the binary's component masses, ⟨m⟩ is the local average mass, and σ is the local velocity dispersion (see Section 2.3.1 of Rodriguez et al. 2022).So the maximum potential energy released in binding the binary is 25⟨m⟩σ 2 , and both quantities scale inversely with r.
Figure 5 complements Figure 4 by showing the distribution of Ẽ versus removal time t rmv .The first ejections occur primarily from strong binary-mediated scattering (blues) in both the NCC'd and CC'd GCs, but especially in the latter due to its higher density.Bursts of BH and neutron star ejections (black and teal, respectively) follow from 10-100 Myr, mostly due to supernovae.The smaller secondary burst of BH ejections arises from the BHs ejected by the supernova of a neutron star companion.A similar secondary burst of neutron star ejections at ≈60 Myr is a result of electron-capture supernovae, which occur a bit later and at lower ejection speeds than the dominant core-collapse supernovae.By t ∼ 100 Myr, escape occurs through a mix of twobody relaxation, three-body binary formation, and strong binary encounters, with many bursts of escapers at common times (vertical streaks) visible due to gravothermal oscilla- tions, a mathematically chaotic phenomenon in which the core sharply contracts and re-expands at frequent irregular intervals throughout the GC's life (e.g., Heggie & Hut 2003).These density spikes especially promote three-body binary formation since the rate for this process scales with the density cubed.A burst of strong binary-mediated ejections from 0.2 ≲ t/Gyr ≲ 1 occurs in the NCC'd GC due its unusually long and deep early core contraction (Figure 3) reversed by BH binary burning.The loss of almost all BHs in the CC'd GC curtails three-body binary formation around t ∼ 10 Gyr and induces core collapse, promoting strong binary-mediated ejections instead.
The typical Ẽ decreases over time-most notably for threebody binary formation since the average mass of bodies in the core drops as the GC ejects its BHs-but peaks again after core collapse in the CC'd GC due to the increased core density.Note in the right corner plot (identical between Figures 4-5) that Ẽ is cumulatively higher across all times and escape mechanisms in the NCC'd GC (solid gray) but the CC'd GC (dashed gray) has a comparable typical Ẽ when restricted to ages near a Hubble time.This is attributable to its high post-collapse core density and the correspondingly stronger relaxation kicks and late burst of strong encounters.

Escape Timescale
We now examine the distribution of escape times, ∆t esc ≡ t esc − t rmv , between removal from CMC (becoming a PE) and first passage beyond r t , at which point we may say the body has 'escaped.'While such 'escapers' can and often do circulate back within the GC's tidal boundary at least once before the GC's eventual dissolution, we examine the ramifications of these return trajectories later, focusing here on the timescale to cross beyond r t for the first time (e.g., FH00; Figure 6 shows survival functions for ∆t esc (left panels) and its normalized form from FH00, t = ∆t esc ω Ẽ2 (right panels).In the top two rows, we separately evolve each PE's trajectory for 14 Gyr past the time t rmv that it first satisfies the raw (top) or α (2nd from top) escape criteria, assuming constant ϕ c = ϕ c (t rmv ).This guarantees conservation of Ẽ while making escape take longer (by neglecting GC mass loss) and allows direct comparison to the results of FH00, who also assumed a constant ϕ c .In the bottom two rows (again, one for each criterion), we instead evolve each trajectory with t rmv < 13 Gyr for only 1 Gyr each in the true evolving CMC ϕ c (t).In this case, we do not know ϕ c (t) beyond the simulation end time (14 Gyr), so we impose the 1 Gyr cutoff to avoid biasing the ∆t esc distribution to smaller values (since many PEs removed at large t rmv will not have had time to escape).Again, solid curves denote the NCC'd GC and dashed the CC'd GC.The thick black curves represent the total survival across all Ẽ while the thinner rainbow curves each correspond to the partial contribution from each of several thin bands in Ẽ spaced 1/4 dex apart from log 10 ( Ẽ) ∈ [−2.5, 0.5] (dark red to violet, respectively; see caption for more detail).Curves are truncated because of the many escapers that do not cross beyond r CMC t within the integration time.Numerous interesting results are apparent in Figure 6, so we start by comparing the overall ∆t esc distributions (black curves; left panels).As pointed out by FH00, the dominance of two-body relaxation (low Ẽ) causes a large fraction (≳ 40%) of PEs under the raw criterion to escape on the Hubble timescale t H , at least when integrated in constant ϕ c (top row).On its own, this would appear problematic for In the top two rows, we integrate each PE's trajectory for a full 14 Gyr beyond the time trmv that it first satisfies the raw (top) or α (2nd from top) escape criteria, assuming constant ϕc = ϕc(trmv).In the bottom two rows (again, one for each criterion), we instead integrate all PEs with trmv < 13 Gyr for only 1 Gyr beyond trmv in CMC's truly time-dependent ϕc(t)-see explanation in text.Again, solid curves denote the NCC'd GC and dashed the CC'd GC.The thick black curve shows the survival function across all escapers while the remaining curves show the contributions from each of 13 thin bins in Ẽ.These are colored in rainbow order (right to left in the left panels) with bin centers Ẽcen uniformly spaced in log-scale 1/4 dex apart Ẽcen = 10 −2.5 (dark red) to Ẽcen = 10 0.5 (violet).Each bin's lower/upper bound is defined narrowly as ( Ẽlow, Ẽupp) = (0.9, 1.1) × Ẽcen.Truncation of curves correspond to exclusion of PEs that do not cross beyond r CMC t within the integration time.Finally, we show for comparison the fitted t distributions (dotted black) from Figure 9/Table 1 of FH00 (see text).
GC modeling, which often neglects phenomena occurring on t H -e.g., evolution of the Galactic potential.But removal under the raw criterion neglects the impact of ongoing weak two-body encounters during ∆t esc .By raising the Ẽ necessary to become a PE in an attempt to roughly account for these interactions, the α criterion dramatically reduces the probability of ∆t esc ≳ t H (2nd row from top).Further accounting for the impact of GC mass loss by integrating PEs in CMC's true evolving ϕ c (t) also hastens escape (lower two rows) but is less influential than the escape criterion because the initial dissolution timescale for most MWGCs surviving today is ≳ t H .For all four combinations of escape criterion and ϕ c assumption, ∆t esc overall is significantly shorter for the NCC'd GCs (solid) than the CC'd GCs (dashed).This is unsurprising since Figure 2 showed the former have higher Ẽ overall (due to more strong ejections from BH-driven threebody binary formation; Figures 4-5 and W23).Each is projected into the xy-plane and integrated for time 1.1∆tesc, just beyond when they cross rt (blue circle; note the projection only makes it appear some of these do not cross rt ).Each panel shows 20 such escapers, one highlighted in red for clarity, belonging to a different ∆tesc interval (see labels), selected to be immediately before/between the first four plateaus visible in the ∆tesc distribution at this Ẽ.Trajectories immediately before the N th plateau typically loop back N − 1 times before escape.We show in blue the L1 (left) and L2 (right) Euler-Lagrange points and in gray the two forbidden realms (of the 80 across all shown trajectories) enclosing the GC the least (dark gray) and most (light gray).This spread is due primarily to the finite Ẽ bin width when selecting the displayed escapers and secondarily to the variation in µ(trmv) across these, since they are removed from CMC at different ages.
A visually striking feature of Figure 6 is the appearance of successive plateaus in the ∆t esc distribution, especially at intermediate energies-the yellow, gray, and teal curves with respective Ẽ ≈ (0.1, 0.18, 0.32).These features are not evident in the results of FH00, Ernst et al. (2008), or Tanikawa & Fukushige (2010), perhaps in part due to their focus on low Ẽ, where ∆t esc is so long that these plateaus are minute.Yet both FH00 and Tanikawa & Fukushige (2010) show ∆t esc for similar bands in Ẽ as high as 0.24, so the phenomenon may be unapparent in their results simply because they cut out low ∆t esc , where the plateaus are most noticeable.Regardless, the plateaus are a true physical feature arising from the geometry of the zero-velocity surface.For 0 < Ẽ < 0.1 (dark red to orange curves), the necks are so narrow that PEs typically loop back many times through the GC before finally 'finding' and escaping through a neck.For Ẽ ≳ 0.5 (light blue to violet), the necks are so large that most PEs cross directly beyond r t .But for 0.1 ≲ Ẽ ≲ 0.5, especially near the center of that interval, escape often requires a small number of additional crossings before the body finally passes through either neck.We show in Figure 7 that bodies escaping immediately before the N th plateau in the ∆t esc distribution typically correspond to N − 1 such additional crossings.
That the plateaus become less obvious as N and ∆t esc increase relates to chaotic scattering theory, in which the basins (regions) of phase space leading to escape through different necks in an underlying conservative potential are separated from each other by a non-attracting fractal boundary known as a chaotic saddle (see, e.g., Ott & Tél 1993;Ott 2002;Tél & Gruiz 2006;Ernst et al. 2008;Seoane & Sanjuán 2013, and references within).Such fractals also appear in the phase space basins corresponding to different bins in the ∆t esc distribution (e.g., de Assis & Terra 2014;Zotos 2015aZotos ,b, 2016;;Zotos & Jung 2017).In general, though there are regions of phase space where trajectories are regular and ∆t esc is ef-fectively infinite (especially retrograde orbits in the CR3BP, e.g., FH00), N and ∆t esc increase where the phase space winds into finer (more chaotic) regions of the fractal boundary.Infinitely many locally smooth patches of each basin (for escape through either neck) exist in the saddle, each with a different characteristic ∆t esc that, when exceeded, induces a sharp drop off in the overall ∆t esc distribution as PEs from the patch exit the GC in a burst, leaving behind a plateau.The shrinking size and winding of these patches in the chaotic saddle as N and ∆t esc increase diminishes the bursts and blurs them together until, for ∆t esc → ∞, the surviving fraction of PEs decays roughly exponentially with ∆t esc .For an excellent discussion with graphics in the context of the Hénon-Heiles potential see Aguirre et al. (2001), especially their Section IV.A/Figure 9 relating to the plateaus.
Figure 6 also shows that for either constant or evolving ϕ c , changing the escape criterion negligibly alters the ∆t esc distributions for almost all specific Ẽ (excluding Ẽ ≈ 0.1 in yellow).This is reasonable; while exact phase space details affect ∆t esc , Ẽ controls the width of the zero-velocity surface's openings.We see more significant differences when comparing ∆t esc at specific Ẽ between the CC'd and NCC'd GCs.∆t esc is higher in the former at all Ẽ shown (especially Ẽ ≈ 0.1 in yellow) except Ẽ ≲ 0.02 (red curves), where the opposite is true.These discrepancies at identical Ẽ must relate to differences in the initial phase space distributions of PEs. Figure 4 showed Ẽ ≈ 0.1 corresponds mostly to relaxation deep in the core (r ≲ r c /10) while Ẽ ≲ 0.02 corresponds mostly to relaxation in the halo (r ≳ r c ).So relative to the NCC'd GCs, relaxation in the CC'd GC yields slower escape from the deep core and faster escape from the halo.
The faster escape in the halo of the CC'd GC may result from radial velocity anisotropy (bias to |v r | > v t ), which develops in the halos of GCs born centrally dense, especially near core collapse, but not in GCs born more diffuse (e.g., Giersz & Heggie 1997;Takahashi et al. 1997;Takahashi & Lee 2000;Baumgardt & Makino 2003;Tiongco et al. 2016;Zocchi et al. 2016;Claydon et al. 2017).PEs on such elongated orbits have long been known to escape more easily for 0 < Ẽ ≪ 1 since they probe the full tidal boundary, precessing to eventually find either neck.The stability of nearcircular retrograde orbits within or near the tidal boundary also promotes preferential escape of radial orbits.Meanwhile, the slower escape from the deep core of the CC'd GC may result from the faster loss of BHs, an important driving source of strong kicks, which increase the effective orbital eccentricity of the kicked body.So the early loss of BHs in the CC'd GC may dampen radial velocity anisotropy in the deep core, slowing escape, but firmer explanation of such subtleties will require further study.
Figure 6's right column shows the ∆t esc distribution in normalized form, t ≡ ∆t esc ω Ẽ2 .This is motivated by the approximate scaling relation ∆t esc ∝ ω −1 Ẽ−2 derived by FH00 for the most common case Ẽ ≪ 1 (see also Tanikawa & Fukushige 2010;Renaud et al. 2011).If this scaling is correct, t distributions at different Ẽ in this limit should be nearly identical.They should also roughly match the fitted t distributions from Figure 9/Table 1 of FH00 (dotted black) since Ẽ is independent of N-their Equation ( 9)-and weakly dependent on the central concentration or assumed ϕ g (Tanikawa & Fukushige 2010).But cumulatively, there are significant differences between FH00's setup and our own; the former assumed a Keplerian (instead of logarithmic) ϕ g and sampled PEs from a static King (1966) ϕ c while we sample PEs from (and for the lower panels, even evolve trajectories within) CMC's time-dependent ϕ c .This ϕ c features large density fluctuations (gravothermal oscillations), velocity anisotropy, and Ẽ drawing from numerous escape mechanisms.
Given the above differences, the relatively close match in the t distribution for Ẽ < 0.1 under constant ϕ c in Figure 6 is encouraging.For these curves (top two rows, red through orange), our t distributions converge and do not show a monotonic trend in the t space (e.g., gradual drift to lower t as Ẽ decreases).This supports the Ẽ−2 scaling of ∆t esc under constant ϕ c in FH00.Yet this is contrary to how t increases monotonically with Ẽ at Ẽ ≳ 1 (blue and violet).The FH00 scaling relation, based on the phase space flow rate near L1/L2, is not valid at such high Ẽ since the corrresponding flow is no longer restricted to the vicinity of L1/L2.This also largely explains the discrepancy between our overall t distributions in black and those of FH00.The tail at high t in the former (literally a sum over the continuous procession of the nearly vertical blue curves) arises from our inclusion of many escapers from strong encounters.The fits to FH00 are worse in the upper panels because many PEs that have not yet escaped are excluded in our curves.If our simulations were run longer, new 'escapers' beyond r t would be added at the top left in each panel, pushing the tail at high t down until it appears more like the thick black curves in the right panel second from the top.

An Empirical Escape Timescale for an Evolving Cluster Potential
When ϕ c (t) is allowed to evolve during the trajectory integration (bottom two rows of Figure 6), the t distributions for different Ẽ ≪ 1 no longer overlap, indicating ∆t esc no longer scales as Ẽ−2 in this more complex case.This is hard to see in Figure 6 because the truncation at 1 Gyr cuts out nearly the entire t distribution for each of the red and orange curves, but it is at least clear that for such low energies, t monotonically decreases with Ẽ.This occurs because the initial dissolution timescale of old MWGCs is ∼t H , so has negligible impact for moderate to high Ẽ, which have ∆t esc ≪ t H .But as Ẽ decreases and ∆t esc lengthens, GC mass loss grows more relevant, limiting ∆t esc .The pattern is much more evident in Figure 8, a duplicate of the lower two rows of Figure 6 that only shows escapers from t rmv < 4 Gyr, allowing us to truncate without bias the ∆t esc distribution at 10 Gyr. Figure 8 shows how GC mass loss causes the ∆t esc distribution at low Ẽ to converge; for Ẽ low enough that escape proceeds on t H , GC evaporation itself limits ∆t esc .
Please note however that the curves in Figure 8 are not directly comparable to Figure 6.Except for ages after core collapse, lowering the maximum removal time from 13 to 4 Gyr increases the typical Ẽ, greatly reducing overall ∆t esc in black.Due mostly to stellar winds in young massive stars and supernovae, GC mass loss is faster at early times (Figure 3), so more rapidly increases Ẽ(t) of PEs as they find their way out of the GC.This truncates the ∆t esc curves more sharply in Figure 8 than in Figure 6 for any specific, but sufficiently low, Ẽ.The distributions at higher Ẽ change less between these Figures since the correspondingly faster escape reduces the impact of GC mass loss on ∆t esc .
Figure 8 shows a very nearly uniform ≈0.47 dex gap between each of our red curves at different Ẽ ≲ 10 −1.5 ≈ 0.03 and a ≈0.40 dex gap between each of our blue/violet curves at Ẽ ≳ 10 −0.25 ≈ 0.56.Since the curves are 1/4 dex apart in Ẽ, these trends suggest that for a realistically evolving ϕ c , t ∝ Ẽ1.9 for 0 < Ẽ ≲ 0.03 and t ∝ Ẽ1.6 for Ẽ ≳ 0.56.So We provide no expression for 0.03 ≲ Ẽ ≲ 0.56 because it is apparent in Figures 6 and 8 that there is no clean power-law scaling due to the elevated importance of the specific initial phase space coordinates at these energies and resulting appearance of more complex features like the plateaus in both the ∆t esc and t distributions.
From the above expressions, we can redefine the normalized escape timescale for the case of an evolving ϕ c as tevol ≡ ∆t esc ω Ẽ0.1 if 0 < Ẽ ≲ 0.03 We verify these empirical relations in Figures 9 and 10, which show that tevol in Equation ( 19) leads to much better convergence in the normalized escape time distributions   18) for low Ẽ ≲ 0.03, the red curves.These converge much better to a common t distribution.at low and high Ẽ, respectively.We also note that we have tested power law exponents in the vicinity of the ones specified above and found them to indeed provide the narrowest convergence (especially the one for high Ẽ).

Formation of Tidal Tails and Stellar Streams
Having empirically re-examined the escape timescale, we turn our attention to the tidal tails and stellar streams from our simulated GCs.Future upgrades to CMC's escape physics will involve tuning via careful comparison to direct N-body models, so for now we simply demonstrate macroscopic features (including new returning tails) and show that escapers from CMC can already reasonably reproduce established tidal phenomena for circular GC orbits in a static, spherical MW potential, so long as PE trajectories are evolved collisionlessly.To give a sense of the timescale and avoid clutter, we cut off each trajectory 2 Gyr after it first crosses to r > r t .The PEs with low Ẽ (top left) escape through the necks in their forbidden realm (gray) near L1/L2, the points on the blue circle indicating r t (t).Increasing Ẽ expands the necks, allowing PEs to cross beyond r t with higher ẏ or ż (in/out of the page).The lower panels show how the forbidden realm's retreat from the xy-plane at Ẽ ≥ 1 enables more immediate return to the GC.At Ẽ ≈ 1 (lower left), in particular, the Coriolis effect causes ∼10% of escapers to temporarily return to r < r t -sometimes many  7, but the panels now distinguish subsets (of 50 escaper trajectories each) belonging to 8 bins in Ẽ ∈ [10 −2 , 10 2 ] from the archetypal NCC'd GC under the raw criterion and integrated to age trmv + ∆tesc + 2 Gyr.Except in the last panel (due to its high Ẽ), this time limit excludes the portions of the trajectories that return to the GC after entirely circumnavigating the Galaxy in the rotating frame, cleaning up the Figure considerably (see also Figure 12, identical to this one but in the inertial center-of-mass frame).Unlike in Figure 7, there is no additional filter for escapers belonging to specific windows in ∆tesc, so we make the bins ten times narrower-between 0.99 and 1.01 the Ẽ specified in each panel.This makes the difference between the least and most enclosing forbidden realm in each panel nearly imperceptible.times, akin to periodic extratidal orbits (e.g., Hénon 1969)before moving beyond several r t from the GC.
In the rotating center-of-mass frame of the GC and MW (Figure 12) the trajectories with Ẽ ≪ 1 complete a single epicycle every ∼10r t (as expected from, e.g., Küpper et al. 2008;Just et al. 2009)-longer for higher Ẽ.Many of the escapers with Ẽ∼10 2 (lower right panel) entirely circumnavigate the MW in this frame within 2 Gyr.Even the tra-jectories at Ẽ ≪ 1 do so in ≈8 Gyr, well within the dissolution timescale of most MWGCs.This, too, agrees well with the epicyclic approximation for low Ẽ, where escapers drift away from the GC along the tails at a speed v d ≈ 2ωr t for a logarithmic ϕ g -e.g., Equation ( 18) of Küpper et al. (2010).In our case, ωR gc ≈ 220 pc Myr −1 and R gc /r t ≈ 80, so v d ≈ 5.5 pc Myr −1 .Since the distance traveled per full or-bit about the MW is 2πR gc ≈ 5 × 10 4 pc, the drift period in the tails (timescale to return to the GC) is T d ≈ 9 Gyr.
The potential impact of return trajectories is apparent in Figure 13, containing several views of the PEs from the CC'd GC under the α criterion at age 12 Gyr (a snapshot from a full 14 Gyr-long movie linked in the caption).The lower right panel shows the projected positions in the orbital (XY ) plane of the inertial center-of-mass frame, and the other panels the three orthographic projections along each axis of the rotating clustercentric frame (see caption).Since weak twobody relaxation dominates escape at t = 12 Gyr, the streams closely follow the GC's circular orbit.This contrasts with ages ≲3 Gyr, when a clumpier, more energetic Ẽ distribution (Figure 5) leads to more irregular, branching streams (see the full movie).But perhaps the most notable feature-novel in the context of star cluster literature-is the appearance (lower right) of robust 'returning tidal tails', which form an X-like structure with the usual outgoing tails.This structure only appears at ages ≳T d , so coloring the PEs by t rmv emphasizes the large age difference between the outgoing (mostly red) and returning (yellow) tails.Crucially, however, the robustness of the latter is largely due to the assumptions of a circular GC orbit in a spherical, unevolving ϕ g .In reality, perturbations from MW substructure and an asymmetric, evolving ϕ g likely disperse such tails to much lower density.We further discuss these caveats in Section 6.1.

Stream Morphology
To examine stream morphology in more detail, we show in Figures 14 and 15 the surface number density Σ of PEs near the tidal boundary and along the full stream, respectively, from the CC'd GC under the α criterion (see the online journal for the other simulations).Σ is time-averaged between ages t/Gyr ∈ [11, 13], achieved by stacking 2001 snapshots of the PEs in that interval, finely binning the PEs by position, and dividing each bin count by 2001 and the bin area in pc 2 .The new clustercentric coordinates (x ′ , y ′ ), still in units of R gc , are flattened to map the full circular orbit to the line segment y ′ ∈ [−π, π], where (x ′ , y ′ ) = 0 corresponds to the GC's center.Specifically, The full stream consists of a leading and trailing tail sandwiching a low-density channel of width ≈0.036R gc ≈ 3.4r t (our GC models have average r t ≈ 85pc in the chosen age range).This channel closely follows the zero-velocity surface for 0 < Ẽ ≪ 1, which has width 2 √ 3r t ≈ 3.46r t -see Equation ( 17) of Just et al. (2009).The width of each of the adjacent tails is ≈2.3r t , in good agreement with the direct Nbody models and epicyclic approximation presented by Just et al. (2009); their Equations ( 17) and ( 34) result in tail width ≈2.1r t .Epicyclic overdensities are readily apparent, spaced ≈12.2r t apart.This is between the epicyclic approximations of 8.9r t (Küpper et al. 2010) and 15.4r t (Just et al. 2009)their respective Equations ( 20) and ( 22), given a logarithmic ϕ g has epicyclic frequency √ 2ω.Their disagreement arises from their differing assumptions on the starting point for escapers' epicyclic trajectories: x ′ = r t or x ′ = √ 3r t , respectively.Our results suggest an optimal approximation is intermediate to these extremes.Note the above distances in terms of r t are independent of GC mass M and R gc in the MW halo (Just et al. 2009), where ϕ g remains logarithmic.
The Figures' exceptional resolution also reveals a subtle feature common to all four simulations: double-ridged density peaks in the tails near each epicyclic overdensity.These arise from a phase difference in the trajectories of low-Ẽ escapers from opposing sides of each tail's neck (i.e., on either side of L1/L2); the corresponding overdensities streaking out from the tidal boundary are visible in the upper panel of Figure 14 (see also Appendix A).This phenomenon disperses each epicyclic overdensity compared to the typical 'particle spray' approximation (e.g., Küpper et al. 2008Küpper et al. , 2010Küpper et al. , 2012;;Just et al. 2009) that escapers all exit the GC at exactly L1/L2 (y ′ = 0).
While Figures 14 and 15 display variations between simulations (see the online Figure Set), these are relatively minor.Most notably, the α criterion reduces the PE number density at clustercentric distances 0.6 ≲ r/r t ≤ 1.This range, spanning the minimum and maximum r to the tidal boundary, is where the GC density deviates most significantly from spherical.The disagreement, resulting from faster escape under the α criterion, means comparison to matching direct Nbody models could determine which energy criterion best reproduces GC properties (e.g., density and velocities) within 0.6 ≲ r/r t ≤ 1.But given the asymmetry here, it is likely that no spherically symmetric criterion (even one based on energy and angular momentum, e.g., Spurzem et al. 2005) will allow our collisionless PE approximation to reproduce all GC properties of interest here on its own.So development of more nuanced escape physics in CMC, perhaps using basis expansion to allow the Monte Carlo method to handle mildly asymmetric ϕ c (e.g., Vasiliev 2015) may be worthwhile for this zone (see also Section 6.4).
Happily, changing between the α and raw criteria has lesser impact on features beyond the tidal boundary.For example, the former slightly widens the tail at the tidal boundary due to the higher typical Ẽ and correspondingly larger necks about L1/L2.Further from the GC, this widening is most perceptible again at the first epicyclic overdensities.But this latter difference is minute, especially given our extreme time-averaged sample size is well beyond that achievable by observations, even stacking tails from many MWGCs.So, unlike for the asymmetric but still collisional region within the tidal boundary, upgrades to escape in CMC should minimally impact the morphology of simulated stellar streams, in which the collisionless approximation is quite accurate.

Stream Density and Internal Velocities
While two-dimensional projections of stream properties are qualitatively revealing, it is easier to study variations between simulations and the impact of the returning tails via quantities averaged along the stream axis.As bulk properties, the stream density and internal velocity profiles are an ideal starting point.To this end, Figure 16 shows the surface number density Σ of PEs (when viewed from the MW's center; top panel), their mean speed ⟨v − v c ⟩ with respect to the GC's circular speed (middle), and the dispersion in this speed σ The epicyclic overdensities are again apparent and Σ is nearly symmetric between the leading and trailing tails.The former is only ≈5% denser since the neck at L1 opens at slightly lower Ẽ than the neck at L2.The density in both tails gradually decreases as they extend further from the GC before reaching a minimum and increasing again about a quarter orbit before returning to the GC.In the earliest portion of each tail's outgoing half-y ′ > 0 (leading) and y ′ < 0 (trailing)-Σ is about twice as high from the CC'd GC than the NCC'd GC, consistent with their evaporation rates in Figure 3.This factor shrinks closer to 3/2 just before the tails return to the GC, consistent with steady leakage of PEs from the stream over time.Finally, since the circular speed ωR at any Galactocentric distance R is constant in a logarithmic ϕ g , then ω is slightly higher in the leading tail.This explains why the crossing point between the tail densities apparent at the far left of the top panel occurs at y ′ /π just slightly > − 1.
Under the epicyclic approximation, the mean speed in the stream relative to the GC's circular speed is ⟨v − v c ⟩ ≈ [(2ω/κ) 2 − 2]ωr t , where κ is the epicyclic frequency-see Equation ( 21) of Küpper et al. (2010).For logarithmic ϕ g , κ = √ 2ω, so ⟨v − v c ⟩ ≈ 0. The central panel of Figure 16 reproduces this expectation along the entire stream, and locally to within ∼1 km s −1 .⟨v − v c ⟩ locally peaks/troughs between the epicyclic overdensities in the leading/trailing tails since most (low-Ẽ) PEs here have a velocity near directly parallel/anti-parallel relative to the GC's velocity at the same y ′ .At each epicyclic overdensity, however, most PEs are briefly moving back toward the GC-or at least counter to the velocity of their epicyclic trajectory's guiding center, which is offset from the GC's orbit.This instead causes ⟨v − v c ⟩ to trough/peak, respectively.These features remain, albeit with lower magnitude, in the combined stream profile (black) since each tail's returning half-y ′ < 0 (leading) and y ′ > 0 (trailing)-is less dense than the outgoing half.
The amplitude in the ⟨v − v c ⟩ oscillations is slightly greater from the CC'd GC than the NCC'd GC.This is counter-intuitive since (to recap Section 5.1) the NCC'd GC has cumulatively higher ejection Ẽ due to longer retention of BHs, which promote strong ejections from three-body binary formation (W23).Though we also noted the CC'd GC has comparable or greater typical Ẽ when limited to t ≳ 10 Gyr, due to the higher post-collapse density's promotion of strong encounters and relaxation, this turns out not to explain the CC'd GC's higher ⟨v − v c ⟩ amplitude in Figure 16.(As we show shortly in Figure 17, restricting the stream to only PEs ejected after 10.5 Gyr instead causes the CC'd GC to have minutely lower ⟨v − v c ⟩ amplitude than the NCC'd GC).Rather, the CC'd GC's faster mass loss-twice that of the NCC'd GGis likely responsible, since the mass loss decreases |ϕ c |, causing the PEs to lose less speed while traveling 'up' the GC's potential well to escape.
While the stream's ⟨v − v c ⟩ profile is very similar between the models, there is a larger difference in the profile of its velocity dispersion σ (lower panel).Since σ should correlate with the local tail width, it is unsurprising that σ peaks between epicyclic overdensities and troughs within them (based on the overdensity locations in Figure 15).Further along the tail, σ both increases and oscillates less since small deviations in the typical epicycle period for different Ẽ compound to randomize the velocity vector far from the GC (apparent in how the epicyclic overdensities distort and overlap far from the GC in Figure 15).On average σ is ≈20% higher from the NCC'd GCs than from the CC'd GCs.This is a direct result of including the returning tails, which allows bodies ejected  16 (including being a stack of 2001 snapshots spanning ages 11-13 Gyr) except for the following changes.To simulate the impact of stream disruption, we now exclude the PEs in each snapshot for which >500 Myr have elapsed since their 'escape' (first crossing beyond r > rt ).We also only show the sum of the leading/trailing tails (i.e., the black curves from Figure 16) and do so for all four simulations.The solid/dashed curves still correspond to the NCC'd/CC'd GCs, and the blue/red curves to the raw/α escape criteria.Finally, because of the limited tail length from the time cutoff, we use finer bins in y ′ -2001 across the entire GC orbit, corresponding to bin widths ≈25 pc.
many Gyr ago to contribute to Figures 14-16, despite only stacking snapshots of the stream in the age range 11-13 Gyr.The NCC'd GC's wider spread in Ẽ before ∼10 Gyr, due to more ejections from BH-driven three-body binary formation, thus inflates σ in the NCC'd GC's stream even at late times.
Due to the impact of the returning tidal tails and their likely disruption by MW substructure in a realistic ϕ g , we reproduce Figure 16 Figure 17's primary difference with respect to Figure 16 is that the amplitude in the ⟨v − v c ⟩ oscillations is ≈2 times higher while the dispersion σ is ≈3 times lower.Eliminating the returning tail increases the former because there is no longer a returning flow opposing the outgoing tail (given our chosen viewpoint from the MW's center).That opposing flow significantly increases ⟨v − v c ⟩-compare, for example, the blue/red curves to the combined black curve in Figure 16, or, for a helpful visual aid, see maps of the two-dimensional projected ⟨v − v c ⟩ in Appendix A. These maps show that neglecting returning tails reduces σ for the same reason, though the lower dispersion in Ẽ at late times (Figure 5) likely helps.
Unlike when including the returning tails in Figure 16, neglecting the returning tails in Figure 17 results in there being little difference in σ between the tails at t ∈ [11, 13] Gyr from the NCC'd GCs versus those from the CC'd GCs.This is because the NCC'd and CC'd GCs have more similar Ẽ distributions at t ∈ [11, 13] Gyr (unlike their Ẽ distributions cumulatively).In particular, looking back at Figure 5 for that age range, we see that though the NCC'd GC still has more high-Ẽ ejections from three-body binary formation (due to longer BH retention; W23), the CC'd GC now has significantly more high-Ẽ ejections from binary-single and binarybinary interactions (due to the higher density after core collapse).As a result, the total number of high-Ẽ escapers contributing to σ within this age window is roughly equivalent between the NCC'd and CC'd GCs, leading to a similar Ẽ distribution and σ within the tidal tails.
Together, Figures 16 and 17 show that the average speed of bodies in the tail ⟨v − v c ⟩ is not a reliable measure of a GC's state of core collapse and by extension its presently retained BH population.Yet ⟨v−v c ⟩ and its dispersion are good indicators of the presence of a returning tidal tail, potentially making these quantities a useful constraint on the MW potential and the number and properties of massive perturberse.g., giant molecular clouds and dark matter subhalos-in the MW halo, since these control stream disruption.
6. DISCUSSION 6.1.Implications and Caveats of Returning Tails A novel finding in this work is the appearance of returning tails in our simulations.When interpreting these features, it is crucial to recognize that our idealizations of a circular GC orbit in an unevolving, spherical MW potential are a bestcase scenario for the stream density.In reality, they should be much more diffuse as asymmetry, substructure, and timedependence-including perturbations from the MW disk and bar, molecular clouds, and dark matter subhalos-should all naturally disrupt streams, blending them into the stellar background.Even detection of the denser outgoing tails remains difficult; the number of GCs observed to have them remains low (≈15; Piatti & Carballo-Bello 2020).In none of these cases are returning tails apparent.
The MW is also not an especially accommodating place for returning tails as they need a drift period T d less than the GC's age to form.In the MW halo, where the circular speed v c is very flat (enclosed mass gc .So bulge substructure, perhaps aided by a rotating bar (e.g., Hattori et al. 2016;Pearson et al. 2017), would have ample time to disrupt the streams.Eccentricity or even precession of the GC orbit-induced by, e.g., a triaxial ϕ g (Capuzzo Dolcetta et al. 2005)-could also misalign the bulk of the returning orbits from the GC at most phases of its orbit, preventing consistent returning tails at the tidal boundary.Even so, the sheer number of escapers (often >10 5 ) and their wide dispersal on the MW's crossing timescale (Figure 13), means that even major disruption or misalignment of the stream is unlikely to prevent some escapers from returning to their original host GC after circumnavigating the MW.Tidal capture also requires low Ẽ (Koon et al. 2011), so re-capture of past escapers from a dispersed stream may still exceed capture from the MW's stellar background.It is merely unclear whether more diffuse returning tails would be observable.
Yet detection of returning tails is not hopeless.Clearly an approach reliant on stellar surface density alone would be hindered by the background, and also the narrowness of the gap between the outgoing and returning tails; from many viewing angles no gap is apparent (e.g., the lower panels in Figures 14-15).But kinematic detection is more promising; thanks to kinematic (and chemodynamic) measurements, many MW streams are already known to be extremely elongated (Mateu 2023).So while no streams todate have been traced over more than a full orbit of the MW, this may be achievable with future surveys.Notably, we found in Section 5.6 that the opposing flow of the outgoing and returning tails boosts the stream's velocity dispersion by ≈3 km s −1 -substantial since the stream's dispersion without accounting for the returning tails is only ≈1 km s −1 .This boost is similar in magnitude to that from perturbations by MW substructure, such as dark matter subhalos or the MW disk (e.g., Carlberg 2009;Carlberg & Agler 2023).It is also much greater than the likely inflation of σ by unresolved binaries in streams.Multi-epoch radial velocity measurements of stars in streams (Li et al. 2019, with the S5 Survey) and in the MW halo (Conroy et al. 2019, with the H3 Survey) indicate that systematic uncertainties from phenomena like binarity inflate observed σ by an amount similar to the known uncertainties from the surveys' pipelines: ≈0.5 km s −1 .So while the true impact of returning tails on σ in streams should be substantially weaker than the 3 km s −1 boost in our ideal case, returning tails may be relevant on a level comparable to observing biases such as unresolved binaries.

Black Hole Mergers from Return Trajectories
In Section 5.4, we noted many escapers at moderate Ẽ ∼ 1 promptly return to r < r t , even multiple times, before moving more than several r t from the GC.An intriguing ramification is that BHs 'ejected' from the GC may pass back through the core and undergo a strong encounter, potentially re-binding to the GC.In the case of BHs ejected by gravitational-wave merger kicks, this may enhance dynamical production of hierarchical BH mergers (e.g., Miller & Hamilton 2002;Rodriguez et al. 2019).Checking all four simulations, we find that within the 14 Gyr simulation runtime, BH escapers pass from r > r t to r ≤ r t an average of 200 times per GC-6 times when limited to BH merger remnants.These drop to 3.5 and 0.25 times, respectively, when considering only core passages from r > r c to r ≤ r c (the density-weighted core radius from Casertano & Hut 1985).Assuming a typical r c ∼ 1 pc containing stellar number density 10 3 pc −3 and a BH remnant with mass 30 M ⊙ and speed 50 km s −1 , then the rate of strong encounters between returning BH merger remnants and typical stars in the core (capable of scattering the BHs back to Ẽ < 0) is only ∼10 −8 per GC per Hubble time.This suggests that returning BH merger remnants negligibly enhance hierarchical BH merger rates in GCs.

How Cluster Density and Black Holes Affect Streams
While external tides dominate stellar stream morphology, factors internal to GCs are also relevant.In particular, the GC's initial number density n (set by r v and N in CMC) determines the timescales for both relaxation and stronger encounters (e.g., Spitzer 1987;Heggie & Hut 2003;Binney & Tremaine 2008).Along with R gc , n also determines how fully the GC's mass distribution fills its tidal boundary.These considerations have competing effects; higher n leads to faster relaxation, hastening evaporation, but also leads to a deeper central potential and so a larger gap between E crit and the typical Ẽ in the GC.Since most bodies achieve Ẽ > 0 in the GC's core this means higher n can also hinder escape, slowing evaporation.The competing effects can lead to seemingly divergent conclusions in the literature on n's impact on evaporation.The speed-up of relaxation dominates in this study, so our denser (CC'd) models evaporate faster, yielding denser streams.This holds for typical CMC models, which have initial n, N, and R gc leading to evolved GCs consistent with most MWGCs today (Kremer et al. 2020).But lower density can lead to faster evaporation in CMC for GCs that more fully fill their tidal boundary (low R gc and high r v or low N; see W23's Figure 1).Relevant to diffuse, low-N GCs in the MW halo, such as Palomar 5, this regime has been more strongly emphasized in recent direct N-body modeling (Gieles et al. 2021;Gieles & Gnedin 2023;Roberts et al. 2024).
BH dynamics is similarly relevant as central heating from BH binary burning supports the GC against collapse, lowering the core density.As our models are tidally under-filling, our denser CC'd GCs evaporate faster and form denser streams despite ejecting their BHs earlier.But in the tidallyfilling regime, BH heating increases the mass loss rate and stream density.Gieles et al. (2021) showed this may explain Palomar 5's unusually dense tidal tails, finding that low-n direct N-body models could reproduce the tails if BHs make up ≈20% of Palomar 5's current mass.More broadly, variation in BH retention may help explain why only some MWGCs have observed tails (Gieles & Gnedin 2023;Roberts et al. 2024).To holistically investigate this prospect and extrapolate to the full MWGC population, future studies of the impact of GC density and BH dynamics on tidal tails should consider both tidally-filling and underfilling GCs.
By extension, tidal tails and streams in the MW halo may be useful probes of not just BH populations in GCs but also more fundamental uncertainties regulating BH formation and retention, such as the stellar initial mass function (IMF) and supernova kick strengths.Notably, GCs born with topheavy IMFs can produce so many BHs-and so much BH burning-that they rapidly dissolve in ∼1 Gyr (e.g., Banerjee & Kroupa 2011;Whitehead et al. 2013;Contenta et al. 2015;Chatterjee et al. 2017;Giersz et al. 2019;Weatherford et al. 2021), potentially producing very dense tidal tails.That the 20% BH mass fraction in Palomar 5 suggested by Gieles et al. (2021) is already higher than achievable from a canonical stellar IMF supports that tails may be a useful constraint on the IMF in low-density GCs.The accelerating evaporation rate near the end of such GCs' lives may further help explain why the density of some tidal tails (e.g., Palomar 5's) rises steeply with decreasing distance from the GC.We plan to investigate such prospects further in future work.

Collisional Dynamics for Potential Escapers
We have compared the raw energy criterion to the Giersz et al. (2008) α criterion, designed to account for scattering of PEs back to Ẽ < 0 prior to escape.Though motivated by conservation of Ẽ in the CR3BP, which does not hold for collisional dynamics, non-circular GC orbits, or evolving ϕ c or ϕ g , such energy-based escape criteria remain useful when the energy change is slow compared to the crossing timescale and ϕ g is not highly substructured near the GC.The intent of our comparison is not to identify which criterion yields more accurate streams, a judgement requiring direct N-body codes, but rather to gauge how much collisional dynamics for PEs during their escape may affect ∆t esc and stream properties.
Summarizing our findings, the α criterion raises the Ẽ of PEs (Figure 2), reducing ∆t esc (Figure 6) and lowering by several times the surface density of PEs in the portion of the GC's halo that is significantly non-spherical (0.6 ≲ r/r t ≲ 1; Figure Set 14).So collisional dynamics for PEs may significantly alter the morphology and kinematics of the GC within the tidal boundary.Yet the two criteria yield very similar evaporation rates Ṅ (Figure 3) and extratidal properties, including the stream density and velocity dispersion profiles.This is expected for slow, steady evaporation (∆t esc ≪ N/ Ṅ), as in our models representative of typical MWGCs.The small extratidal impact implies collisional dynamics can reasonably be neglected for PEs when simulating tails and streams from MWGCs, though perhaps not for r < r t , or small N and fast Ṅ (e.g., near GC dissolution).
To the extent that better accounting for PE collisional dynamics may still have small effects, keep in mind that weak two-body scattering is a competition between cooling (dynamical friction) and heating (relaxation)-e.g., Section 7.8 of Binney & Tremaine (2008).By scattering some PEs back to Ẽ < 0 before they can escape, cooling should reduce the evaporation rate, which the α criterion achieves by raising the threshold for removal from the GC dynamics to Ẽ ≳ 0.1.But this does not account for how cooling would also lower the speed of the PEs that do escape; by raising the threshold Ẽ, the criterion does the opposite.This may unintentionally help account for heating of PEs during escape, but whether the α criterion does so accurately is unclear.The accuracy should at least vary with PE mass, since cooling dominates at high mass and heating at low mass (most stars).
A more physically-motivated alternative is delayed escape, in which PEs continue participating in collisional dynamics before removal from the GC simulation.This option has been implemented into MOCCA by Giersz et al. (2013).MOCCA identifies PEs every timestep ∆T via the raw criterion.Based on FH00's results, it estimates the probability the PE will escape during ∆T as P(∆T ) ≡ 1 − (1 + bω Ẽ2 ∆T ) −c , where b ≈ 3 and c ≈ 0.8 are tuned to best fit the mass loss rate from direct N-body codes.Random sampling with this probability determines which PEs to remove each ∆T .Yet the scaling ∆t esc ∝ Ẽ−2 from FH00 neglects GC dynamics and mass loss and only applies for Ẽ ≪ 1.So ∆t esc from MOCCA's algorithm does not necessarily scale properly with Ẽ.Our results accounting for GC evolution show that ∆t esc ∼ ( Ẽ−0.1 , Ẽ−0.4 ) for PEs with (low, high) Ẽ, with some ambiguity in between due to chaotic scattering.In principle, these new scalings can be swapped with the Ẽ−2 relation in MOCCA's algorithm, and the coefficients b and c re-tuned against direct N-body codes.
While MOCCA's delayed escape prescription could mildly improve accuracy for escaper kinematics, it has not been implemented into CMC since the impact on the GC mass loss rate is small.And like the simpler energy criteria with immediate escape, it still limits the Hénon method's capacity to simulate more complex aspects of stellar stream morphology arising from evolution or substructure in the tidal field.More general alternatives are therefore necessary.

Generalizing to Non-spherical Clusters and Evolving Tides
There are many sources of evolving tides that we do not account for in this work.The orbits of most MWGCs are eccentric and inclined relative to the MW disk (Baumgardt et al. 2019), traits that induce tidal shock heating from passage near the MW's center or through its disk, respectively.This hastens GC dissolution and causes the GC's bound mass to fluctuate (e.g., Gnedin & Ostriker 1997;Baumgardt & Makino 2003;Webb et al. 2013Webb et al. , 2014a,b;,b;Madrid et al. 2014).The eccentricity also causes the tidal tails to fan (e.g., Küpper et al. 2008Küpper et al. , 2010) ) and more realistic, non-spherical MW potentials add further complexity.For instance, a triaxial ϕ g (e.g., Capuzzo Dolcetta et al. 2005) produces further stream gaps/overdensities and causes the GC orbit to precess, while a rotating Galactic bar induces further fanning and asymmetry in the lengths of the leading and trailing tails (Hattori et al. 2016;Pearson et al. 2017).Modifications to gravity, an alternative to dark matter, can cause similar asymmetry (Thomas et al. 2018;Kroupa et al. 2022).Finally, as noted in Section 1, finer MW substructures such as dark matter subhalos, molecular clouds, other GCs, and infalling MW satellites can all heat and strip stars from streams, leaving behind significant gaps and kinks.These introduce fine timedependence to the tides experienced by a GC, but tidal evolution on longer timescales is also likely, especially for the many GCs in the MW halo suspected to have been accreted from disrupted MW satellites (e.g., Helmi 2020).
Our motivations for developing the Hénon method for stream modeling are its detail for internal stream populations and its speed-and thus its capacity, in principle, to generate large model grids of streams to explore the impacts of the many above competing factors.To realize this prospect, improvements to CMC's escape criteria and tidal physics must be sufficiently general while not dramatically slowing CMC.
One option is to allow non-spherical ϕ c and dynamics via basis expansion and position-dependent velocity diffusion coefficients (Vasiliev 2015).While this only enables mild (factor of 2) asymmetries, so cannot on its own resolve tidal tails, it accounts for ϕ c 's asymmetry and allows PEs to continue participating in dynamics until escape.Collisionless trajectory integration could then follow as in our approach.But critical drawbacks accompany these advantages: (1) The method requires following each body's orbit and applying diffusion on the dynamical timescale, making it far slower than the Hénon method.(2) The asymmetry makes assigning pairs of nearest neighbors highly non-trivial, hindering implementation of strong encounters.(3) Updating CMC to use the method or adding many years of CMC development (stellar evolution, strong encounters, escape) to Vasiliev's code would be time-consuming.The benefits for stream modeling are not likely enough to outweigh these drawbacks.As we argued above, explicit treatment of collisional dynamics for PEs should not substantially alter stream properties and ϕ c 's mild asymmetries in the GC halo are not as relevant to streams as external tides or dynamics in the core, where most diffusion to Ẽ > 0 occurs (Spitzer & Shapiro 1972;W23).
To generalize the Hénon method to evolving tides, Sollima & Mastrobuono Battisti (2014) developed an escape prescription for non-circular GC orbits in two distinct cases of a static, cylindrically symmetric MW potential: a pointmass ϕ g and ϕ g with bulge, disk, and halo components.At each timestep ∆T , for each body in the GC, their algorithm isotropically samples a random escape vector, along which it computes the distance (and ϕ eff at) the tidal boundary.It flags as PEs bodies with sufficient energy and angular momentum to cross the boundary during ∆T , based on each body's orbital period in the GC.(Note this is mildly improper as it assumes escape occurs within several crossing times).The orbits of the GC and all PEs are then integrated within the external ϕ g for a full orbital period of the GC in the MW.PEs are only removed from the simulation if they pass beyond the maximum r t of the GC during this integration, allowing PEs to continue participating in dynamics prior to removal.Drawbacks are that the PE trajectories do not account for ϕ c and that the above selection of PEs is based on laborious calculations (in their appendix) for their specific choices of ϕ g , an approach hard to generalize to more substructured ϕ g .
More recently, Rodriguez et al. (2023) developed new CMC escape prescriptions for non-circular GC orbits and applied them to an evolving, substructured ϕ g from a FIRE-2 MHD simulation of a MW-like galaxy.They first extracted time series data on the tidal tensor (spatial 2nd derivatives of the local ϕ g ) from the orbits of particles representing GCs in the galaxy simulation, then ran CMC computing r t every timestep based on that data (loaded as an input file).When loading tidal tensor data from Gala integrations of eccentric GC orbits in a static ϕ g (with bulge, disk, and halo components), CMC reproduced well the mass evolution of equivalent direct N-body models (Webb et al. 2014a)-at least when using r apo > r t rather than an energy criterion.But mass loss was slower/faster than in the direct N-body models for low/higheccentricity GC orbits.These respective deviations are at least partly due to Ẽ > 0 not universally implying r apo > r t , the maximum distance to the tidal boundary, and the algorithm's inability to re-capture some escaped bodies when the boundary re-expands after each perigalacticon.
In the near term, when applying CMC to model streams from GCs on eccentric or inclined orbits, or in evolving, substructured ϕ g , we will blend our collisionless integration of escape trajectories-similar to Sollima & Mastrobuono Battisti (2014)-with the tidal tensor technique of Rodriguez et al. (2023).Specifically, we intend to identify PEs like Rodriguez et al. (2023), immediately removing them from CMC (neglecting collisional dynamics for PEs during their escape due to its small impact on stream properties).We can then directly integrate the PE trajectories each timestep in post-processing as in this work.This is advantageous since it will not entail building into CMC a separate collisionless integrator for escaping/escaped bodies.Unlike in this work, the orbit integration would occur in an inertial Galactocentric frame since the pseudo-forces (centrifugal, Coriolis, Eulerian, etc.) in the non-inertial GC frame are not generally known for an evolving, substructured ϕ g .
In the longer term, we aspire to add a galactic dynamics integrator to CMC to allow collisional dynamics (and delayed escape) for PEs and re-capture of escaped bodies.Due to the lack of a general definition for Ẽ in an evolving, substructured ϕ g , and to avoid missing any bodies capable of escape, we would select PEs via a more liberal version of the apocenter criterion than used by Rodriguez et al. (2023)-e.g., r apo /r t > 2/3.This corresponds roughly to the minimum distance to the tidal boundary, and the radius beyond which the GC is significantly non-spherical.Rather than being removed immediately from CMC, PEs flagged in this manner would then be randomly projected into full six-dimensional phase space coordinates (as in this work) and their trajectories integrated in the full combined ϕ c + ϕ g with the integration time set to CMC's current timestep.At the end of the integration, CMC would remove the PE from its collisional dynamics if r/r t > 2/3, but if r/r t ≤ 2/3 the body would stay in the GC to take part in dynamics in the next timestep.(Repeating this procedure every timestep is appropriate since CMC operates on the relaxation timescale, over which bodies' orbits are randomized.)Collisionless integration of the escapers would continue from there, separately but still internal to CMC.At the end of every timestep, escapers that return within r/r t < 2/3 could then be transferred back to the main collisional dynamics routine.This would allow re-capture, unlike Sollima & Mastrobuono Battisti (2014) or Rodriguez et al. (2023), but implementation would not be trivial as it requires expanding CMC's parallelization scheme to escapers.

SUMMARY AND FUTURE WORK
We have applied for the first time orbit-averaged star cluster models to study the formation of tidal tails and stellar streams from GCs.Specifically, we use CMC-a stateof-the-art, publicly-available implementation of the Hénon (1971a,b) Monte Carlo method.Though this method assumes spatial spherical symmetry in its collisional dynamics, we show that treating energetically unbound bodies (potential escapers; PEs) as collisionless enables formation in postprocessing of asymmetric tidal phenomena such as tidal tails with exquisite detail.The benefits of the Hénon method are that it is much faster than direct summation N-body methods for GCs of typical N and density, making it far better suited to large-parameter-space modeling of tidal tails and stellar streams from GCs-an essential feature due to the vast array of considerations internal and external to the GC relevant to stream properties.While even faster but more approximate stream simulation is possible via common particlespray techniques, CMC can generate streams with many of the advantages of a full GC-modeling code.These include highly accurate detail on time of escape, escaper kinematics, and details on streams' internal populations (e.g., masses, star types, binary properties, etc.).CMC may be especially ideal for rapid but reliable exploration of tidal tails from large grids of non-standard GCs, such as those with atypical stellar IMFs, extreme density, unusual binary fractions, or intermediate-mass BHs (work in progress).Key limitations are that the Hénon method's assumptions of evolution on the relaxation timescale, virial equilibirum, and high N ≳ 10 5 (for realistically broad stellar IMFs), make it unsuited for modeling tidal shocking or the final stage of GC dissolution.Yet as discussed in Section 6, further refinements can allow the Hénon method to simulate streams from GCs on non-circular orbits, even in an evolving, highly substructured external potential.
The main findings of this work, the second in our series on escape from GCs, are as follows: 1. We demonstrate that the Hénon method can accurately reproduce known features of tidal tails and stellar streams from GCs on circular orbits within an unevolving, spherical Galactic potential.This is achieved by collisionlessly integrating the trajectories of potential escapers in the full tidal field.
2. We examine for the first time the in-cluster survival timescale (escape timescale ∆t esc ) of PEs in a realistically evolving GC potential, and also account roughly for the impact of ongoing collisional dynamics on PEs via the Giersz et al. (2008) α escape criterion.Escape in this case occurs on a timescale of ≈100 Myr instead of ≈10 Gyr in the static, collisionless case.
4. We highlight discreteness in the ∆t esc distribution arising from chaotic scattering.Different characteristic ∆t esc within distinct locally smooth regions of the phase space of PEs introduce successive plateaus in the ∆t esc distribution for 0.03 ≲ Ẽ ≲ 0.56.This hinders a clean power-law scaling of ∆t esc on Ẽ in this interval.
5. We analyze for the first time the impact of return trajectories circumnavigating the Galactic center, finding they produce robust returning tidal tails on a timescale of ∼10 Gyr for GCs at R gc = 8, kpc in our idealized circumstances.Though a realistically evolving Galaxy with significant substructure is likely to disperse such tails, they may eventually be observable in proper motion space and could excellently constrain the history and substructure of the Galaxy over longer timescales than typical streams.Returning tails may be more relevant for GCs in dwarf galaxies, where the timescale for low-Ẽ escapers to return to the GC is much shorter.
6.In our ideal case, the returning tails increase velocity dispersion in stellar streams by several km s −1 , an effect similar in magnitude to the perturbative influence on streams by giant molecular clouds, the Galactic disk, and dark matter subhalos.Though the velocity dispersion enhancement is likely much smaller for streams dispersed by a realistically evolving and substructured Galaxy, the boost to velocity dispersion may still be comparable to the ≲ 0.5 km s −1 boost from other observing biases such as unresolved binaries.
In future work, we will further refine escape prescriptions in CMC to model streams from GCs on non-circular orbits in an evolving, substuctured Galaxy.To do so, we will combine our collisionless integration of escape trajectories with the general computation of an approximate spherical tidal boundary relying on tidal tensors (Rodriguez et al. 2023).We will also study in more detail the internal stellar populations inhabiting the tidal tails and streams from our simulated GCs.This will include an analysis of observability by continuing stellar evolution for escapers and converting their luminosities to common observing magnitudes (e.g., with Gaia), following the technique of Rui et al. (2021a).Other focii will be binaries in our streams (relevant to velocity dispersion biases) and mass segregation both parallel and perpendicular to the stream axis, which respectively would bias the observed length and width of streams due to lower observing completeness at lower mass.The former case (Webb & Bovy 2022) naturally develops since GCs tend to eject low-mass stars faster, and therefore earlier on average.Mass segregation perpendicular to the stream has not been considered before, but may arise from velocity-dependent ejection speeds of stars-e.g., due to the strength of the velocity kick scaling inversely with mass, or due to the higher rate of strong encounters for high-mass stars segregated in the GC's core.Such analysis would improve our understanding of observing biases relevant to stream morphology and kinematics, useful to reducing uncertainties when applying stream observations to, e.g., place firm constraints on the nature of dark matter.In summary, CMC's speed and aforementioned capabilities/planned improvements should enable us to construct for the first time large grids of stream models using a ded-icated GC-modeling code, along with detailed comparisons to stream observations.

Figure 1 .
Figure 1.Slices in the xy-plane of the rotating, clustercentric frame, illustrating features of the CR3BP for the case of Keplerian ϕg and ϕc, µ = 1/101, and in units where G = ω = 1.The panels distinguish different choices of excess relative energy Ẽ, in terms of the values Ẽi at each Euler-Lagrange point(Li)  when v = 0.These points and the centers of the Galaxy (G) and cluster (C) are labeled in the first panel but appear in all.The forbidden realm (gray) is bounded by the zero-velocity surface and cannot be entered by bodies with the given Ẽ or lower.The blue curves in the last four panels exemplify trajectories at the given Ẽ initiated with (x, y, z) = (µ, 0, 0), velocity v0 ∥ ŷ, and integrated for 20 full orbits of the cluster in the Galaxy.Upper left: Ẽ < 0, so bodies cannot transit between the cluster interior, Galactic interior (within Rgc), or Galactic exterior (beyond Rgc).Upper center: Ẽ = 0, so the zero-velocity surface can expand no further before opening a neck between the cluster and Galactic interiors at L1.The closed portion of the surface passing through L1 is the tidal boundary.Upper right: Ẽ1 < Ẽ < Ẽ2 opens the neck at L1.The example trajectory-nearly periodic in the rotating frame-transits back-and-forth through the cluster/Galactic interiors and illustrates an escaper taking many tc to escape before eventually returning to the cluster.Lower left: Ẽ2 < Ẽ < Ẽ3 opens a second neck at L2, connecting all three domains.The example trajectory illustrates an escaper that immediately finds the neck to the Galactic interior and briefly passes again through the cluster after nearly 20 full Galactocentric orbits.Lower center: Ẽ3 < Ẽ < Ẽ4 opens a third neck at L3.The example trajectory transits between all three domains, with the escaper temporarily returning to the cluster after only four orbits about the Galaxy (to the cluster's three in that time).Lower right: Ẽ = Ẽ4 = 1, at which point the entire xy-plane is available to escapers, though the zero-velocity surface still exists out-of-plane.The example trajectory alternates between orbiting in the Galactic interior/exterior, passing back through the cluster during each transit.Though at higher-Ẽ and not nearly as regular, this trajectory is loosely analogous to those of some Jovian comets.

Figure 2 .
Figure2.Upper panel: The cumulative density function (CDF) for Ẽ at key steps transforming CMC's Ẽ (solid gray; approximate due to the assumption of spherical symmetry) into Gala's Ẽ (solid black; the true Ẽ in the full, asymmetric tidal field).Since Gala uses analyticallydefined potentials, deviation from a perfect analytic fit to the numerical CMC potential can dramatically bias the Ẽ fed into Gala.For illustrative purposes, the dotted red, yellow, and teal curves exemplify how such biases arise, with decreasing severity, from fitting Kepler, Plummer, and 2-component Plummer potentials, respectively.Even fitting a 3-component Plummer potential (our final choice; dashed light blue) still biases Ẽ to higher values.We undo this fitting-induced bias by minutely changing the velocities of each escaper via Equation (17), yielding the dashed dark blue curve-a perfect realignment of the post-fit Ẽ to match the original CMC distribution.Finally, the projection of radially-symmetric CMC positions and velocities into full 6-dimensional phase space (and precise numerical computation of Ecrit, rather than the expansion as in the CMC definition) smears out the low end of the Ẽ distribution when accounting for the asymmetry of the true tidal field in Gala (black curve).Lower panel: The final Ẽ CDF as computed in the Gala potential-the black curve from the upper panel-but shown for four different models: the archetypal NCC'd and CC'd GCs under both energy criteria-Equations (6) and (9).The α criterion results in significantly higher Ẽ.

Figure 3 .
Figure 3. Evolution of the GC models over time.Top left: number of particles relative to their initial number.Top right: cluster mass relative to the initial mass.Bottom left: core radius over half-mass radius (rolling average).Bottom right: 99% Lagrange radius over tidal radius (rolling average).In all panels, solid/dashed curves correspond to the NCC'd/CC'd GCs, black/blue to the raw/α escape criterion, and red to the near-identical simulations in W23 (which used the α criterion).

Figure 4
Figure 4 shows the distribution of Ẽ versus the clustercentric position r rmv /r t when each PE first satisfies Ẽ > 0. Corner plots show the cumulative density functions (CDFs) in Ẽ and r rmv /r t while the lower left and central panels show the corresponding scatter plots for the NCC'd and CC'd GCs (solid and dashed curves in the CDFs, respectively).Each use the raw criterion since the equivalent plot under the α criterion simply shifts PEs with Ẽ ≲ 0.1 to Ẽ ≳ 0.1.Colors distinguish escape mechanisms, indicated in the legends and caption.As a reproduction of Figure 2 in W23, with only slightly updated models and the vertical axis now Ẽ instead of velocity, we only highlight the Figure's key features.For more detailed discussion, including the algorithmic definitions of escape mechanisms, see W23.First, two-body relaxation (yellow) dominates overall, producing escapers with Ẽ ≪ 1.About half of escapers from relaxation originate within the typical core radius at removal r c (t rmv ), indicated by the vertical line and shaded interval in each scatter plot.This reflects that even bodies with Ẽ just below 0 in the GC halo typically first cross to Ẽ > 0 only after first plunging back through the core, where the higher density greatly enhances relaxation's efficiency (e.g.,Spitzer &

Figure 4 .
Figure 4. Escapers from the archetypal NCC'd (lower left panel) and CC'd (central panel) GC models, distributed according to position rrmv/rt and excess relative energy Ẽ upon removal from CMC (t = trmv).The corner plots show the corresponding CDFs for rrmv/rt and Ẽ, with solid (dashed) curves corresponding to the NCC'd (CC'd) model.The gray curves in the CDFs include all escapers while other colors distinguish different escape mechanisms.Regardless of mechanism, escapers (single or binary) containing a BH (black) or neutron star (but no BH; teal) are shown separately.All other escapers are categorized by mechanism: those caused by the induced kick from a binary companion's supernova (magenta), three-body binary formation (from three singles; red), binary-single (light blue) and binary-binary (dark blue) strong encounters, and two-body relaxation (yellow).The legends display the total number of escapers and subtotals for each category.The vertical lines and surrounding shaded intervals indicate the median and 10th-90th percentile range of the theoretical density-weighted core radius rc(trmv) from Casertano & Hut (1985), normalized by rt (trmv).

Figure 5 .
Figure 5. Same as Figure 4 but the horizontal axis is now the time of removal from CMC, trmv.

Figure 6 .
Figure6.Survival functions (1-CDF) for the escape timescale ∆tesc (left panels) and its normalized equivalent t = ∆tescω Ẽ2 (right panels).In the top two rows, we integrate each PE's trajectory for a full 14 Gyr beyond the time trmv that it first satisfies the raw (top) or α (2nd from top) escape criteria, assuming constant ϕc = ϕc(trmv).In the bottom two rows (again, one for each criterion), we instead integrate all PEs with trmv < 13 Gyr for only 1 Gyr beyond trmv in CMC's truly time-dependent ϕc(t)-see explanation in text.Again, solid curves denote the NCC'd GC and dashed the CC'd GC.The thick black curve shows the survival function across all escapers while the remaining curves show the contributions from each of 13 thin bins in Ẽ.These are colored in rainbow order (right to left in the left panels) with bin centers Ẽcen uniformly spaced in log-scale 1/4 dex apart Ẽcen = 10 −2.5 (dark red) to Ẽcen = 10 0.5 (violet).Each bin's lower/upper bound is defined narrowly as ( Ẽlow, Ẽupp) = (0.9, 1.1) × Ẽcen.Truncation of curves correspond to exclusion of PEs that do not cross beyond r CMC

Figure 7 .
Figure7.Escaper trajectories selected from the solid teal curve in the upper left panel of Figure6-Ẽ ∈ [0.9, 1.1] × 10 −1/2 -assuming a constant ϕc = ϕc(trmv) from the removal time trmv of each escaper.Each is projected into the xy-plane and integrated for time 1.1∆tesc, just beyond when they cross rt (blue circle; note the projection only makes it appear some of these do not cross rt ).Each panel shows 20 such escapers, one highlighted in red for clarity, belonging to a different ∆tesc interval (see labels), selected to be immediately before/between the first four plateaus visible in the ∆tesc distribution at this Ẽ.Trajectories immediately before the N th plateau typically loop back N − 1 times before escape.We show in blue the L1 (left) and L2 (right) Euler-Lagrange points and in gray the two forbidden realms (of the 80 across all shown trajectories) enclosing the GC the least (dark gray) and most (light gray).This spread is due primarily to the finite Ẽ bin width when selecting the displayed escapers and secondarily to the variation in µ(trmv) across these, since they are removed from CMC at different ages.

Figure 8 .
Figure 8. Same as the bottom two rows of Figure 6, but we instead integrate all PEs with trmv < 4 Gyr for 10 Gyr beyond trmv.

Figure 9 .
Figure 9. Same as the right panels in Figure 8 but with t redefined from Equation (18) for low Ẽ ≲ 0.03, the red curves.These converge much better to a common t distribution.

Figure 10 .
Figure 10.Same as Figure 8 but with t redefined from Equation (18) for high Ẽ ≳ 0.56, the blue and violet curves.These, too, converge much better to a common t distribution.

Figure 11
Figure 11 projects into the clustercentric xy-plane sample PE trajectories (from the raw criterion's NCC'd GC) integrated under constant ϕ c = ϕ c (t rmv ), distinct for each PE.Each panel shows 50 trajectories, one highlighted in red as a visual aid, belonging to different narrow bins spaced 1/2 dex apart in the range Ẽ ∈ [10 −2 , 10 2 ]-see caption.To give a sense of the timescale and avoid clutter, we cut off each trajectory 2 Gyr after it first crosses to r > r t .The PEs with low Ẽ (top left) escape through the necks in their forbidden realm (gray) near L1/L2, the points on the blue circle indicating r t (t).Increasing Ẽ expands the necks, allowing PEs to cross beyond r t with higher ẏ or ż (in/out of the page).The lower panels show how the forbidden realm's retreat from the xy-plane at Ẽ ≥ 1 enables more immediate return to the GC.At Ẽ ≈ 1 (lower left), in particular, the Coriolis effect causes ∼10% of escapers to temporarily return to r < r t -sometimes many

Figure 11 .
Figure11.As in Figure7, but the panels now distinguish subsets (of 50 escaper trajectories each) belonging to 8 bins in Ẽ ∈ [10 −2 , 10 2 ] from the archetypal NCC'd GC under the raw criterion and integrated to age trmv + ∆tesc + 2 Gyr.Except in the last panel (due to its high Ẽ), this time limit excludes the portions of the trajectories that return to the GC after entirely circumnavigating the Galaxy in the rotating frame, cleaning up the Figure considerably (see also Figure12, identical to this one but in the inertial center-of-mass frame).Unlike in Figure7, there is no additional filter for escapers belonging to specific windows in ∆tesc, so we make the bins ten times narrower-between 0.99 and 1.01 the Ẽ

Figure 12 .
Figure 12.As Figure 11, but shown in the rotating center-of-mass frame of the Galaxy and GC.

Figure 13 .
Figure13.Two-dimensionally projected positions of PEs, colored by time since removal from CMC, at age t = 12 Gyr (but see a full 14 Gyr movie accessible via DOI: 10.5281/zenodo.10714860and youtu.be/zJKCvAf6U3E).We show the denser CC'd GC under the α criterion simply to optimize visualization, as the higher density produces more PEs and the lower ∆tesc under the α criterion amplifies the color contrast.The lower left panel shows the projected positions in the static center-of-mass coordinates (Section 4.1), a face-on view of the GC orbit.The other panels show the three orthographic projections along each of the cardinal directions in the rotating clustercentric coordinates.The views in the upper row are edge-on to the GC's orbit, looking along (upper left) and perpendicular to (upper right) the ray connecting the cluster center to the Galactic center, while the view in the lower right panel is face-on to the cluster orbit.In the three orthographic panels, the blue circles have radii r/rt = [1, 2, 3, 4, 5].In the lower left panel, the red circles have radii R/kpc =[2, 4, 6, 8, 10, 12], and the blue circle radius r/rt = 5 with guiding center at R = Rgc.The true cluster center's slight offset from this circle's center illustrates the subtle decrease in the distance between the GC and center-of-mass as the GC loses mass.
cos θ and y ′ ≡ (1−µ)θ, where θ ≡ arctan2[y/(x+1−µ)].So in each Figure the GC's velocity points to the right and the upper panel is a face-on view to the orbit (with the MW's center down the page at x ′ = −1) while the lower panel is a panoramic edge-on view from the MW's center.The extreme sample size (≳10 8 individual stellar positions) achieved from time-averaging escaper trajectories from such large GC simulations makes these Figures the highest-resolution tidal tail/stellar stream density profiles we were able to find in the literature.

Figure 14 .
Figure 14.Projected surface number density Σ of PEs from the CC'd GC under the α criterion (see the full Figure set in the online journal for the other simulations).Σ is time-averaged between ages t/Gyr ∈ [11, 13], achieved by stacking 2001 snapshots of PEs in that interval, finely binning PEs by location, and dividing each bin count by 2001 and the bin area in pc 2 .As described in the text, the new clustercentric coordinates (x ′ , y ′ ), still in units of Rgc, are flattened to map the full circular orbit to the line segment y ′ ∈ [−π, π], where (x ′ , y ′ , z ′ ) = 0 corresponds to the GC center.So the GC's velocity is to the right.The upper panel is a face-on view to the GC orbit (with the Galactic center down the page at x ′ = −1) while the lower panel is a panoramic edge-on view from the Galactic center to any point along the GC orbit.

Figure 15 .
Figure 15.As Figure 14 (including being a Figure set in the online journal) but with the horizontal axis compressed to fully span the entire GC orbit.The first several epicyclic overdensities are visible in the lower panel, spaced of order 10rt apart.In the upper panel, the strip from x ′ Rgc/rt ∈ [−2, 2] has much lower density than in the tails/streams due to the presence of the forbidden realm there for low Ẽ (most escapers).
bottom), all from 501 uniform bins along the stellar stream's flattened y ′ -axis.Based on Figure 15, we define the stream at any time t as all PEs with |x ′ (t)| < 0.05 and |z(t)| < r t (t)/R gc , including for reference the PEs within the tidal boundary.The results are again timeaveraged across t/Gyr ∈ [11, 13] and the solid/dashed curves correspond to the NCC'd/CC'd GCs under the raw criterion.The blue/red curves distinguish the contributions from the leading/trailing tails to their combined profile (black).

Figure 16 .
Figure 16.The surface number density of PEs (viewed from the Galactic center; upper panel), their mean speed ⟨v − vc⟩ relative to the GC's circular speed (middle panel), and the dispersion in this speed σ ≡ ⟨(v − vc) 2 ⟩ − ⟨v − vc⟩ 2 (lower panel), along the stellar stream's flattened y ′ -axis.The profiles are time-averaged across t ∈ [11, 13] Gyr and split into 501 bins of uniform width ≈100 pc.To be counted in the tail/stream, a PE at time t must have |x ′ (t)| < 0.05 and |z(t)| < rt (t)/Rgc.As usual, solid/dashed curves indicate the NCC'd/CC'd GCs-in this case, under the raw escape criterion.The blue/red curves correspond to the leading/trailing tails and the black curves to their sum.

Figure 17 .
Figure17.As Figure16(including being a stack of 2001 snapshots spanning ages 11-13 Gyr) except for the following changes.To simulate the impact of stream disruption, we now exclude the PEs in each snapshot for which >500 Myr have elapsed since their 'escape' (first crossing beyond r > rt ).We also only show the sum of the leading/trailing tails (i.e., the black curves from Figure16) and do so for all four simulations.The solid/dashed curves still correspond to the NCC'd/CC'd GCs, and the blue/red curves to the raw/α escape criteria.Finally, because of the limited tail length from the time cutoff, we use finer bins in y ′ -2001 across the entire GC orbit, corresponding to bin widths ≈25 pc.
's black curves (both tails combined) in Figure 17, this time showing all four simulations and eliminating the contribution of the returning tails.We do so by excluding from each snapshot in the Figure's 11-13 Gyr stack any escapers for which >500 Myr have elapsed since their first crossing to r > r t .This cutoff produces streams with angular span in Galactic longitude of ≈40 deg, roughly average for streams associated with MWGCs (e.g.Mateu 2023)-hence the truncated horizontal axis.
Figure19further excludes the returning tails by cutting out PEs 500 Myr after their first 'escape' (crossing to r > r t ).Doing so clearly increases ⟨v − v c ⟩ from the perspective of the Galactic center (lower panels) and reduces the dispersion in ⟨v − v c ⟩ (spread in color) at any y ′ because there is no longer a returning tail to oppose the local flow of the outgoing tail.

Figure 18 .
Figure 18.As Figure 14, but mapping the time-averaged local mean stream speed ⟨v − vc⟩ relative to the GC circular speed instead of the surface density, and counting only PEs with |x ′ (t)| < 0.05 and z(t) < rt (t)/Rgc.Black regions indicate bins outside this interval or where no escapers are present in any of the 2001 snapshots across ages t ∈ [11, 13] Gyr.

Figure 19 .
Figure 19.As Figure 18, but excluding the returning tails by cutting out PEs 500 Myr after their first 'escape' (crossing to r > rt ).