A supra-massive population of stellar-mass black holes in the globular cluster Palomar 5
Mark Gieles, Denis Erkal, Fabio Antonini, Eduardo Balbinot, Jorge Peñarrubia
AA supra massive population of stellar-massblack holes in the globular cluster Palomar 5
M. Gieles , , D. Erkal , F. Antonini , Eduardo Balbinot and Jorge Pe˜narrubia , ICREA, Pg. Llu´ıs Companys 23, E08010 Barcelona, Spain Institut de Ci`encies del Cosmos (ICCUB), Universitat de Barcelona (IEEC-UB), Mart´ı i Franqu`es1, E08028 Barcelona, Spain Department of Physics, University of Surrey, Guildford, GU2 7XH, Surrey, UK Gravity Exploration Institute, School of Physics and Astronomy, Cardiff University, Cardiff,CF24 3AA, UK Kapteyn Astronomical Institute, University of Groningen, Postbus 800, NL-9700AV Groningen,The Netherlands Institute for Astronomy, University of Edinburgh, Royal Observatory, Blackford Hill, EdinburghEH9 3HJ, UK Centre for Statistics, University of Edinburgh, School of Mathematics, Edinburgh EH9 3FD, UK
Palomar 5 is one of the sparsest star clusters in the Galactic halo and is best-known for itsspectacular tidal tails, spanning over 20 degrees across the sky. With N -body simulations weshow that both distinguishing features can result from a stellar-mass black hole population,comprising ∼ of the present-day cluster mass. In this scenario, Palomar 5 formed witha ‘normal’ black hole mass fraction of a few per cent, but stars were lost at a higher rate thanblack holes, such that the black hole fraction gradually increased. This inflated the cluster,enhancing tidal stripping and tail formation. A gigayear from now, the cluster will dissolveas a 100% black hole cluster. Initially denser clusters end up with lower black hole fractions,smaller sizes, and no observable tails at present day. Black hole-dominated, extended starclusters are therefore the likely progenitors of the recently discovered thin stellar streams inthe Galactic halo. In recent years, a few dozen thin ( (cid:46) pc) stellar tidal streams have been discovered in theMilky Way halo . Their elemental abundances and distribution in the Galaxy provide importantconstraints on the hierarchical growth of the Milky Way and its dark matter distribution . Theirnarrow widths imply that their progenitor stellar systems had a low velocity dispersion and weredark matter-free star clusters rather than dark matter-dominated dwarf galaxies. However, for noneof these streams has a progenitor been found and the star cluster nature is questioned by two recentfindings: firstly, the inferred mass-loss rate of the GD-1 stream is several times higher than whatis found in ‘vanilla’ cluster evolution models and secondly, only mild tidal distortions and no tidaltidal tails were found for several globular clusters (GCs) with extremely radial orbits passingthrough the strong tidal field near the Galactic center where tidal stripping is efficient. Theseresults raise the questions what was driving the escape rate of the streams’ progenitors and whythis mechanism is not operating in all star clusters.The metal-poor GC Palomar 5 (hereafter Pal 5) is one of the few known star cluster with1 a r X i v : . [ a s t r o - ph . GA ] F e b xtended tidal tails associated with it , spanning (cid:38) degrees on the sky, making it a Rosettastone for understanding tidal tail/stream formation. The cluster has an unusually large half-lightradius of ∼ pc and combined with its relatively low mass of ∼ M (cid:12) , its average density isamong the lowest of all Milky Way GCs: ∼ . M (cid:12) / pc , comparable to the stellar density in thesolar neighbourhood. A low density facilitates tidal stripping and the formation of tidal tails , butit is not known whether this low density is the result of nature, or nurture.It has been proposed that Pal 5 simply formed with a low density and has always beena collisionless system, meaning that two-body interactions were not important in its evolution.However, some properties of Pal 5 are reminiscent of other GCs, such as a spread in sodiumabundances and a flat stellar mass function . These features have been attributed to high initialdensities
15, 16 and collisional evolution and suggest that Pal 5 in fact is, or was, a collisional stellarsystem, like the rest of the Milky Way GCs . In this study we aim to reconcile the low density ofPal 5 with collisional evolution.Since the discovery of gravitational waves , model predictions for the mass function ofstellar-mass black holes (BHs) and their natal kicks have changed . It is now believed that alarge fraction of BHs that form from massive stars have masses above M (cid:12) and do not receivea natal kick, as the result of fallback of material, damping the momentum kick resulting fromasymmetries in the supernova explosion . Models of cluster evolution show that the presence of aBH population accelerates the relaxation driven expansion
20, 21 and escape rate
22, 23 . Observationalmotivation for considering the effect of BHs on GC evolution stems from the discovery of accretingBH candidates in several GCs with deep radio observations
24, 25 and a BH candidate in a detachedbinary was reported in NGC 3201 . Given these recent insights, we here ask whether Pal 5 wasonce much denser, and that the present-day structure and prominent tidal tails are the result of aBH population. Results
We perform star-by-star, gravitational N -body simulations with NBODY
GPU . All clustersare evolved for 11.5 Gyr on the orbit of Pal 5 in a three-component Milky Way (bulge, disc, halo)and the simulations include the effect of stellar and binary evolution. No primordial binaries wereincluded. We consider two prescriptions for BH natal kicks. First, we consider the most up-to-dateBH recipes , in which approximately 73% of the mass of the BH population is retained after natalkicks, almost independently of the initial escape velocity of the cluster (see Methods). Then, inthe second suite of models we test the collisionless hypothesis and draw BH kick velocities froma Maxwellian with the same dispersion as that of neutron stars (NSs). Because we need lowerdensities for these models, all BHs are ejected by the natal kicks in almost all models (in fourmodels a single BH was retained). The latter approach is similar to that of Dehnen et al. , withthe added effect of stellar evolution and using a direct integrator (i.e. without force softening),such that the effect of two-body interactions are included. Although the initial two-body relaxationtimescale of these models is long ( ∼ Gyr), two-body relaxation could be important in the latestages of evolution when the cluster approaches dissolution. We will be referring to these two setsof models as wBH and noBH, respectively. The initial parameters and the results of both sets ofsimulations are summarised in Tables 1 and 2.We vary the initial number of stars ( N ) and the initial mass density within the half-mass2adius: ρ h0 ≡ M / (8 πr ) , where M (cid:39) . M (cid:12) × N is the initial mass and r h0 is the initialhalf-mass radius. We search for a model that best reproduces the observed number of stars ofPal 5, N cluster = 1550 , and its half-light radius, R eff = 3 . arcmin, which equals . pc at theadopted distance of 19.98 kpc (see Methods). R eff is defined as the distance to the cluster centercontaining half the number of observed stars. We first run a coarse grid of models, followed by afiner grid close to the parameters that give the best match in the coarse grid. The wBH model thatbest reproduces N cluster and R eff (wBH-1) has N = 2 . × and ρ h0 = 80 M (cid:12) / pc . This clusterlost 92% of its initial mass of . × M (cid:12) by stellar evolution and escapers and the densitydecreased nearly three orders of magnitude because of stellar evolution and dynamical heating byBHs. We note that the half-mass radius r h (cid:39) . pc, as determined from the 3-dimensional massdistribution, is similar to R eff = 18 . pc, which is determined from the projected distribution ofobservable stars. If mass follows light, R eff (cid:39) . r h and for mass segregated clusters withoutBHs it can be as small as . r h . The fact that R eff (cid:39) r h for wBH-1 implies that the observablestars are less concentrated than the mass profile. This is the result of the BH population whichsinks to the center via dynamical friction against the lower mass stars
20, 21, 30 , where they remain ina quasi-equilibrium density distribution that is more centrally concentrated than the stars . Thesurface density profile of the observable stars in the cluster and the properties of the stream arein good agreement with the observations (Figure 1). The observed over and under-densities in thetails are not reproduced by our model (see Figure 1 in the Supplementary Material). This is becausethey are likely the result of interactions with dark matter subhalo or the Galactic bar and giantmolecular clouds in the disc , which are not included in our model. The observed line-of-sightvelocity dispersion of stars in the tails is . ± . km/s , which is well reproduced by wBH-1( . ± . km/s, for giants in the same region on the sky as the observations). The most strikingproperty of wBH-1 is its large BH fraction at present: f BH = 22% . We define f BH as the total massin BHs that are bound to the cluster over the total bound cluster mass, i.e. f BH = M BH /M . TheBH population is made up of 124 BHs with an average mass of . M (cid:12) (i.e. M BH = 2178 M (cid:12) ),currently residing within R eff (see Figure 1). This f BH is more than twice as large as what isexpected from a canonical stellar initial mass function (IMF) and stellar evolution alone, and isthe result of the efficient loss of stars over the tidal boundary, while the BHs were mostly retainedbecause they are in the center.In all wBH models, f BH increases in the first ∼ Myr because of BH formation and stellarevolution mass loss. Then, f BH decreases in the following ∼ Myr because of BH ejectionsfrom the core . This happens because binary BHs form and interact with other BHs when two-body relaxation becomes important. These binaries become more bound in these interactions,eventually ejecting BHs and themselves from the cluster. What happens next depends on the initialdensity ρ h0 : in dense clusters, most BHs are ejected before the cluster dissolves and these GCs haveescape rates and R eff comparable to what is found in models of clusters without BHs (Figure 3).Clusters with lower initial densities have larger relaxation times, resulting in fewer BH ejections,while tidal stripping of stars is efficient, leading to an increasing f BH until (Figure 3).When clusters have f BH (cid:38) . , their R eff are ∼ − pc, comparable to what is found for abouthalf of the GCs in the outer halo of the Milky Way , suggesting that these ‘fluffy’ GCs are BH richand candidates to produce prominent stellar streams. This idea is further supported by the strong3orrelation between R eff and f BH and the fraction of stars in the stream with f BH (Figures 2 and 4).Theory suggests that in idealised single-mass star clusters that fill their tidal radius there existsa critical f BH (cid:39) at which the mass-loss rate of stars and BHs is the same and f BH remainsconstant at while the cluster loses mass. For higher f BH , stellar mass is lost at a higher rateby tidal striping than BH mass is lost by ejections from the core, such that f BH increases, and viceversa for f BH < . This implies that clusters can evolve to 100% BH clusters if they form with f BH (cid:38) and can remain above it during their evolution. Our models suggest that for multimassmodels this critical fraction is lower: ∼ . . Since for a canonical IMF the initial f BH (cid:39) − ,depending on metallicity, it is possible for some GCs to remain above the critical f BH and evolveto 100% BH clusters , as we find in our wBH-1 model. Because f BH always evolves away from10%, the distribution of f BH values of a GC population becomes bimodal. Whether cluster evolvetowards BH-free clusters, or 100% BH clusters depends on their initial density relative to the tidaldensity. Because f BH affects the density, a unimodal initial density distribution can evolve towardsa bimodal preseny-day density distribution, as is observed in the halo .We now discuss the noBH models. The best-fit model (noBH-1) has N = 3 . × and ρ h0 = 9 . M (cid:12) / pc , i.e. approximately twice as massive and an order of magnitude less densethan wBH-1. The resulting observable parameters N cluster ( R eff ) are within of the valuesof Pal 5. The resulting cluster density profile and stream properties are similar to those of wBH-1(Figure 1 in Supplementary material), hence based on these observables it is not possible to prefereither of the assumptions of the BH kicks. The wBH-1 model predicts a higher central velocitydispersion of 580 m/s vs. 350 m/s for noBH-1 (see Methods). The inferred central dispersionfrom the literature compilation by Baumgardt & Hilker is +150 − m/s, i.e. favouring the BHhypothesis. We quantified the degree of fine-tuning of N and ρ h0 that is required to obtain the bestmodel in both cases. For the wBH models, we find that the uncertainty in the present-day propertiescan be covered by a relatively large range of initial densities, while for the noBH models we findthat variations in the initial density are amplified by more than an order of magnitude in variationsin the present-day properties (see Methods). This means that a relatively large range of initialconditions in wBH models lead to similar present-day properties, while for the noBH models ahigh degree of fine-tuning in the initial density is required to obtain noBH-1. In addition, noBH-1completely dissolves at 11.8 Gyr, i.e. 300 Myr after we observe the cluster, while wBH-1 survivesfor another 1.4 Gyr. By assuming simple power-law distributions for the initial cluster massesand densities, we estimate that the probability of finding a cluster with the properties of wBH-1(noBH-1) is / . × ) . Given that the Milky Way has ∼ GCs, the wBH-1 modelprovides the more likely explanation for Pal 5. The velocity dispersion, the fine-tuning argumentsand the timing arguments all favour the BH hypothesis. The final argument in support of the BHhypothesis is that the higher initial density and the resulting collisional nature of Pal 5 make iteasier to understand the flat stellar mass function, its multiple populations and its relation to therest of the Milky Way GC population.In the observed mass range ( . − . M (cid:12) ), the mass function of Pal 5 is flatter ( dN/dm ∝ m α , with α (cid:39) − . ) than what is expected from a canonical IMF ( α = − . ), suggesting that Pal5 has preferentially lost low-mass stars. We note that this result has been questioned and the massfunction slope may actually be close to the initial value (Iskren Georgiev, private communication).4n Figure 5 we show the mass function of bound stars and remnants of wBH-1. For main sequencestars > . M (cid:12) it has a slope of α = − . +0 . − . , which is flatter than the IMF, but steeper thanthe observed slope. The noBH-1 model has a slope comparable to wBH-1, implying that in thelater parts of its evolution the cluster became collisional. The similarity in mass function slopebetween wBH-1 and noBH-1 means that we can not use the observed mass function to distinguishbetween the two scenarios. However, in the noBH models the only way to reconcile the modelswith the observations is to start with a flatter IMF, because the alternative (i.e. reducing the initialrelaxation timescale by increasing the initial density) leads to a R eff that is too small. The BHscenario, however, leaves the possibility to increase the initial cluster density and still end up withthe same present-day BH fraction, R eff and tidal tails. For example, tidal heating by interstellargas clouds in the early evolution is probably important in the evolution of GCs
38, 39 , but this effectis not included in our models. This tidal heating leads to a decrease in the cluster mass
40, 41 andfor clusters with a low central concentration the cluster density also decreases . BHs sink to thecluster center by dynamical friction on a fraction of a relaxation time and if this is shorter thanthe tidal heating timescale, then the BHs are less affected by this mechanism of tidal strippingthan the stars, such that f BH increases, counteracting the reduction of f BH from BH ejectionsfrom the core due to the higher density. This would be a pathway to start with higher densities.Alternatively, something may have happened to Pal 5 at a later stage of its evolution. It is likelythat Pal 5 formed in a dwarf galaxy that was accreted onto the Milky Way, which is supportedby its tentative association with the Helmi streams . The cluster could have lost a fraction of itsloosely bound stars as a result of the removal of its host galaxy , thereby increasing f BH . Finally,a flatter IMF at high masses would also lead to a higher f BH . This also affects the evolution, anda flatter IMF would lead to faster expansion higher mass loss rate
22, 23, 45 . In all these scenarios,denser initial conditions and therefore more equipartition among low-mass stars, whilst endingwith the same f BH at present, is a possibility. Understanding the interplay between early massloss, mass segregation, the (high-mass) IMF and accretion on the Milky Way to find limits on theallowed initial density is an interesting topic for a follow-up study. The initial density is a criticalingredient for our understanding of GC formation
15, 16, 46 , evolution , BH natal kicks and binaryBH mergers .Our results have implications for our understanding of GC evolution. In ‘vanilla’ modelsof cluster evolution without BHs, the relaxation driven escape rate on the orbit of Pal 5 is about ∼ M (cid:12) / Myr (Figure 3). We note that our noBH models have much higher escape rates, but thisis because they start with lower initial densities than what is usually done . It is well-establishedthat the mass loss rate in ‘vanilla’ models is insufficient to ‘turn-over’ a power-law initial clustermass function with index − , as is observed for young star clusters in the nearby Universe , intothe observed peaked (logarithmic) mass distribution with a typical mass of × M (cid:12) : too manylow-mass GCs survive beyond (cid:38) kpc . The required escape rate is about an order magnitudelarger: M (cid:12) / Myr . From our models we can conclude that Pal 5 is currently losing massat that rate (Figure 3), hence if a large fraction of GCs go through a similar evolutionary phase,then relaxation driven evaporation is more important in shaping the GC mass function than usuallyassumed, reducing the need for additional GC disruption mechanisms
38, 39 or a peaked initial clustermass function. Because variations in f BH lead to an order of magnitude variation in the escape rate5Figure 3), we conclude that the effect of BHs is of comparable importance in shaping the GC massfunction as the details of the Galactic orbit. Discussion
We now consider our results in the context of the Milky Way GC population. About of theMilky Way GCs have R eff (cid:38) pc, which led some authors to suggest that two-body relaxationis not important and that their evolution is collisionless
12, 53 . These fluffy GCs are found predom-inantly at large Galactocentric radii and have low masses ( (cid:46) M (cid:12) ). More specifically, beyond (cid:38) kpc from the Galactic center, half of the GCs have R eff = 10 − pc . Our results suggest thatthese large and low-density GCs are BH rich. Their low densities, relative to the tidal densities ontheir orbits can be used to constrain their initial densities and masses in a similar way as we did forPal 5. In wBH-1, about half of the observable stars are lost in the last 3 Gyr and in this period thestream became visible above the background. Combined with the remaining lifetime of ∼ Gyr,Pal 5 has an observable stream for 30% of its lifetime. There are low-density GCs beyond 8 kpcfrom the Galactic center and for about half of them tidal tails or tidal features have been found(Table 1 in the Supplementary Material). Of these, there are 12 that are at distances comparable toPal 5, or closer. If we assume that the escape rate is the same for these GCs, combined with thefact the GC mass function is uniform at low masses, then about 30% of these GCs are in the final30% of their evolution. From this we expect that ∼ GCs to have prominent tidal tails, providingan explanation for the rarity of Pal 5 and the low-number of known GCs with streams. We cannow make an order of magnitude estimate of how many streams the MW fluffy GC population hasgenerated in the past. The ∼ fluffy, nearby, GCs all have masses below M (cid:12) , hence the massfunction of fluffy clusters is dN/dM (cid:39) − M − (cid:12) . If these GCs all lose mass at a constant rate of M (cid:12) / Myr , then this mass function is constant in time and we expect that these GCs contributed × − streams per Myr, or about 20 streams in the last 10 Gyr. This estimate supports theidea that the ∼ known cold streams in the halo resulted from BH-rich, extended GCs thatdissolved in the Milky Way halo.We conclude with a discussion on observational tests of the BH hypothesis. We computed therate of microlensing events of background quasars in the field of Pal 5 and find that the event rateis too low ( ∼ − / yr). We then looked at the effect of BHs on stellar kinematics and found that itis possible to infer the BHs from the kinematics of the stars. For a BH fraction of 22%, the centralvelocity dispersion of giant stars within R eff is m/s, which is ∼ m/s higher than for acluster without BHs (Figure 6). The significance of the available velocity dispersion measurementof +150 − m/s can be improved by increasing the sample of stars, and constraining the propertiesof binary stars, including their orbital elements. We estimate an upper limit for the periods ofbinaries with giant stars of P (cid:39) . yr because longer period binaries were ionised by interactionswith BHs. Finding a binary with a larger period would challenge our BH hypothesis, while theabsence of such binaries could in combination with additional N -body model with primordialbinaries be used as support for it. Finally, BHs with stellar binary companions may also exist andcould be found from their large velocity variations. A multi-epoch observing campaign to obtainline-of-sight velocities of Pal 5 is can therefore be used to establish with high precision the centralvelocity dispersion and the properties of the binaries (see Methods for details). This would providethe critical test of the hypothesis that Pal 5 hosts a (cid:38) × M (cid:12) population of stellar-mass BHs.6 D E C [ d e g ] observable starsblack holes (b) Pal 5 stream R eff log projected radius [arcmin] l o g nu m b e r d e n s i t y [ a r c m i n ] (a) Pal 5 cluster N -body (wBH-1)Observations (Ibata et al. 2017)01log density [arcmin ] (c) (d) Observations (Erkal et al. 2017) N -body (wBH-1) Figure 1: Comparison between the N -body model wBH-1 and observations of Pal 5 and its stream:(a) Density profile of observable stars in wBH-1, which provides an excellent match to the obser-vations. There is also good agreement between the stellar stream track (b), the stream density andits gradient (c) and stream width (d) between observations and wBH-1, implying that the rate of es-cape of stars from Pal 5 and their velocity dispersion in the last few Gyrs are correctly reproducedby wBH-1. The blowout of the cluster shows that almost all of the 124 bound BHs are within R eff .7 R e ff [ p c ] Pal 5 (Ibata et al. 2017) f BH N t r a ili n g / N c l u s t e r Pal 5 (Erkal et al. 2017)
Figure 2: Effective radius ( R eff , top) and the relative number of stars in the trailing tail (bottom)for all surviving wBH clusters when their mass is . × M (cid:12) (the mass of wBH-1 at 11.5 Gyr)as a function of f BH . Observed values
11, 32 are shown as horizontal lines with the shaded regionindicating the 1 σ confidence interval. Values of wBH-1 are shown with large stars. Here N trailing is defined as the number of observable stars in the trailing tail, within 4 degree of right ascensionof Pal 5, which in the model corresponds to stars that escaped in the preceding ∼ Gyr. Becausethese quantities vary along the orbit at these low masses, we plot with error bars the range withinthe nearest peri and apo. The strong correlation of both quantities with f BH suggests that for agiven orbit and cluster mass, the BH content sets both R eff of the cluster and the prominence of thetidal tails. 8 l o g ( M ) [ M M y r ] c o lli s i o n a l N - b o d y w i t h o u t B H s B a u m g a r d t & M a k i n o ( ) s t e ll a r e v o l u t i o n collisionless N -body Dehnen et al. (2004) R e ff [ p c ] h0 R e ff M / M [ M ]210 l o g f B H M B H = M Figure 3: Results of wBH-1 (blue, thick lines), compared to other wBH models with N = 100 k (green lines) and k (red lines) with different initial densities. Different panels show the variationof diagnostic quantities as a function of bound cluster mass ( x -axis), which decreases in time. Top:Total (positive) mass-loss rate ( − ˙ M ), which is initially high as the result of stellar evolution. Afterabout 30% of the mass is lost, ˙ M is dominated by escaping stars and BHs. The results of thecollisionless N -body model of Pal 5 and that of the collisional N -body models without BHs are overplotted. Because of the BHs, wBH-1 loses mass at a similar rate as the lower densitycollisionless models. Middle: Evolution of R eff with the age of 11.5 Gyr indicated with dots. Themodels do not evolve towards a constant luminosity density (i.e. R eff ∝ M / ), as predicted forsingle-mass clusters . Instead, R eff remains approximately constant while the cluster evolves tolower mass and there is a factor of five spread in R eff at a given mass due to variations in f BH (bottom). Dense clusters eject all BHs, while low-density clusters lose stellar mass at a higher ratethan BH mass, and evolve along tracks of nearly constant BH mass towards a 100% BH cluster.9 b [ d e g ] wBH-10 f BH = 0.32%40.042.545.047.5 b [ d e g ] wBH-5 f BH = 3.1%40.042.545.047.5 b [ d e g ] wBH-7 f BH = 15% 15105051015 l [deg]40.042.545.047.5 b [ d e g ] wBH-1 f BH = 22% Figure 4: Density map in Galactic coordinates of unbound stars of four different wBH modelswith different black hole fractions at 11.5 Gyr. A density of background stars of . / arcmin wasadded to each model. Models with f BH (cid:46) do not have prominent tidal tails.10 .0 0.5 0.0 0.5 1.0 1.5 2.0log m / M l o g M j / M i n i t i a l m a ss f un c t i o n stars whitedwarfs neutronstars blackholesstellar evolutionstellar evolution + escaperssupernovakicks dynamicalejections Figure 5: Mass function of bound stars and remnants of wBH-1 at 11.5 Gyr. The mass in eachbin ( M j ) is plotted for the IMF (red line), the mass function after only stellar evolution has beenapplied (light shaded), and for the final results of the N -body simulation (dark shaded). There is aflattening of the stellar mass function due to preferential escape of low-mass stars.11 l o g / a r c m i n (a)wBH-1 (b)noBH-1King (1966) N -body1.5 1.0 0.5 0.0 0.5 1.0log R / R eff L O S [ k m / s ] (c) 1.5 1.0 0.5 0.0 0.5 1.0log R / R eff (d) Figure 6: Fits of King models plus a constant background (green lines and shaded region repre-senting the σ uncertainty) to the surface density profile of observable stars (blue circles with errorbars) for wBH-1 (a) and noBH-1 (b). Panels (c) and (d) show the resulting velocity dispersion ofthe model fits scaled to a total mass of × M (cid:12) and half-mass radius of pc and assuming thatthe total mass is due to all the stars, white dwarfs and neutron stars, i.e. no BHs. Blue points witherror bars in (c) and (d) are the dispersions in the N -body models of the giants within R eff scaledto the same mass and half-mass radius for wBH-1 and noBH-1, respectively. For noBH-1 the dis-persion of the giants is . ± . , consistent with the King model prediction. For wBH-1, thisdispersion is . ± . km/s, or 3 σ above the central dispersion prediction from the King model.With precise line-of-sight velocities and a photometric mass estimate this difference is detectable.12 ethodsMilky Way model and orbit of Pal 5. We adopt a time-independent, axisymmetric Milky Way,consisting of a dark matter (DM) halo, a disc and a bulge. We derive the Milky Way parametersand the orbit of Pal 5 from a fit to the Pal 5 stream. This stream fit is nearly identical to that inErkal et al. which used the MWPotential2014 potential of Bovy , but we use updated priorsbased on more recent measurements of distance and proper motion of Pal 5 and the distance tothe Galactic center . The fit treats Pal 5’s present-day distance, radial velocity, proper motions,the distance to the Galactic center, and the DM-halo mass and scale radius as free parameters.The DM halo is described by a spherical NFW profile , with potential Φ halo ( R G ) = − GM NFW R G ln (cid:18) R G a NFW (cid:19) , (1)where R G is the Galactocentric radius, M NFW = 4 . × M (cid:12) is a mass scale and a NFW =16 . kpc is the scale radius. For these parameters the virial mass and concentration are M vir =8 . × M (cid:12) and c = 15 . , respectively. We note that these parameters are close to those ofBovy . The disc potential is that of an axisymmetric Miyamoto-Nagai disc Φ disc ( R G ) = − GM MN (cid:104) X + Y + (cid:0) a MN + √ Z + b (cid:1) (cid:105) / , (2)where M MN = 6 . × M (cid:12) is the mass of the disc, a MN = 3 kpc is the disc scale lengthand b = 0 . kpc is the scale height and X, Y, Z are Cartesian Galactocentric coordinates (i.e. R = X + Y + Z ). The bulge is described by a Hernquist potential Φ bulge ( R G ) = − GM H R G + a H , (3)with M H = 5 × M (cid:12) the bulge mass and a H = 0 . kpc its length scale. The bulge modelis slightly different from what was used by Erkal et al. , but within the orbit of Pal 5 the totalenclosed mass of both bulge models is the same.Cartesian Galactocentric phase-space coordinates of Pal 5 in the potential described abovewere given by the fit: [5 . , . , . kpc and [ − . , − . , − . km/s. The Sun isfound to be at [ − . , , kpc, with a velocity of [11 . , . , . km/s, which puts Pal 5 at adistance from the Sun of 19.98 kpc. To find the initial position and velocity of Pal 5’s orbit, weflipped the sign of the velocity vector and integrated the orbit backward in time. We adopt an ageof 11.5 Gyr and fix this for all models. The typical uncertainty on the age determinations is ∼ Gyr and varying the age would somewhat increase the uncertainty on our initial parameters,but not the parameters of our best model. The initial Galactocentric position and velocity of Pal 5are [1 . , . , . kpc and [43 . , − . , . km/s, respectively. Cluster initial conditions and N -body code. The initial positions and velocities of the stars weresampled from an isotropic Plummer model , truncated at 20 scale radii. Initial stellar masses weresampled from a Kroupa IMF in the range . − M (cid:12) and a metallicity of Z = 6 × − wasadopted, i.e. [Fe / H] (cid:39) − . . All simulations were run with the direct (i.e. no softening) N -body13ode NBODY
GPU
27, 66, 67 , which deploys a 4 th -order Hermite integrator with an Ahmad-Cohenneighbour scheme
68, 69 . It has recipes for stellar and binary evolution
70, 71 , with recent updatesfor BH masses and kicks and it deals with close encounters with Kustaanheimo-Stiefel (KS)regularisation. We use the Graphics Processing Unit (GPU)-enabled version and the simulationswere run on a server with GeForce RTX 2080 Ti GPUs at ICCUB. A few modifications to the codewere made for this project. The isothermal halo was replaced by the NFW halo, with the force andforce derivatives derived from equation (1). Stars were stripped from the simulation when theyreached a distance of 40 times the instantaneous half-mass radius of the bound particles and foreach escaper the time, position and velocity in the Galactic frame were stored. The escapers werethen integrated as tracer particles in the Galactic potential until 11.5 Gyr with a separate integratorto construct the tidal tails. We did not include the contribution of the cluster to the equation ofmotion of the tail stars. Near the end of the simulation the fractional contribution of the cluster tothe total force is approximately × − for stars that just left the cluster, and smaller for stars atlarger distances, justifying this assumption. All models were evolved until complete dissolution,and for the models that dissolved after 11.5 Gyr, a snapshot was saved at exactly 11.5 Gyr, i.e.when the cluster is at the position of Pal 5 today. Black hole recipes.
The fast evolution codes
SSE for stars and BSE for binaries in NBODY
GPU were recently modified with updated prescriptions for wind mass loss, compact object (i.e. neu-tron stars and BHs) formation and supernova kicks . We adopt here the rapid supernova mecha-nism in which the explosion is assumed to occur within the first 250ms after bounce , which cor-responds to nsflag=3 in SSE / BSE . The BH natal kick velocities are drawn from a Maxwellianwith dispersion σ = 265 km/s , and in the wBH models they are subsequently lowered by theamount of fallback such that momentum is conserved (i.e. kmech = 1 ). As a result, 63%(73%)of the number(mass) of BHs does not receive a kick for the IMF and metallicity we used. TheBHs that receive a kick form from stars with ZAMS masses in the range − M (cid:12) , resultingin BH masses in the range − M (cid:12) , and because of the low escape velocities of our modelclusters (10-20 km/s) they are almost all lost. As a result, the lowest mass BH in the cluster justafter supernovae and kicks has a mass of M (cid:12) and the average BH mass is M (cid:12) .In natal kick models that consider the effect of fallback of mass on the BH , the exact frac-tion of BHs that do not receive a kick depends on the prescription for compact-object formation.For our IMF and metallicity we find for the S TAR T RACK , rapid, and delayed explosion mecha-nisms described in section 4 of Fryer et al. , that about , and of the total BH mass isin BHs that form without a natal kick, respectively. Our result of is therefore an intermediatevalue. A lower(higher) initial f BH requires a lower(higher) initial cluster density to end up with thesame cluster properties at the present day. Pal 5 parameters.
To find the number of observed stars ( N cluster ) and half-light radius ( R eff ) of Pal5, we fit King models to the density profile of figure 10 of Ibata et al. using the LIMEPY models and the Markov Chain Monte Carlo (MCMC) code EMCEE . We fit for the dimensionless centralpotential: W = 4 . ± . , the number of observed stars: N cluster = 1550 ± , the half-massradius: r h = 4 . ± . arcmin and the background: . ± . arcmin − . For these parameters R eff = 3 . ± . arcmin, which for the adopted distance to Pal 5 of 19.98 kpc corresponds to R eff = 18 . ± . pc. 14 ata analysis. For each simulation, snapshots were saved approximately every 50 Myr. For eachsnapshot, we find the stars and remnants that are energetically bound to the cluster, i.e. which havea specific energy E i = 0 . v i + φ i < , where v i is the velocity in the cluster’s centre of mass frameand φ i is the potential due to the mass of the other bound stars, which we determine iteratively.From the bound particles we determine the total cluster mass M , the mass of the BH population M BH , the half-mass radius r h , N cluster and R eff .To compare the N -body model to observations, we use isochrones to convert masses of starsin different evolutionary phases to magnitudes. We use SDSS g -band magnitudes from MISTisochrones
78, 79 for 11.5 Gyr, Z = 6 × − ( [Fe / H] = − . ), [ α/ Fe] = 0 and rotational velocitiesof 0.4 times critical. For the observations shown in Figure 1, CFHT photometry was used whichis a slightly different photometric system than SDSS. A colour transformation is provided by Ibataet al. . For the magnitude limits we adopt here, the colors of Pal 5 stars vary between . (cid:46) ( g − r ) CHFT (cid:46) . in the relevant magnitude range, resulting in corrections in the g -band between . − . . We therefore adopt g SDSS = g CFHT + 0 . to convert between the two systems. For theobserved number density profile, a magnitude range of < g CFHT < was used by Ibata et al. .We apply the corresponding magnitude range . < g SDSS < . and combined with the adopteddistance (DM=16.5 mag) this implies a mass range of main sequence stars of . − . M (cid:12) ,which is what we use to select stars in the N -body model to construct the surface density profile.For the stream we selected stars with . < g CFHT < . mag (i.e. . < g SDSS < . ) asin Ibata et al. . This magnitude cut implies mass limits that depend on the distance of stars in thetidal tails. Finding the best model.
We varied the initial number of stars N and the initial half-mass density ρ h0 and try to reproduce N cluster and R eff . For the wBH models, we first ran a coarse grid of modelswith N = [50 k, k, k ] and ρ h0 = [30 , , , ] M (cid:12) / pc . We ran each combination of N and ρ h0 , apart from N = 200 k and ρ h0 = 10 M (cid:12) / pc . For each model we computed thefractional differences between the N -body results and the observations for N cluster and R eff at 11.5Gyr: δ N and δ R , respectively, and then define the best model as the one for which δ = (cid:112) δ N + δ R islowest. For each model we also find the time to reach N cluster = 1550 ( T ). We introduce T toestablish a goodness of fit for models that dissolve before 11.5 Gyr. Model IDs are increasing withthe value of δ and for clusters that dissolved the ID is in order of increasing difference between T and 11.5 Gyr. The model with the smallest δ in the coarse grid has N = 200 k, ρ h0 =100 M (cid:12) / pc and reproduces the two observables within 10%. From the coarse grid we estimatethat a slightly larger N and lower ρ h0 would reduce the difference. We then ran a finer gridwith 6 more models with N in the range k − k and ρ h0 in the range − M (cid:12) / pc .The model with the smallest δ , wBH-1, has N = 210 k and ρ h0 = 80 M (cid:12) / pc . It reproduces bothobservables within 3% and it has a present-day mass of . × M (cid:12) , r h = 18 . pc and f BH = 0 . .For the noBH models, we first ran five models with N = [100 k, k, k, k, k ] and ρ h0 = 10 M (cid:12) / pc , two models with N = [100 k, k ] and ρ h0 = 30 M (cid:12) / pc and two modelswith N = 100 k and ρ h0 = [100 , M (cid:12) / pc . Models with ρ h0 = 10 M (cid:12) / pc gave resultssimilar to the observations. We then ran an additional 11 models with densities ρ h0 (cid:46) M (cid:12) / pc and k ≤ N ≤ k and found that the model with the lowest δ has N = 350 k and ρ h0 =9 . M (cid:12) / pc . All model results are summarised in Tables 1 and 2.15e determine how sensitive the final parameters are to variations in the initial parameters.For a quantity A , we express the difference between its value for model i ( A i ) and the best model( A ) as ∆ log A = log( A i /A ) , where A can be an initial property ( N or ρ h0 ) or an observableproperty. For the latter we use N cluster and the number density within the half-light radius: ρ eff =3 N cluster / (8 πR ) . We can write the variation in the final properties in terms of the initial propertiesas (cid:18) ∆ log N cluster ∆ log ρ eff (cid:19) = Σ (cid:18) ∆ log N ∆ log ρ h0 (cid:19) . (4)Here Σ is a matrix that contains the constants that relate variations in initial parameter to variationsin the final parameters. We find the four elements of Σ from the two models that are nearest to theobservations (i.e. wBH-2, wBH-3 and noBH-2, noBH-3) Σ wBH = (cid:18) − . − . − . − . (cid:19) , Σ noBH = (cid:18) .
47 38 . − .
71 122 (cid:19) . (5)Absolute values of 1 mean that a fractional change in an initial parameter leads to the same frac-tional change in the final parameter. The absolute values > in the left column of Σ wBH meanthat the final properties are most sensitive to N , which is because of the collisional nature of thewBH models. The > value in the right column of Σ noBH show that the final parameters are mostsensitive to the initial density, which is because of the collisionless nature of these models. Takingthe inverse of the Σ matrices, and assuming small variations, i.e. ∆ log A (cid:39) log(1 + (cid:15) A ) ∝ (cid:15) A ,with (cid:15) A << , we can write (cid:18) (cid:15) N (cid:15) ρ h0 (cid:19) wBH (cid:39) (cid:18) . − . − .
26 2 . (cid:19) (cid:18) (cid:15) N cluster (cid:15) ρ eff (cid:19) (6)and (cid:18) (cid:15) N (cid:15) ρ h0 (cid:19) noBH (cid:39) (cid:18) . − . . . (cid:19) (cid:18) (cid:15) N cluster (cid:15) ρ eff (cid:19) . (7)This shows that the level of fine-tuning to obtain the correct N is similar in both models, albeitmore sensitive to variations in ρ eff for the noBH models. However, we find different behaviourfor the initial density, ρ h0 : for the wBH models, variations in N cluster and ρ eff allow for largervariations in ρ h0 , meaning that a relatively large range of initial densities can contribute to the errorbars on the present-day properties. The results of the noBH models are extremely sensitive to theinitial density, because variations in the observed properties correspond to variations of less thana per cent in the initial density. We can also estimate what fraction of the parameter space of theinitial conditions is covered by the uncertainties in N and ρ h0 . We assume that the initial clusterproperties are sampled from power-law distributions with indices − for N and − for ρ h0 in the ranges ≤ N ≤ × and ≤ ρ h0 / ( M (cid:12) pc − ) ≤ . Then we find that the initialconditions of wBH-1(noBH-1) cover a fraction 1/200(1/ . × ) of the area. Given that the MilkyWay has ∼ GCs, this exercise shows that finding Pal 5 is probable in the wBH scenario, whilethe probability in the noBH scenario is − . Rate of growth of the tidal tails and their visibility.
From the escaping stars we find that half ofthe observable stars in the tails of wBH-1 were ejected in the final 3 Gyr. To estimate the number16f MW field stars in Figure 4 that share the same locus in colour-magnitude as Pal 5, we use theCFHT data and colour-magnitude selection criteria from Ibata et al. . Additionally, we select onlystars more than 0.4 deg away from the best-fit stream track from Erkal et al. and adopt theirmagnitude limits of . < g CFHT < . mag. This sample is dominated by MW field stars andhas an average density of 0.142 stars/arcmin . Note that Pal 5 stream runs roughly at a constant b in the region explored in Figure 5, thus variations in the MW field density with position arenegligible. Predictions for observations.
To estimate whether the BH population can be detected from thekinematics, we look at the velocity dispersion of the wBH-1 model. Within R eff ( ∼ . arcmin)from the centre, there are 40 giant stars, with a line-of-sight velocity dispersion of . ± . km/s.Because our wBH-1 model has a steeper mass function than Pal 5, the mass of wBH-1 is likelytoo high. We therefore scale our model results to a conservative mass of × M (cid:12) (excludingBHs) and r h = 25 pc. This reduces the predicted dispersion of the giants to . ± . km/s. Inthe noBH-1 model we find 51 giants within R eff , which have a (scaled) line-of-sight dispersion of . ± . km/s. We then fit King models and a constant background using LIMEPY to thenumber density profile of bright main sequence stars shown in Figure 1(a). In Figure 6(a,b) weshow the results. For wBH-1 we find a dimensionless central concentration of W = 5 . ± . , r h = 27 . ± . pc, which results in R eff = 20 . ± . pc. For noBH-1 we find W = 5 . ± . , r h = 31 . ± . pc, which results in R eff = 23 . ± . pc. We then derive the line-of-sightvelocity dispersion of the King models by adopting a total mass (excluding the BHs) of: × M (cid:12) and r h = 25 pc. The resulting velocity dispersion profiles are shown in Figure 6(c,d). Thecentral dispersion of the King models is ∼ . km/s. The dispersion of the giants in the noBH-1 model agrees with this, which means that the derived dispersion provides an accurate measureof the total mass when assuming that the surface density traces the total mass. For the wBH-1model, however, the dispersion of the giants is about 50% (or 200 m/s) higher. A small part ofthis difference can be explained by the fact that wBH-1 is a factor of / (1 − . (cid:39) . moremassive, because of the BHs. This higher mass increases the dispersion by a factor of √ . (cid:39) . .The additional factor of 1.3 needed to explain the dispersion of the giants is because the BHsare centrally concentrated, inflating the central dispersion more than when we assume that massfollows light. Although it is challenging to find a 200 m/s velocity difference, it is feasible withexisting high-resolution spectographs ( R ∼
20 000 ) on m-class telescopes with a multi-epochobserving strategy. Baumgardt & Hilker present a compilation of line-of-sight velocities of MilkyWay GCs. There are 32 stars in their Pal 5 sample, and the dispersion of the inner 15 stars is . +0 . − . km/s. From a comparison of N -body models to the kinematics and surface brightnessprofiles the authors derive a central dispersion of 0.6 km/s. For such a low dispersion the orbitalmotions of binary stars are important. The authors have repeat observations for about half of theirstars, and they reject stars with large velocity variations. Their result supports the BH hypothesis,but a more thorough analysis of the binary content is desirable. Pal 5 has approximately stars brighter than g < mag within 3 arcmin from the center for which (additional) line-of-sight velocities can be obtained. For individual measurement errors of 300 m/s and a velocitydispersion of 400 m/s, the estimated uncertainty in the velocity dispersion for 50(100) stars is63(4) m/s, i.e. smaller than the difference between wBH-1 and noBH-1. Similar uncertainties17an be obtained even when simultaneously fitting on the velocity dispersion and the properties ofbinary stars .An additional critical prediction of the BH scenario is for the properties of the binary stars.A BH population affects binary properties. Soft binaries are ionised when they interact with stars and the binaries with the lowest binding energy are therefore indicators of the most energetic clustermembers, including invisible remnants. From wBH-1 we find that the average kinetic energy ofthe BHs is (cid:104) K BH (cid:105) = (cid:104) . mv (cid:105) (cid:39) M (cid:12) (km / s) , while for all the other stars, white dwarfs andneutron stars we find (cid:104) K ∗ (cid:105) (cid:39) . M (cid:12) (km / s) . Adopting circular orbits and equal-mass binarycomponents of . M (cid:12) (i.e. the mass of giants for which we can obtain velocities) we find thatthe orbital period is capped at P max (cid:39) . × (120) yr by the stars(BHs). However, findingno binaries with P > yr does not confirm the presence of BHs, because the cluster wasinitially more massive and compact, such that soft binaries had shorter periods. From the initialmass and half-mass radius of wBH-1(noBH-1) we find that the initial (cid:104) K (cid:105) was a factor of higher, implying that in the case there are BHs, the maximum binary period is P max (cid:39) . year,while in the noBH hypothesis the maximum period is P max (cid:39) yr. This suggests that thepresence of binaries with . (cid:46) P/ yr (cid:46) would be an argument against the presence ofa BH population, while the absence of such binaries combined with more detailed predictionsfrom N -body simulations with primordial binaries could be used as support for the BH case. For P = [1 , , yr, the orbital velocities are [12 , . , . km/s, hence these binaries would beeasily detectable with a baseline of weeks/months and moderate spectral resolution. In addition,multi-epoch kinematics serves to look for stars with a (detached) BH binary companion such asfound in NGC3201
26, 84 . Although dynamical formation of BH binaries is more common, a BHwith stellar companion can form in an exchange interaction .We note that several studies have put forward observational signatures of BH populations inclusters by using dynamical models. These studies used either the degree of mass segregation
30, 86–88 or the central surface brightness to make predictions for the size of the BH populations in a subsetof Milky Way GCs. None of these studies included Pal 5, so we can not make a comparison. Wecan check whether these methods would be able to infer the correct BH population given our modelproperties. From the scaling between f BH and the ratio of core over half-light radius ( R core /R eff )shown in figure 7 of Weatherford et al. and the ratio R core /R eff (cid:39) . in wBH-1, we would inferthat Pal 5 has a BH population of (only) f BH (cid:39) . − . The reason the inferred f BH from theirrelation is so different could be because R core /R eff saturates or goes down again for vary large f BH where they have no models, or because of a systematic difference in R core between the MonteCarlo and N -body methods
90, 91 . Microlensing event rate.
For a BH mass of M (cid:12) at a distance of 20 kpc, the Einstein angle is θ E (cid:39) . mas. In wBH-1 there are 124 BHs, mostly within R eff (cid:39) arcmin, resulting in a surfacedensity of Σ BH (cid:39) . × − mas − . The proper motion of Pal 5 is µ (cid:39) . mas/yr, such that theevent rate for a single background source (i.e. an AGN) is θ E Σ BH µ (cid:39) − /yr. In the recent AGNcatalogue from the Gaia and WISE surveys we find 2 AGNs within R eff of Pal 5, too low to get toa reasonable event rate. The event rate would be larger if we also consider background halo stars,but even with the estimated density of background galaxies in LSST ( ∼ / arcmin ), the lensingevent rate would be too low ( ∼ − / yr). 18 ata availability. All N -body data are available upon request.1. Bernard, E. J. et al. A Synoptic Map of Halo Substructures from the Pan-STARRS1 3 π Survey.
Mon. Not. R. Astron. Soc. , 1759–1768 (2016). .2. Shipp, N. et al.
Stellar Streams Discovered in the Dark Energy Survey.
Astrophysical Journal , 114 (2018). .3. Malhan, K., Ibata, R. A. & Martin, N. F. Ghostly tributaries to the Milky Way: charting thehalo’s stellar streams with the Gaia DR2 catalogue.
Mon. Not. R. Astron. Soc. , 3442–3455(2018). .4. Ibata, R. A., Malhan, K. & Martin, N. F. The Streams of the Gaping Abyss: A Populationof Entangled Stellar Streams Surrounding the Inner Galaxy.
Astrophysical Journal , 152(2019). .5. Koposov, S. E., Rix, H.-W. & Hogg, D. W. Constraining the Milky Way Potential with aSix-Dimensional Phase-Space Map of the GD-1 Stellar Stream.
Astrophysical Journal ,260–273 (2010). .6. de Boer, T. J. L., Erkal, D. & Gieles, M. A closer look at the spur, blob, wiggle, and gaps inGD-1.
Mon. Not. R. Astron. Soc. , 5315–5332 (2020). .7. Baumgardt, H. & Makino, J. Dynamical evolution of star clusters in tidal fields.
Mon. Not. R.Astron. Soc. , 227–246 (2003).8. Kuzma, P. B., Da Costa, G. S. & Mackey, A. D. The outer envelopes of globular clusters. II.NGC 1851, NGC 5824 and NGC 1261 ∗ . Mon. Not. R. Astron. Soc. , 2881–2898 (2018). .9. Myeong, G. C., Evans, N. W., Belokurov, V., Sand ers, J. L. & Koposov, S. E. The SausageGlobular Clusters.
Astrophysical Journal Letters , L28 (2018). .10. Odenkirchen, M. et al.
Detection of Massive Tidal Tails around the Globular Cluster Palomar5 with Sloan Digital Sky Survey Commissioning Data.
Astrophysical Journal Letters ,L165–L169 (2001). astro-ph/0012311 .11. Ibata, R. A., Lewis, G. F., Thomas, G., Martin, N. F. & Chapman, S. Feeling the Pull: A Studyof Natural Galactic Accelerometers. II. Kinematics and Mass of the Delicate Stellar Stream ofthe Palomar 5 Globular Cluster.
Astrophysical Journal , 120 (2017). .12. Dehnen, W., Odenkirchen, M., Grebel, E. K. & Rix, H.-W. Modeling the Disruption of theGlobular Cluster Palomar 5 by Galactic Tides.
Astronomical Journal , 2753–2770 (2004).13. Gieles, M., Heggie, D. C. & Zhao, H. The life cycle of star clusters in a tidal field.
Mon. Not.R. Astron. Soc. , 2509–2524 (2011). .194. Smith, G. H., Sneden, C. & Kraft, R. P. A Study of Abundances of Four Giants in the Low-Mass Globular Cluster Palomar 5.
Astronomical Journal , 1502–1508 (2002).15. Elmegreen, B. G. Globular Cluster Formation at High Density: A Model for Elemental En-richment with Fast Recycling of Massive-star Debris.
Astrophysical Journal , 80 (2017). .16. Gieles, M. et al.
Concurrent formation of supermassive stars and globular clusters: im-plications for early self-enrichment.
Mon. Not. R. Astron. Soc. , 2461–2479 (2018). .17. Abbott, B. P. et al.
Observation of Gravitational Waves from a Binary Black Hole Merger.
Physical Review Letters , 061102 (2016). .18. Fryer, C. L. et al.
Compact Remnant Mass Function: Dependence on the Explosion Mecha-nism and Metallicity.
Astrophysical Journal , 91 (2012). .19. Hobbs, G., Lorimer, D. R., Lyne, A. G. & Kramer, M. A statistical study of 233 pulsar propermotions.
Mon. Not. R. Astron. Soc. , 974–992 (2005). astro-ph/0504584 .20. Merritt, D., Piatek, S., Portegies Zwart, S. & Hemsendorf, M. Core Formation by a Populationof Massive Remnants.
Astrophysical Journal Letters , L25–L28 (2004). astro-ph/0403331 .21. Mackey, A. D., Wilkinson, M. I., Davies, M. B. & Gilmore, G. F. Black holes and coreexpansion in massive star clusters.
Mon. Not. R. Astron. Soc. , 65–95 (2008). .22. Giersz, M. et al.
MOCCA survey data base- I. Dissolution of tidally filling star clus-ters harbouring black hole subsystems.
Mon. Not. R. Astron. Soc. , 2412–2423 (2019). .23. Wang, L. The survival of star clusters with black hole subsystems.
Mon. Not. R. Astron. Soc. , 2413–2423 (2020). .24. Strader, J., Chomiuk, L., Maccarone, T. J., Miller-Jones, J. C. A. & Seth, A. C. Two stellar-mass black holes in the globular cluster M22.
Nature , 71–73 (2012). .25. Chomiuk, L. et al.
A Radio-selected Black Hole X-Ray Binary Candidate in the Milky WayGlobular Cluster M62.
Astrophysical Journal , 69 (2013). .26. Giesers, B. et al.
A detached stellar-mass black hole candidate in the globular cluster NGC3201.
Mon. Not. R. Astron. Soc. , L15–L19 (2018).27. Wang, L. et al.
NBODY6++GPU: ready for the gravitational million-body problem.
Mon.Not. R. Astron. Soc. , 4070–4080 (2015). .208. Banerjee, S. et al.
BSE versus StarTrack: Implementations of new wind, remnant-formation,and natal-kick schemes in NBODY7 and their astrophysical consequences.
Astron. & Astro-phys. , A41 (2020). .29. Hurley, J. R. Ratios of star cluster core and half-mass radii: a cautionary note on intermediate-mass black holes in star clusters.
Mon. Not. R. Astron. Soc. , 93–99 (2007). .30. Peuten, M., Zocchi, A., Gieles, M., Gualandris, A. & H´enault-Brunet, V. A stellar-mass blackhole population in the globular cluster NGC 6101?
Mon. Not. R. Astron. Soc. , 2333–2342(2016). .31. Breen, P. G. & Heggie, D. C. Dynamical evolution of black hole subsystems in idealized starclusters.
Mon. Not. R. Astron. Soc. , 2779–2797 (2013). .32. Erkal, D., Koposov, S. E. & Belokurov, V. A sharper view of Pal 5’s tails: discovery of streamperturbations with a novel non-parametric technique.
Mon. Not. R. Astron. Soc. , 60–84(2017). .33. Banik, N. & Bovy, J. Effects of baryonic and dark matter substructure on the Pal 5 stream.
Mon. Not. R. Astron. Soc. , 2009–2020 (2019). .34. Kuzma, P. B., Da Costa, G. S., Keller, S. C. & Maunder, E. Palomar 5 and its tidal tails: asearch for new members in the tidal stream.
Mon. Not. R. Astron. Soc. , 3297–3309 (2015). .35. Banerjee, S. & Kroupa, P. A New Type of Compact Stellar Population: Dark Star Clusters.
Astrophysical Journal Letters , L12 (2011). .36. Baumgardt, H., Parmentier, G., Gieles, M. & Vesperini, E. Evidence for two populations ofGalactic globular clusters from the ratio of their half-mass to Jacobi radii.
Mon. Not. R. Astron.Soc. , 1832–1838 (2010). .37. Baumgardt, H. & Hilker, M. A catalogue of masses, structural parameters, and velocity disper-sion profiles of 112 Milky Way globular clusters.
Mon. Not. R. Astron. Soc. , 1520–1557(2018). .38. Elmegreen, B. G. The Globular Cluster Mass Function as a Remnant of Violent Birth.
Astro-physical Journal Letters , L184–L188 (2010). .39. Kruijssen, J. M. D. Globular clusters as the relics of regular star formation in ‘normal’ high-redshift galaxies.
Mon. Not. R. Astron. Soc. , 1658–1686 (2015). .40. Spitzer, J., Lyman. Disruption of Galactic Clusters.
Astrophysical Journal , 17 (1958).41. Gieles, M. et al.
Star cluster disruption by giant molecular clouds.
Mon. Not. R. Astron. Soc. , 793–804 (2006). astro-ph/0606451 .212. Gieles, M. & Renaud, F. If it does not kill them, it makes them stronger: collisional evolutionof star clusters with tidal shocks.
Mon. Not. R. Astron. Soc. , L103–L107 (2016). .43. Massari, D., Koppelman, H. H. & Helmi, A. Origin of the system of globular clusters in theMilky Way.
Astron. & Astrophys. , L4 (2019). .44. Bianchini, P., Renaud, F., Gieles, M. & Varri, A. L. The inefficiency of satellite accretion informing extended star clusters.
Mon. Not. R. Astron. Soc. , L40–L44 (2015). .45. Chatterjee, S., Rodriguez, C. L. & Rasio, F. A. Binary Black Holes in Dense Star Clusters: Ex-ploring the Theoretical Uncertainties.
Astrophysical Journal , 68 (2017). .46. Kim, J.-h. et al.
Formation of globular cluster candidates in merging proto-galaxies at highredshift: a view from the FIRE cosmological simulations.
Mon. Not. R. Astron. Soc. ,4232–4244 (2018). .47. Rodriguez, C. L., Chatterjee, S. & Rasio, F. A. Binary black hole mergers from globularclusters: Masses, merger rates, and the impact of stellar evolution.
Phys. Review D , 084029(2016). .48. Kremer, K. et al. Modeling Dense Star Clusters in the Milky Way and Beyond with the CMCCluster Catalog.
Astrophysical Journal Supplement , 48 (2020). .49. Antonini, F. & Gieles, M. Merger rate of black hole binaries from globular clusters: Theoreti-cal error bars and comparison to gravitational wave data from GWTC-2.
Phys. Review D ,123016 (2020). .50. Krumholz, M. R., McKee, C. F. & Bland -Hawthorn, J. Star Clusters Across Cosmic Time.
Annu. Rev. Astron. Astrophys , 227–303 (2019). .51. Vesperini, E. Evolution of globular cluster systems in elliptical galaxies - II. Power-law initialmass function. Mon. Not. R. Astron. Soc. , 247–256 (2001). astro-ph/0010111 .52. Fall, S. M. & Zhang, Q. Dynamical Evolution of the Mass Function of Globular Star Clusters.
Astrophysical Journal , 751–765 (2001).53. Sollima, A., Mart´ınez-Delgado, D., Valls-Gabaud, D. & Pe ˜narrubia, J. Discovery of TidalTails Around the Distant Globular Cluster Palomar 14.
Astrophysical Journal , 47 (2011). .54. Bovy, J. galpy: A python Library for Galactic Dynamics.
Astrophysical Journal Supplement , 29 (2015). .55. Price-Whelan, A. M. et al.
Kinematics of the Palomar 5 Stellar Stream from RR Lyrae Stars.
Astronomical Journal , 223 (2019). .226. Vasiliev, E. Proper motions and dynamics of the Milky Way globular cluster system from GaiaDR2.
Mon. Not. R. Astron. Soc. , 2832–2850 (2019). .57. Gravity Collaboration et al.
A geometric distance measurement to the Galactic center blackhole with 0.3% uncertainty.
Astron. & Astrophys. , L10 (2019). .58. Navarro, J. F., Frenk, C. S. & White, S. D. M. The Structure of Cold Dark Matter Halos.
Astrophysical Journal , 563 (1996). astro-ph/9508025 .59. Miyamoto, M. & Nagai, R. Three-dimensional models for the distribution of mass in galaxies.
PASJ , 533–543 (1975).60. Hernquist, L. An analytical model for spherical galaxies and bulges. Astrophysical Journal , 359–364 (1990).61. Martell, S. L., Smith, G. H. & Grillmair, C. J. A New Age Measurement for Palomar 5. In
American Astronomical Society Meeting Abstracts , vol. 201 of
American Astronomical SocietyMeeting Abstracts , 07.11 (2002).62. Dotter, A., Sarajedini, A. & Anderson, J. Globular Clusters in the Outer Galactic Halo: NewHubble Space Telescope/Advanced Camera for Surveys Imaging of Six Globular Clusters andthe Galactic Globular Cluster Age-metallicity Relation.
Astrophysical Journal , 74 (2011). .63. Xu, X. et al.
New Determination of Fundamental Properties of Palomar 5 Using Deep DESIImaging Data.
Astronomical Journal , 12 (2021). .64. Plummer, H. C. On the problem of distribution in globular star clusters.
Mon. Not. R. Astron.Soc. , 460–470 (1911).65. Kroupa, P. On the variation of the initial mass function. Mon. Not. R. Astron. Soc. ,231–246 (2001).66. Aarseth, S. J. From NBODY1 to NBODY6: The Growth of an Industry.
PASP , 1333–1346(1999).67. Aarseth, S. J.
Gravitational N-Body Simulations (Cambridge University Press, November2003., 2003).68. Ahmad, A. & Cohen, L. A numerical integration scheme for the N-body gravitational problem.
Journal of Computational Physics , 389–402 (1973).69. Makino, J. & Aarseth, S. J. On a Hermite integrator with Ahmad-Cohen scheme for gravita-tional many-body problems. PASJ , 141–151 (1992).70. Hurley, J. R., Pols, O. R. & Tout, C. A. Comprehensive analytic formulae for stellar evolutionas a function of mass and metallicity. Mon. Not. R. Astron. Soc. , 543–569 (2000). arXiv:astro-ph/0001295 . 231. Hurley, J. R., Tout, C. A. & Pols, O. R. Evolution of binary stars and the effect of tides onbinary populations.
Mon. Not. R. Astron. Soc. , 897–928 (2002). arXiv:astro-ph/0201220 .72. Nitadori, K. & Aarseth, S. J. Accelerating NBODY6 with graphics processing units.
Mon.Not. R. Astron. Soc. , 545–552 (2012). .73. Belczynski, K. et al.
Compact Object Modeling with the StarTrack Population Synthesis Code.
Astrophysical Journal Supplement , 223–260 (2008). astro-ph/0511811 .74. King, I. R. The structure of star clusters. III. Some simple dynamical models.
AstronomicalJournal , 64 (1966).75. Ibata, R. et al. Do globular clusters possess dark matter haloes? A case study in NGC 2419.
Mon. Not. R. Astron. Soc. , 3648–3659 (2013). .76. Gieles, M. & Zocchi, A. A family of lowered isothermal models.
Mon. Not. R. Astron. Soc. , 576–592 (2015). .77. Foreman-Mackey, D., Hogg, D. W., Lang, D. & Goodman, J. emcee: The MCMC Hammer.
PASP , 306–312 (2013). .78. Choi, J. et al.
Mesa Isochrones and Stellar Tracks (MIST). I. Solar-scaled Models.
Astrophys-ical Journal , 102 (2016). .79. Dotter, A. MESA Isochrones and Stellar Tracks (MIST) 0: Methods for the Construction ofStellar Isochrones.
Astrophysical Journal Supplement , 8 (2016). .80. Ibata, R. A., Lewis, G. F. & Martin, N. F. Feeling the Pull: a Study of Natural Galactic Ac-celerometers. I. Photometry of the Delicate Stellar Stream of the Palomar 5 Globular Cluster.
Astrophysical Journal , 1 (2016). .81. Kravtsov, A. V. & Gnedin, O. Y. Formation of Globular Clusters in Hierarchical Cosmology.
Astrophysical Journal , 650–665 (2005). arXiv:astro-ph/0305199 .82. Cottaar, M., Meyer, M. R. & Parker, R. J. Characterizing a cluster’s dynamic state using asingle epoch of radial velocities.
Astron. & Astrophys. , A35 (2012). .83. Heggie, D. C. Binary evolution in stellar dynamics.
Mon. Not. R. Astron. Soc. , 729–787(1975).84. Giesers, B. et al.
A stellar census in globular clusters with MUSE: Binaries in NGC 3201.
Astron. & Astrophys. , A3 (2019). .85. Kremer, K., Ye, C. S., Chatterjee, S., Rodriguez, C. L. & Rasio, F. A. How Black Holes ShapeGlobular Clusters: Modeling NGC 3201.
Astrophysical Journal Letters , L15 (2018). . 246. Alessandrini, E., Lanzoni, B., Ferraro, F. R., Miocchi, P. & Vesperini, E. Investigating theMass Segregation Process in Globular Clusters with Blue Straggler Stars: The Impact of DarkRemnants.
Astrophysical Journal , 252 (2016). .87. Weatherford, N. C., Chatterjee, S., Rodriguez, C. L. & Rasio, F. A. Predicting Stellar-massBlack Hole Populations in Globular Clusters.
Astrophysical Journal , 13 (2018). .88. Weatherford, N. C., Chatterjee, S., Kremer, K. & Rasio, F. A. A Dynamical Survey of Stellar-mass Black Holes in 50 Milky Way Globular Clusters.
Astrophysical Journal , 162 (2020). .89. Askar, A., Arca Sedda, M. & Giersz, M. MOCCA-SURVEY Database I: Galactic globularclusters harbouring a black hole subsystem.
Mon. Not. R. Astron. Soc. , 1844–1854 (2018). .90. Rodriguez, C. L. et al.
Million-body star cluster simulations: comparisons between MonteCarlo and direct N-body.
Mon. Not. R. Astron. Soc. , 2109–2118 (2016). .91. Rodriguez, C. L. et al.
A new hybrid technique for modeling dense star clusters.
Computa-tional Astrophysics and Cosmology , 5 (2018). .92. Shu, Y. et al. Catalogues of active galactic nuclei from Gaia and unWISE data.
Mon. Not. R.Astron. Soc. , 4741–4759 (2019). .93. Astropy Collaboration et al.
Astropy: A community Python package for astronomy.
Astron.& Astrophys. , A33 (2013). .94. Astropy Collaboration et al.
The Astropy Project: Building an Open-science Project andStatus of the v2.0 Core Package.
Astronomical Journal , 123 (2018). .95. H´enon, M. Sur l’´evolution dynamique des amas globulaires.
Ann. Astrophys. , 369 (1961).25 cknowledgements MG and EB acknowledge financial support from the European Research Council(ERC StG-335936, CLUSTERS) and MG acknowledges support from the Ministry of Science and Inno-vation through a Europa Excelencia grant (EUR2020-112157). EB acknowledges financial support from aVici grant from the Netherlands Organisation for Scientific Research (NWO). MG thanks Mr Gaby P´erezForcadell for installing the GPU server at the ICCUB on which all the simulations were run. The authorsthank Rodrigo Ibata for sharing the data of Pal 5’s surface density profile, Łukasz Wyrzykowski for discus-sions on microlensing and Sverre Aarseth, Keigo Nitadori and Long Wang for maintaining
NBODY
NBODY
GPU and making the codes publicly available. MG and FA thank Long Wang and SambaranBanerjee for discussions on the recent
SSE and
BSE updates and the implementation in
NBODY
GPU .This research made use of
ASTROPY , a community-developed core Python package for Astronomy
93, 94
Competing Interests
The authors have no competing financial interests.
Author contributions
MG ran all N -body simulations, analysed them and was in charge of the writing.DE was in charge of stream modelling and deriving the orbit of Pal 5 and the parameters of the MW model.FA contributed to the BH physics of the N -body models and paper. EB converted stream models to observedquantities and JP contributed to the binary properties. All authors assisted in the development, analysis andwriting of the paper. Correspondence
1) (2) (3) (4) (5) (6) (7) (8) (9) (10) (11) (12) (13) (14) N ρ h0 M r h0 T N cluster R eff M r h f BH δ N δ R δ Model[ ] [ M (cid:12) / pc ] [ M (cid:12) ] [pc] [Gyr] [pc] [ M (cid:12) ] [pc] [%] [%] [%] [%]50 30 0.319 5.04 2.74 – – – – – – – – wBH-1750 100 0.319 3.38 4.06 – – – – – – – – wBH-1650 300 0.319 2.34 5.58 – – – – – – – – wBH-1450 1000 0.319 1.57 8.34 1052 4.05 0.483 5.29 0.701 -32.1 -78.3 84.7 wBH-9100 30 0.638 6.35 5.40 – – – – – – – – wBH-15100 100 0.638 4.25 7.85 – – – – – – – – wBH-13100 300 0.638 2.95 11.36 1448 9.8 0.732 12.8 3.09 -6.58 -47.6 48 wBH-5100 1000 0.638 1.97 17.99 3541 6.26 1.84 9.04 0.321 +128 -66.5 145 wBH-10200 30 1.28 7.98 8.12 – – – – – – – – wBH-12200 100 1.28 5.34 11.43 1399 17.7 0.86 20 20.1 -9.74 -5.3 11.1 wBH-3200 300 1.28 3.7 17.64 7210 11.1 4.05 14.9 2.75 +365 -40.8 367 wBH-11210 80 1.34 5.85 11.47 1497 18.1 0.949 18.8 22.3 -3.42 -3.11 4.62 wBH-1215 80 1.37 5.9 11.13 459 19.9 0.419 15.3 44.6 -70.4 +6.15 70.7 wBH-8220 70 1.4 6.21 11.45 1411 18.5 0.953 19.8 21.9 -8.92 -0.948 8.97 wBH-2220 80 1.4 5.94 11.31 791 20.8 0.604 19.4 33 -49 +11.5 50.3 wBH-6225 70 1.44 6.26 11.29 1091 16.6 0.77 16.5 29.1 -29.6 -11.4 31.7 wBH-4225 80 1.44 5.98 11.71 2372 16.8 1.45 19.8 15.1 +53.1 -10 54 wBH-7Table 1: Overview of the 17 wBH N -body simulations of Pal 5.The different columns present: (1) initial number of stars; (2) ini-tial half-mass density; (3) initial mass; (4) initial half-mass ra-dius; (5) age when the cluster had the observed number of stars( N cluster = 1550 ); (6) observable number of stars at 11.5 Gyr; (7)half-light radii at 11.5 Gyr; (8) bound mass at 11.5 Gyr; (9) half-mass radius of the bound stars at 11.5 Gyr; (10) BH mass fractionat an age of 11.5 Gyr; (11) fractional difference with N cluster ; (12)fractional difference with observed R eff = 18 . ± . pc; (13) over-all fractional difference δ = (cid:113) δ N + δ R ; (14) Model ID, sorted inlevel of agreement with the data with wBH-1 being the closestmatch.
1) (2) (3) (4) (5) (6) (7) (8) (9) (10) (11) (12) (13) (14) N ρ h0 M r h0 T N cluster R eff M r h f BH δ N δ R δ ID[ ] [ M (cid:12) / pc ] [ M (cid:12) ] [pc] [Gyr] [pc] [ M (cid:12) ] [pc] [%] [%] [%] [%]100 10 0.638 9.15 7.84 609 4.42 0.289 5.76 0 -60.7 -76.4 97.5 noBH-14100 30 0.638 6.35 15.09 2682 4.71 1.49 7.24 0 +73.1 -74.8 105 noBH-15100 100 0.638 4.25 16.71 3607 3.57 2.05 6.33 0 +133 -80.9 155 noBH-18100 300 0.638 2.95 18.62 4219 4.05 2.16 6.65 0 +172 -78.3 189 noBH-19200 9.5 1.28 11.7 10.79 1062 8.76 0.501 11.7 0 -31.5 -53.2 61.8 noBH-10200 10 1.28 11.5 11.65 1675 9.03 0.812 12 0 +8.11 -51.7 52.4 noBH-6200 30 1.28 7.98 2.35 10320 9.52 6.49 12.7 0 +566 -49.1 568 noBH-20300 9.5 1.91 13.4 11.11 1065 13.5 0.508 17.5 0.321 -31.3 -27.7 41.8 noBH-5300 9.625 1.91 13.3 11.38 1349 12.2 0.67 16 0.244 -13 -34.7 37.1 noBH-4300 9.75 1.91 13.3 11.45 1500 12.2 0.736 16.1 0.222 -3.23 -35 35.1 noBH-3300 10 1.91 13.2 12.07 2316 12.1 1.2 16.1 0.136 +49.5 -35.5 60.9 noBH-9350 9.625 2.23 14 11.37 1136 19.6 0.583 25.3 0 -26.7 +4.84 27.1 noBH-1350 9.75 2.23 14 11.63 1864 13.7 0.966 17.9 0 +20.3 -26.7 33.6 noBH-2400 9.5 2.55 14.7 11.14 81 22.5 0.0436 33.5 0 -94.8 +20.2 96.9 noBH-12400 9.625 2.55 14.7 11.40 1105 27.1 0.539 34.7 0 -28.7 +45.1 53.4 noBH-7400 9.75 2.55 14.6 11.79 2516 14.3 1.28 19.1 0 +62.3 -23.7 66.7 noBH-11400 10 2.55 14.5 12.16 3237 14.1 1.7 19.1 0 +109 -24.4 112 noBH-16500 9.5 3.19 15.9 11.14 112 24 0.0573 30.8 0 -92.8 +28.2 97 noBH-13500 9.75 3.19 15.7 11.66 2429 18.7 1.24 24.6 0 +56.7 +0.0658 56.7 noBH-8500 10 3.19 15.6 11.92 3650 15.8 1.86 21.1 0 +136 -15.3 136 noBH-17Table 2: As in Table 1 but now for the 20 noBH models. upplementary material ] 225.0227.5230.0232.5235.0237.5240.0242.5 RA [deg]
Figure 1: Comparison between stream properties of wBH-1 and noBH-1 and the observed streamfrom Erkal et al. . The stream width of both N -body models is similar over the range includedin the observations and shows some systematic deviations from the observed width. The densityprofile of the N -body models is smoother than the observed profile, which shows signatures ofover/under-densities. The decline in the density of the trailing arm (large RA) is faster in wBH-1than in noBH-1, which agrees more with the observed decline.29 ompact GCs ‘fluffy’ GCsname distance [kpc] tail data name distance [kpc] tail dataNGC 3201 4.9 long tidal tails
2, 3
NGC 288 8.9 long tidal tails
2, 4, 5
NGC 6205 7.1 tidal tails Pal 1 11.1 tidal tails NGC 6341 8.3 long tidal tails NGC 6101 15.4 -NGC 362 8.6 tidal feature NGC 5466 16.0 tidal tails
9, 10
NGC 6779 9.4 - NGC 5053 17.4 tidal feature NGC 2808 9.6 - IC 4499 18.8 -NGC 5272 10.2 - BH 176 18.9 -NGC 4590 10.3 long tidal tails Pal 12 19.0 tidal tails NGC 7078 10.4 - NGC 6426 20.6 -NGC 2298 10.8 tidal features Rup 106 21.2 -NGC 7089 11.5 - ESO 280 21.4 -NGC 5286 11.7 - Ter 7 22.8 -NGC 1851 12.1 long tidal tails
2, 5
Pal 5 23.2 long tidal tails
2, 13, 14
NGC 1904 12.9 tidal tails IC 1257 25.0 -NGC 6934 15.6 - Pal 13 26.0 long tidal tails NGC 1261 16.3 long tidal tails
2, 5
Ter 8 26.3 -NGC 5024 17.9 - NGC 7492 26.3 tidal tails NGC 6981 17.0 - Arp 2 28.6 -NGC 4147 19.4 tidal feature AM 4 32.2 -NGC 6864 20.9 - Pyxis 39.4 -NGC 5634 25.2 - Pal 15 45.1 tidal tails Pal 2 27.2 - Pal 14 76.5 tidal tails NGC 6229 30.5 - Eridanus 90.1 tidal tails NGC 5824 32.1 - Pal 3 92.5 tidal feature NGC 5694 35.0 - Pal 4 108.7 tidal feature NGC 7006 41.2 - AM 1 123.3 -Table 1: Summary of results of stream searches for GCs at dis-tances > kpc from the Galactic centre. The classification ‘com-pact’ (left) and ‘fluffy’ (right) is from Baumgardt et al. . All clus-ters are sorted in distance from the Sun. Among the compact clus-ters, no tidal tails nor tidal features were found for clusters that aremore than 20 kpc away, while they were found for about half ofthe fluffy clusters beyond 20 kpc.
1. Erkal, D., Koposov, S. E. & Belokurov, V. A sharper view of Pal 5’s tails: discovery of streamperturbations with a novel non-parametric technique.
Mon. Not. R. Astron. Soc. , 60–84(2017). .2. Ibata, R. et al.
Charting the Galactic acceleration field I. A search for stellar streams with GaiaDR2 and EDR3 with follow-up from ESPaDOnS and UVES. arXiv e-prints arXiv:2012.05245(2020). . 30. Palau, C. G. & Miralda-Escud´e, J. The tidal stream generated by the globular cluster NGC3201. arXiv e-prints arXiv:2010.14381 (2020). .4. Kaderali, S., Hunt, J. A. S., Webb, J. J., Price-Jones, N. & Carlberg, R. Rediscovering thetidal tails of NGC 288 with Gaia DR2.
Mon. Not. R. Astron. Soc. , L114–L118 (2019). .5. Shipp, N. et al.
Stellar Streams Discovered in the Dark Energy Survey.
Astrophysical Journal , 114 (2018). .6. Leon, S., Meylan, G. & Combes, F. Tidal tails around 20 Galactic globular clusters. Observa-tional evidence for gravitational disk/bulge shocking.
Astron. & Astrophys. , 907–931 (2000). astro-ph/0006100 .7. Niederste-Ostholt, M. et al.
The tidal tails of the ultrafaint globular cluster Palomar 1.
Mon.Not. R. Astron. Soc. , L66–L70 (2010). .8. Carballo-Bello, J. A. Using Gaia DR2 to detect extratidal structures around the Galactic glob-ular cluster NGC 362.
Mon. Not. R. Astron. Soc. , 1667–1671 (2019). .9. Belokurov, V., Evans, N. W., Irwin, M. J., Hewett, P. C. & Wilkinson, M. I. The Discovery ofTidal Tails around the Globular Cluster NGC 5466.
Astrophysical Journal Letters , L29–L32(2006). astro-ph/0511767 .10. Bernard, E. J. et al.
A Synoptic Map of Halo Substructures from the Pan-STARRS1 3 π Survey.
Mon. Not. R. Astron. Soc. , 1759–1768 (2016). .11. Jordi, K. & Grebel, E. K. Search for extratidal features around 17 globular clusters in theSloan Digital Sky Survey.
Astron. & Astrophys. , A71 (2010). .12. Balbinot, E., Santiago, B. X., da Costa, L. N., Makler, M. & Maia, M. A. G. The tidal tails ofNGC 2298.
Mon. Not. R. Astron. Soc. , 393–402 (2011). .13. Odenkirchen, M. et al.
Detection of Massive Tidal Tails around the Globular Cluster Palomar5 with Sloan Digital Sky Survey Commissioning Data.
Astrophysical Journal Letters , L165–L169 (2001). astro-ph/0012311 .14. Bonaca, A. et al.
Variations in the Width, Density, and Direction of the Palomar 5 Tidal Tails.
Astrophysical Journal , 70 (2020). .15. Shipp, N., Price-Whelan, A. M., Tavangar, K., Mateu, C. & Drlica-Wagner, A. Discovery ofExtended Tidal Tails around the Globular Cluster Palomar 13.
Astronomical Journal , 244(2020). .16. Navarrete, C., Belokurov, V. & Koposov, S. E. The Discovery of Tidal Tails around theGlobular Cluster NGC 7492 with Pan-STARRS1.
Astrophysical Journal Letters , L23 (2017). . 317. Myeong, G. C., Jerjen, H., Mackey, D. & Da Costa, G. S. Tidal Tails around the Outer HaloGlobular Clusters Eridanus and Palomar 15.
Astrophysical Journal Letters , L25 (2017). .18. Sollima, A., Mart´ınez-Delgado, D., Valls-Gabaud, D. & Pe ˜narrubia, J. Discovery of TidalTails Around the Distant Globular Cluster Palomar 14.
Astrophysical Journal , 47 (2011). .19. Sohn, Y.-J. et al.
Wide-Field Stellar Distributions around the Remote Young Galactic GlobularClusters Palomar 3 and Palomar 4.
Astronomical Journal , 803–814 (2003).20. Baumgardt, H., Parmentier, G., Gieles, M. & Vesperini, E. Evidence for two populations ofGalactic globular clusters from the ratio of their half-mass to Jacobi radii.
Mon. Not. R. Astron.Soc. , 1832–1838 (2010).0909.5696