The influence of black holes on the binary population of the globular cluster Palomar 5

The discovery of stellar-mass black holes (BHs) in globular clusters (GCs) raises the possibility of long-term retention of BHs within GCs. These BHs influence various astrophysical processes, including merger-driven gravitational waves and the formation of X-ray binaries. They also impact cluster dynamics by heating and creating low-density cores. Previous N-body models suggested that Palomar 5, a low-density GC with long tidal tails, may contain more than 100 BHs. To test this scenario, we conduct N-body simulations of Palomar 5 with primordial binaries to explore the influence of BHs on binary populations and the stellar mass function. Our results show that primordial binaries have minimal effect on the long-term evolution. In dense clusters with BHs, the fraction of wide binaries with periods>$10^5$ days decreases, and the disruption rate is independent of the initial period distribution. Multi-epoch spectroscopic observations of line-of-sight velocity changes can detect most bright binaries with periods below $10^4$ days, significantly improving velocity dispersion measurements. Four BH-MS binaries in the model with BHs suggests their possible detection through the same observation method. Including primordial binaries leads to a flatter inferred mass function because of spatially unresolved binaries, leading to a better match of the observations than models without binaries, particularly in Palomar 5's inner region. Future observations should focus on the cluster velocity dispersion and binaries with periods of $10^4-10^5$ days in Palomar 5's inner and tail regions to constrain BH existence.

Palomar 5 (Pal 5) is among the Galactic GCs renowned for its long tidal streams and unusually low central density (e.g.Rockosi et al. 2002;Odenkirchen et al. 2001Odenkirchen et al. , 2002Odenkirchen et al. , 2003;;Koch et al. 2004;Odenkirchen et al. 2009;Carlberg et al. 2012;Kuzma et al. 2015;Ishigaki et al. 2016;Price-Whelan et al. 2019;Bonaca et al. 2020;Starkman et al. 2020), which suggests the possible presence of a substantial number of BHs in the cluster (Gieles et al. 2021, hereafter G21).Understanding the properties of the BH population in Pal 5 is also crucial for explaining the pronounced nature of its stream.G21 employed self-consistent -body models that resolve individual stars to propose the existence of a large population of BHs in the cluster core (20% of the total mass), enhancing tidal disruption.However, the BH hypothesis needs further confirmation, because the observed density profiles of the cluster and the stream could also be reproduced by an -body model of a BH-free cluster with a low initial density.
The binary population of Pal 5 plays a crucial role in resolving this degeneracy.According to the Heggie (1975)-Hills (1975) law, close encounters with binaries can result in two opposing evolutionary trends: wide/soft binaries become less bound and decay with a few close encounters, while tight/hard binaries become tighter due to the increased kinetic energy of the intruder and the centre-of-mass of the binary.The boundary between these two types depends on the local kinetic energy of particles where the binary resides.G21 argue that the kinetic energy of BHs is higher than that of stars in a cluster without BHs with similar half-light radius.It is therefore expected that fewer soft binaries could survive in the case the cluster contains BHs, which is a prediction that can be tested with observations.Furthermore, due to the large distance of Pal 5, most binaries cannot be resolved spatially by current state-of-art observational instruments.Because unresolved binaries might influence the determination of velocity dispersion and present-day mass functions, it is worthwhile to investigate how primordial binaries and BHs collectively affect the line-of-sight velocity measurement and mass function and whether it can be used to indirectly constrain the existence of BHs.
In this study, we perform -body simulations of several Pal 5-like clusters with and without BHs and incorporating a large number of binaries, to examine the impact of BHs on binary disruption and the long-term evolution of Pal 5 and its tidal tails.Section 2 describes the -body simulation method, data analysis tools, and the observational data of Pal 5 utilized in this study.Section 3 presents the results of our -body models, comparing the structural evolution, surface number density, binary properties, and present-day mass function with models from G21 and observational data.Section 4 discusses the limitations of our models and outlines prospects for future observations.Finally, Section 5 concludes this work.

N-body code
We conducted simulations of Pal 5-like clusters using the highperformance -body code petar (Wang et al. 2020b).To achieve high parallel performance, the framework for developing parallel particle simulation codes (fdps) is implemented in petar (Iwasawa et al. 2016(Iwasawa et al. , 2020)).The code incorporates the particle-tree and particleparticle method (P 3 T) (Oshino et al. 2011), which enables the separate integration of long-range and short-range interactions between particles.For accurate integration of the weak long-range interactions, the code uses a Barnes & Hut (1986) particle-tree method with a 2nd-order Leap-frog integrator, which has a computational cost of  ( log()).To accurately follow orbital motions of binaries, hyperbolic encounters, and the evolution of hierarchical few-body systems, the 4th-order Hermite method along with the slowdownalgorithmic regularization (SDAR) method is used (Wang et al. 2020a).One of the major advantages of the petar code is its capability to include a large fraction of binaries, up to 100%, in the simulation of stellar systems without significant performance loss.This feature enables us to carry out the models presented in this work.
In our simulations, we included binaries with a wide period distribution (see Section 2.5), requiring the use of Leap-frog, Hermite, and SDAR integrators for integrating binary orbits.While Leap-frog and SDAR are symplectic methods that conserve energy and angular momentum, the Hermite integrator does not.We employ sufficiently small time steps for the Hermite integrator to ensure that the artificial drift of semi-major axes and eccentricities remains insignificant throughout the entire evolutionary time of all our models.The key parameters for switching the integrator and controlling the accuracy of one simulation in this work are provided below: • Changeover inner radius: 0.0027 pc • Changeover outer radius: 0.027 pc • SDAR separation criterion: 0.000216 pc • Tree time step: 0.0009765625 Myr • Hermite time step coefficient : 0.1 See Wang et al. (2020b) for the details on the definition of these parameters.
The population synthesis code for single and binary stellar evolution, sse and bse, are implemented in petar (Hurley et al. 2000(Hurley et al. , 2002)).Furthermore, the code utilizes an updated version from Banerjee et al. ( 2020) that incorporates semi-empirical stellar wind prescriptions from Belczynski et al. (2010); Vink et al. (2011), a "rapid" supernova model for remnant formation and material fallback from Fryer et al. (2012), and the pulsation pair-instability supernova (PPSN) model from Belczynski et al. (2016).By including or excluding fallback we control the retention of BHs in our simulations.

Mock photometry
To convert snapshots from the -body models to photometric data for different filters used in observations, we use the code galevnb (Pang et al. 2016), which selects corresponding spectral templates from the library of Lejeune et al. (1997Lejeune et al. ( , 1998) ) according to the fundamental stellar properties, such as stellar mass, temperature, luminosity and metallicity from -body simulations.By convolving the spectra with the filter response curve from a given filter, we obtain the observational magnitudes of specific filters of main-stream telescopes, such as Hubble Space Telescope (HST) and the future Chinese Survey Space Telescope (CSST) for individual stars in the body models.In this way, we produce mock observations for -body models, which allows a direction comparison with observational data.This is useful to compare the density or surface brightness profiles, unresolved binaries and stellar mass functions between observations and the models.In this study, the line-of-sight velocity of unresolved binaries is calculated using the Johnson I-band filter (as described in Section 3.2.4).For creating the color-magnitude diagram, we employ the HST F555W and F814W filters, along with the CSST g and i fil-ters.To convert luminosity to mass for unresolved binaries, we utilize the HST F555W filter.Further details can be found in Section 3.4.

Observational data
To validate our -body model and ensure its accuracy in reproducing the surface number density Σ() and mass function of Pal 5, we compare it with observational data.We utilize the data from Ibata et al. (2017) for the surface number density and the masses of stars obtained from two HST observations with Program IDs 6788 (PI: Smith;Grillmair & Smith 2001) and 14535 (PI: Kuepper) as reported in Baumgardt et al. (2023).
The observed surface number density Σ() encompasses stars with g-band magnitudes ranging from 19 to 23, with photometry obtained from the Canada-France-Hawaii Telescope.The corresponding mass range of these stars is 0.625 to 0.815  ⊙ , determined using the magnitude-mass conversion provided by G21.
Regarding the masses of stars derived from the HST data, Baumgardt et al. (2023) employed Dartmouth isochrones to fit the CMDs of the clusters and employed them to convert magnitudes into masses.Further details can be found in their work.

Star cluster models
To reproduce Pal 5's observed surface density and present-day position in the Galaxy, we generate the initial conditions of -body models by referring to the wBH-1 and noBH-1 models in G21, which have the closest property to the observational data assuming Pal 5 contains a cluster of BHs and no BH, respectively.
For the wBH-1 model, natal kick velocities of BHs after supernovae are affected by the material fallback from Fryer et al. (2012).A large fraction of BHs are retained in the clusters and finally sink to the centre via dynamical friction.The existence of a BH subsystem can significantly affect the structure and evolution of star clusters.As a result, the cluster has a loose core of luminous stars.The wBH-1 model has an initial half-mass radius,  h,0 = 5.85 pc, and an initial number of stars,  0 = 2.1 × 10 5 .
In contrast, the noBH-1 model assumes BHs have the same high kick velocities as neutron stars and almost none are retained after supernova explosions.Without BHs, the core collapse of luminous stars result in a dense core.In order to reproduce the observed surface brightness profile, G21 find that the cluster must therefore have had a much lower density initially.Thus, for the noBH-1 model,  h,0 = 14 pc and  0 = 3.5 × 10 5 .
We conducted five -body models with varying setups of primordial binaries and the presence of BHs.The initial conditions for these five models are summarized in Table 1.We assigned labels to the models to indicate the existence of primordial binaries and BHs.
For BH treatment, models with the label "BH" refer to the wBH-1 model from G21, where the mass fallback scaling for kick velocities is applied so that a part of the BHs has low kick velocities and stays in the clusters.They also have the same  0 and  h,0 as those of the noBH-1 model.
Models with the label "noBH" refer to the noBH-1 model from G21.In these models, all BHs have high kick velocities similar to the neutron stars after asymmetric supernovae.The velocity distribution follows a (1D) Maxwellian distribution with a dispersion of 265 km/s.As a result, we found no BHs are retained in our noBH models.
The prefix "noBin" and "Bin" represent without and with primordial binaries, respectively.For "Bin" models, all stars are in binaries initially.For massive binaries with the component mass > 5 M ⊙ , except the Bin-noBH-F model, all other "Bin" models have the period and mass ratio distributions follow the observational constraints of OB binaries from Sana et al. (2012).
For low-mass binaries, except the Bin-BH-Alt model, all other "Bin" models assume the properties of primordial binaries following the model from Kroupa (1995a,b) and Belloni et al. (2017) (naming as Kroupa binary model).The orbital parameters of this model are derived from the inverse dynamical population synthesis of binaries in the Galactic field.This model assumes an universal property of primordial binaries and all stars forming in star clusters.In addition, a correction of the period and eccentricity distributions from Belloni et al. (2017) is included to better fit the observational data of GCs.
For the Bin-BH-Alt model, we assume a different setup of lowmass primordial binaries (referred to as FlatLog model) as a comparison with the Kroupa binary model.The semi-major axes follow a flat distribution in the logarithmic scale where the minimum and maximum value are 3 solar radius and 2 pc, respectively.The eccentricity and mass ratio distributions are the same as those of the Kroupa binary model.
The period and eccentricity distributions are shown in Figure 2.For both binary models, the initial distribution of periods covers a wide region with 9 orders of magnitudes.The initial eccentricities exhibit a sharp peak at  = 0 and a broader peak at  = 0.8, respectively.All binaries with peri-centre separation less than the sum of the stellar radii of the two components are excluded.Thus, an empty region is visible in the period-eccentricity distribution of Figure 2. In addition, the eccentricity distributions of the Kroupa and FlatLog are different after adjustment.
These binary setups cover a wide range of binary orbital periods, where a large fraction of binaries are unstable in the cluster environment.After a short time (about one crossing time), the binary fraction significantly reduces.Referring to Pal 5, the binary fraction of our setup may be overestimated.The benefit is that we can investigate how long-term dynamical evolution of the clusters with and without BHs affect both the tight and wide binaries.
The Bin-noBH-F model has the same  0 and  h,0 as those in the noBH-1 model.However, after finishing the simulation, we found that the Bin-noBH-F model cannot reproduce the final structure of the noBH-1 model at 11.5 Gyr and it has sufferred complete tidal disruption before 10 Gyr.The suffix "F" in the name of the model indicates that this is a failed model.Thus, we conducted another model "Bin-noBH" by reducing  h,0 to 13.2 pc.This small modification results in a cluster similar to Pal 5 after 11.5 Gyr.
In addition, we excluded massive binaries in the Bin-noBH-F model to prevent non-supernovae BH formation in a binary, but we observed that such events did not occur.Therefore, in the Bin-noBH model, we added the Sana distribution to massive binaries to ensure consistency with the Bin-BH models.
The common setup for all models is also summarized in Table 1.All models were evolved for a duration of 12.0 Gyr.At 11.5 Gyr, the clusters are located at the same Galactic position as Pal 5.However, since the model did not precisely reproduce the surface number density of Pal5, we continue to evolve the cluster further to determine the age (referred to as  mat ) when the model matches the observation more closely, as detailed in Section 3.1.5.We assumed a spherically symmetric Plummer profile (Plummer 1911) with no primordial mass segregation.The initial mass function (IMF) of stars followed the two-component power-law shape described by Kroupa (2001).We adopted the same mass range of 0.1 − 100 ⊙ as used in G21, and the power-law indices () and mass ranges are described as: In this study, we adopted a cluster metallicity of  = 0.0006, which is consistent with the value reported in Smith et al. (2002) of [Fe/H] ≈ −1.4 dex for Pal 5.The initial star cluster models were generated using the updated version (Wang et al. 2019) of the mcluster code (Küpper et al. 2011).This update includes the implementation of the Kroupa binary model generator, as shown in Figure 2.

Structural evolution
First, we present the evolution of the cluster structure and compare our results to the models from G21 and the observational data.Generally, although the existence of binaries does not significantly affect the structural evolution, the small difference can be amplified by the Galactic tidal field and result in early dissolution of the Bin-noBH-F model.In addition, the existence of primordial binaries reduces the BH populations and results in shorter relaxation times in the early evolution.The stochastic formation of BBHs also affects the expansion of the cluster and eventually influences the disruption of the cluster.The surface number density of -body models roughly agree with observations with a larger central density.

Half-mass relaxation time
The two-body relaxation time is an important timescale of stellar dynamics, which reflects the speed of changes in the density and mass segregation of a cluster and its tidal dissolution.The onecomponent half-mass relaxation time ( rh1 ) defined in Spitzer (1987) has the form as where  is number of stars,  h is the half-mass radius,  is the average mass of stars,  is the gravitational constant, and ln Λ is the Comlumb logarithm.When BHs exist, the binary heating is dominated by BBHs,  rh1 leads to an underestimation of the relaxation timescale of the system.Wang (2020) found that a proper two-component relaxation time ( rh ) can be obtained by dividing a correction factor , defined as and where the suffixes 1 and 2 represent the quantities for non-BH and BH components, respectively.Figure 3 illustrates the evolution of  rh and .The three BH models exhibit significantly shorter  rh compared to the noBH models.During the first 100 Myr, the noBin-BH model displays a longer  rh compared to the Bin-BH and Bin-BH-Alt models because the Bin models treat binaries as single objects when calculating  rh .Consequently, the Bin-BH and Bin-BH-Alt models experience relatively faster expansion of  h and faster mass segregation of BHs (see Section 3.1.2).Subsequently, the trend reverses, and the  rh of the noBin-BH model becomes shorter than that of the Bin-BH and Bin-BH-Alt models due to the difference in the number of BHs (see Section 3.1.3).As a result, the  h of the noBin-BH model expands faster than that of the other two models.After 8 Gyr, the  rh of all three BH models starts to decrease due to mass loss via tidal evaporation.
The values of  for the BH models exceed 5, indicating that BHs significantly impact the relaxation process of the clusters.Further discussion of  h is provided in Section 3.1.2.
In contrast, the two noBH models exhibit much longer  rh .There is a rapid increase in  rh during the first 100 Myr, primarily due to the strong stellar winds from massive stars and the escape of BHs.Consequently, although the morphology appears similar at 11.5 Gyr for models with and without BHs, the relaxation processes differ significantly.These differences can lead to variations in the properties of binaries.In Section 3.2, we analyze the impact of these differences and discuss their implications for binary systems.It is important to note that assuming  = 1 for the noBH models is not accurate, as there is still an order of magnitude difference between the minimum and maximum masses of stars.

Half-mass radius
Figure 4 illustrates the evolution of  h for all models, including the ones from G21 for comparison.We observe that the presence of primordial binaries has a weak impact on the evolution of  h , consistent with the theoretical findings of Wang et al. (2022).When BHs exist, the long-term structural evolution of star clusters is primarily controlled by binary heating driven by the dynamical interactions between BBHs and the surrounding objects at the cluster center.The majority of primordial binaries have much smaller masses compared to BBHs, and therefore have a negligible impact on the binary heating until most BHs have escaped from the cluster.A small subset of massive primordial binaries can eventually evolve into BBHs.However, even in the absence of these massive binaries, a star cluster can generate BBHs through chaotic three-body interactions when the central density of the cluster reaches a threshold after the core collapse of BHs (see Section 3.1.4).Consequently, we only observe minor differences of  h between the Bin-BH, Bin-BH-Alt, and wBH-1 models during the first 10 Gyr of evolution.This can be explained by the differences in relaxation times ( rh ) discussed in Section 3.1.1.The galactic potential also affects  h , but since all models share the same orbit, the influence is similar.However, after 10 Gyr, the Bin-BH-Alt model exhibits a similar  h to that of the wBH-1 model, but its  h shows significant variations, indicating an energy imbalance and the onset of a disruptive tidal phase.In contrast, both the Bin-BH and wBH-1 models remain stable until 12 Gyr.This differing behavior is attributed to stochastic BBH heating, as explained in Section 3.1.4.
The BH models with binaries (Bin-BH) and without binaries (noBin-BH) exhibit different timescales for the mass segregation of black holes, as indicated by the initial rapid contraction of  h,BH .In the Bin-BH model,  h,BH undergoes faster contraction during the early stages of evolution compared to the noBin-BH model.This disparity can be attributed to the difference in  rh , as the timescale for mass segregation is proportional to  rh .
When comparing the noBH models with binaries (Bin-noBH-F) and the model from G21 without binaries (noBH-1), significant differences in the evolution of  h emerge after 8 Gyr.The Bin-noBH-F model experiences tidal disruption at around 9 Gyr, whereas the noBH-1 model survives until 11.5 Gyr.G21 noted that the final properties of the noBH models are more sensitive to changes in the initial conditions, and in fact argued that this 'fine tuning' problem disfavours the noBH scenario.An offset of  h,0 needs to be introduced in the Bin-noBH model to achieve consistent  h at 11.5 Gyr.
Two factors may explain the need for this offset.Firstly, in the absence of BHs, binary heating is primarily generated by low-mass binaries.Consequently, the influence of primordial binaries is more pronounced compared to models with BHs.Secondly, due to the larger  h,0 , the cluster becomes more sensitive to the galactic tide.The presence of primordial binaries affects the relaxation time of the system, as the dynamical effect of tight binaries is equivalent to that of single objects, resulting in a shorter relaxation time for the sys-tem.Consequently, the system dissolves faster, necessitating a denser initial cluster to allow the cluster's survival, as seen in the noBH-1 model.Additionally, the differences caused by the stochastic scatter of  h resulting from the random seeds used to generate the initial conditions may also be amplified by the galactic tide, contributing to the divergent evolution.

Mass loss
The upper panels of Figure 5 show the evolution of the total mass ( ()) of our models.Data of the wBH-1 and the noBH-1 from G21 are also shown as references.The mass loss has two channels: wind mass loss driven by stellar evolution and escapers via stellar dynamics of star clusters.To have a consistent definition of , all models use the same criterion to select escapers.First, we calculate the bound energy of stars and centre-of-the-mass of binaries without external potential and then select escapers with energy >0.
Here we compare the three cases: For models with no primordial binary and with BHs,  () of our noBin-BH model agrees with the wBH-1 model from G21.The final mass of the noBin-BH model at 11.5 Gyr is slightly larger than that of the wBH-1 model.
For models with primordial binaries and with BHs, compared to the wBH-1 model, the Bin-BH and the Bin-BH-Alt models lose mass faster during the first few hundred Myr, but mass loss of the Bin-BH model becomes slower near the end of the simulation.Finally, the Bin-BH and the wBH-1 models agree with each other, while the Bin-BH-Alt model dissolves after about 11 Gyr.
For models with no BHs, the Bin-noBH-F model with primordial binaries loses mass faster than the noBH-1 model with no binaries.The Bin-noBH model, with a smaller  h,0 , experiences a relatively slower mass loss, and its  () remains slightly above that of the noBH-1 model at 11.5 Gyr.In general, the evolution of  () and  h are similar for all three cases.

Black holes
BHs significantly affect the long-term dynamical evolution.We investigate the mass fraction of BHs (  BH ) and the bound mass of BHs ( BH ) in Figure 5.The evolution of  BH in the noBin-BH and the wBH-1 models agree with each other in the first 8 Gyr.Then,  BH increases more slowly in the noBin-BH model and is half that in the wBH-1 model at 11.5 Gyr. BH of the noBin-BH model is slightly smaller than that of the wBH-1 model initially and such a difference is inherited in the long-term evolution.Finally, as a large fraction of stars escape, such initial differences lead to a large difference of  BH at the end.
For the Bin-BH and the Bin-BH-Alt models,  BH are significantly smaller than that of the noBin-BH model during the early evolution.This difference is due to the stellar evolution of massive binaries.Based on the orbital parameters of binaries from Sana et al. (2012), the progenitors of BHs (massive stars) are all in binaries.A fraction of the tight binaries suffers mass transfer and mergers.The BHs formed from these binaries can have different distribution of masses.The maximum  BH of the Bin-BH model is about 250  ⊙ less than that of the noBin-BH model.Then, after the mass segregation of BHs (a few hundreds Myr), binary heating of BBHs start to kick out BHs from the cluster, and result in larger difference of  BH during the long-term evolution.Although the Bin-BH (Bin-BH-Alt) and the noBin-BH models show a large difference of  BH , their evolution of  and  h is similar before 10 Gyr.This was also observed in Wang et al. (2022).The evolution of the semi-major axes () of BBHs reflects both binary heating and mergers driven by gravitational wave (GW) radiation.Figure 6 provides a comparison of this evolution for the three BH models.Despite the absence of primordial binaries in the noBin-BH model, we can still observe the formation of BBHs and their orbital contraction.The frequency of BBH formation and the overall trend of  are similar for all three models, except that the two models with primordial binaries exhibit a higher number of BBHs formed from these binaries during the first 1000 Gyr.Some of these BBHs with  < 1 AU undergo orbital shrinking due to GW radiation, ultimately merging to form more massive BHs.These newly formed BHs lead to the creation of massive BBHs with masses exceeding 100  ⊙ .The presence of these massive BBHs can have a substantial impact on the evolution of the star cluster, influencing its dynamical and structural properties.
In particular, for the Bin-BH-Alt model, the formation of a massive BBH around 8 Gyr coincides with a faster expansion of  h compared to the Bin-BH model, ultimately leading to an earlier disruption of the Bin-BH-Alt model.Hence, the divergent evolution of the Bin-BH and Bin-BH-Alt models after 8 Gyr is attributed to the stochastic formation of BBHs.
It is important to note that our models do not account for the high-velocity kicks experienced by newly formed black holes due to asymmetric GW radiation following mergers.Therefore, the formation of such massive BBHs might not be as common as our models suggest.Consequently, the stochastic effect of massive BBH heating could be overestimated in our cases.

Surface number density profiles
The determination of  h and  relies on the selection criteria for identifying cluster members.When comparing the -body models with observational data from Pal 5, it is challenging to use the exact same selection criterion for both.A more appropriate approach is to compare the surface number density (Σ()), where  represents the angular distance from the cluster center in the International Celestial Reference System (ICRS).
Figure 7 illustrates the Σ() profiles for our -body models and the observational data of Pal 5 obtained from Ibata et al. (2017).To ensure consistency with the observations, only main-sequence stars with masses ranging from 0.625 ⊙ to 0.815 ⊙ are considered in the -body data (see G21 for details).
No stars are removed during the simulation, allowing for the tracking of the tidal tail evolution.The centre-of-mass position of the star clusters in the Galaxy at exactly 11.5 Gyr does not perfectly align with that of Pal 5.This is due to the long-term evolution of star cluster, where the center of the cluster drifts as a result of asymmetric mass loss due to stellar winds, supernovae, and the escape of stars.Therefore, we select snapshots from the simulations that have the closest centre-of-mass distance to that of Pal 5 whenever a comparison is required in the subsequent analysis.We then correct the positions and velocities of the stars by applying the offset between the centre-of-mass of the -body models and the observational data.The results of this correction are presented in the upper panel of Figure 7. Due to the complete disruption of the Bin-noBH-F model, it is not possible to determine the centre-of-mass position for this particular model.Therefore, it is excluded from some analysis and comparisons.
The vertical lines in Figure 7, representing the half surface number radii ( hn ), indicate that all models except the Bin-BH-Alt model are more centrally concentrated than the observed Pal 5.In Figure 5, it is shown that these models retain more mass at 11.5 Gyr compared to the models presented in G21.
The Bin-noBH and Bin-BH models exhibit similar Σ() profiles, but this similarity is coincidental since they had different initial density profiles and evolved in opposite ways, as demonstrated in Figure 4.
Given the time-consuming nature of the simulations, it is challenging to precisely reproduce the models of G21 and the observational data.To enhance the comparison with the observational data, we selected snapshots at different ages that match the observed Σ() profile.These results are displayed in the bottom panel of Figure 7.Although the tidal streams differ substantially, we can still compare  the internal properties of binaries and mass functions using these snapshots.

Binding energy of binaries
While the BH and noBH models may exhibit a similar Σ() profile, as demonstrated in Figure 7, their relaxation processes differ.This discrepancy can lead to different properties of binaries at 11.5 Gyr.
In star clusters, perturbations from incoming objects can significantly alter the orbits of binaries.According to the Heggie (1975)-Hills (1975) law, wide or soft binaries are prone to disruption after experiencing a few close encounters with intruding objects.Conversely, tight or hard binaries tend to become even tighter after these encounters.
The hard-soft boundary of binding energy ( hs ) at the distance to the cluster center () is determined by the local velocity dispersion: where 0.5⟨ 2 ⟩ is the average kinetic energy of stars and binaries at , and  is the velocity.The hard-soft boundary of binaries evolves as the structure of the cluster changes over time.Initially, during the first 100 Myr of star cluster evolution, there is a rapid reduction in the hard-soft boundary.This is due to the expansion of  h caused by the strong stellar wind mass loss from massive stars, as shown in Figure 4.After 100 Myr, the evolution of  h slows down, and the hard-soft boundary,  hs , evolves more gradually.The Bin-BH and Bin-noBH models have different initial  hs () curves as shown in Figure 4, but their final  hs () curves at 11.5 Gyr converge to a similar shape.This indicates that the distribution of binary binding energy at 11.5 Gyr may reflect the different evolutionary histories of  hs .
To further analyze the distribution of binary binding energy, Figure 8 presents a comparison of the contour plot of  b versus  at approximately 11.5 Gyr for the Bin-BH and Bin-noBH models.Across a wide range of  values, spanning from the center of the cluster to the distant tidal tail, two distinct peaks can be observed.The first peak, located around 10-30 pc, represents the population of binaries inside the cluster.The second peak, with  > 3000 pc, corresponds to binaries that have escaped from the cluster and are distributed along the tidal tail.
We focus on the discussion of binaries within the cluster and examine the hard-soft boundaries,  hs (), at three different ages: 0 Myr, 100 Myr, and 11.5 Gyr.These boundaries are plotted as reference curves.To calculate  hs (), we divide the cluster into 10 radial bins, ensuring an equal number of objects per bin.Binaries are treated as unresolved objects in this analysis.The maximum value of  is set to be at 90% of the Lagrangian radius, providing a radial range that reflects the cluster's size at the three ages.
The results show that  hs () does not exhibit strong variations along .The two models, Bin-BH and Bin-noBH, have similar  hs () curves, except for an offset in the radial region at 0 Myr and 100 Myr.The peak of  b falls between the  hs () curves at 100 Myr and 11.5 Gyr.This suggests that during the first 100 Myr, not all soft binaries with  b <  hs are immediately disrupted, and many of them can survive and become hard binaries by 11.5 Gyr.
Therefore, the final distribution of  b does not clearly reflect the initial conditions of the two models, as anticipated by G21.However, the Bin-noBH model has a relatively larger number of binaries compared to the Bin-BH model.This difference suggests that the overall rate of binary disruption depends on the evolutionary history of the cluster density.

Period distribution
To analyze the binary disruption rate in relation to cluster dynamics, we examine the period distributions normalized by the bound mass of the cluster ( hb ) for three models: Bin-BH, Bin-noBH, and Bin-BH-Alt, as depicted in Figure 9.The period distributions at the initial phase (0 Gyr) and the median age (5 Gyr) are compared.
In the Bin-BH and Bin-noBH models, the initial period distributions are the same, but they exhibit different density profiles.At 5 Gyr, the Bin-noBH model retains more wide binaries compared to the Bin-BH model.The hard-soft boundaries of periods, estimated for stars within  h , do not exhibit significant differences between the two models.However, the peak of the period distribution in the Bin-BH model is closer to the hard-soft boundary at zero age, whereas in the Bin-noBH model, it aligns with the boundary at 5 Gyr.This disparity suggests that the disruption rate of binaries is not solely determined by the hard-soft boundary.During long-term evolution, the Bin-BH model, which is denser and contains BH subsystems, experiences a higher rate of disruption for wide binaries, resulting in the peak of the period distribution being closer to the boundary.In contrast, the Bin-noBH model preserves more wide binaries, and the peak of the period distribution reflects the boundary at 5 Gyr for the cluster.
Comparing the Bin-BH and Bin-BH-Alt models, they share a similar density evolution but differ in the assumptions of their primordial binaries.The ratio of  hb at 5 Gyr to the initial phase,  hb (5 Gyr)/ hb (0), exhibits an identical trend for both models.This finding implies that the binary disruption is not highly sensitive to the assumption of the initial period distribution.Consequently, it is possible to infer the initial binary properties through inverse derivation if the evolution history of the cluster density is known (see Kroupa 1995a;Marks et al. 2011;Marks & Kroupa 2012).Moreover, by utilizing the derived ratio, we can extrapolate the evolution of the period distribution of binaries for any arbitrary assumption regarding the primordial binary populations.This provides a valuable tool for understanding the long-term dynamical evolution of binary systems within star clusters and can aid in studying the impact of different initial binary properties on the binary disruption rate and cluster dynamics.

Radial distribution
Figure 10 compares the radial distribution of the binary fraction (  bin ) for the Bin-BH and Bin-noBH models at 11.5 Gyr.
In the upper panel, the real  bin is plotted as a function of the 3D radial distance from the cluster center.Both models exhibit a similar trend, with a systematic offset of  bin along .The central region of the cluster shows a higher  bin compared to the outer halo.At the distant tail of the cluster,  bin experiences a significant increase.This can be attributed to binaries that escaped from the cluster during the early stages of evolution, as they suffer fewer dynamical perturbations and have a higher chance of survival.
The lower panel of Figure 10 presents the predicted observed binary fraction as a function of projected distance.To identify binaries from the color-magnitude diagram, we assume that unresolved binaries with B-band magnitudes between 20.5 and 23 mag and a mass ratio above 0.6 can be detected.The B-band magnitudes for stars are generated by using galevnb.Notably,  bin (obs) for both models is nearly identical within a projected distance up to 30 arcmin, unlike the real  bin for all binaries.The observed binary fraction  bin (obs) falls in the range of 0.2 to 0.3.

Half-year evolution of line-of-sight velocities
With high-resolution multi-epoch spectroscopic observations, it is possible to identify binaries by comparing the line-of-sight velocity changes (|Δ LOS |) over a span of approximately six months.
The line-of-sight velocity  LOS of an unresolved binary is the combination of two  LOS of two components and is dominated by the brighter component.Thus, the |Δ LOS | values exhibit considerable variation during the multiple epochs of observation.These variations are determined by the periods, eccentricities, inclinations, and orbital phases of the binaries.Notably, larger variations are observed for short-period binaries, which could potentially aid in distinguishing these binaries from other effects that cause changes in velocity.The baseline of approximately half a year is sensitive to a maximum period of ∼ 10 4 days.
We estimate  LOS of binaries by taking the I-band flux-weighted average of the  LOS of the two components.In Figure 11, we present the |Δ LOS | versus period plot for observable unresolved binaries with |Δ LOS | > 0.3 km/s and  < 10 arcmin after multiple epochs, respectively.We specifically select binaries with at least one bright (post-main-sequence) star component, and some binaries include white dwarfs.These bright stars have a luminosity in the HST 555 filter brighter than 20 mag.The three models (Bin-noBH, Bin-BH, and Bin-BH-Alt) exhibit observable binaries across a wide range of period distributions, spanning from 1 to 10 4 days.The snapshots at  mat (see the bottom panel of Figure 7) are chosen as the first epoch of observation.The choices of time intervals between epochs were chosen to be roughly equal space in half a year time interval, and the exact values are defined by the time step algorithm of the petar code.The number of detectable binaries is similar for all three models, with the Bin-noBH model exhibiting slightly more binaries with periods above 3000 days.This trend aligns with the period distributions shown in Figure 9, although some stochastic scatter may be present.
To assess the completeness of detectable binaries via multi-epoch observations of |Δ LOS |, we compare the number counts of detectable binaries and all bright binaries as a function of periods, as shown in Figure 12.For all models, periods up to 10 4 days are detectable and all binaries with periods below 10 3 days can be detected with multiple epochs.From Figure 11, one binary in the Bin-BH model with a period between 10 3 − 10 4 days has only one epoch that shows |Δ LOS | > 0.3 km/s.A few binaries above 10 3 days in the Bin-noBH models have epochs where |Δ LOS | < 0.3 km/s, indicating that they might be missed if the observational epochs are limited to two.
The observed  LOS of unresolved binaries does not represent the  LOS of the center-of-mass of the binaries, which complicates the determination of the physically useful line-of-sight velocity dispersion ( LOS ).A complete sample of detectable bright binaries with periods below 10 4 days can mitigate this effect and significantly improve the determination of ( LOS ).When binaries are detectable from multiepoch observations, we can exclude them from the computation of   dimensional velocity dispersion  1D within  h , assuming a virial equilibrium state of the cluster: This normalization allows us to account for any differences in the overall dynamical state of the clusters and facilitates a more meaningful comparison of the  LOS .
The presence of BHs affects the  LOS in the cluster center.To illustrate the difference between models with and without BHs, we calculate the  LOS of single stars within a projected distance of  < 3 arcmin ( LOS,S,hn ), which corresponds to the  hn (17 pc).All three models exhibit similar values of  LOS,S,hn .Additionally, the  LOS values of single stars within a projected distance of  < 10 arcmin (58pc), which includes stars outside the effective radius of the cluster, are similar to  LOS,S,hn , except the Bin-noBH model, which has a lower value.
Since the normalization factor  1D is different for the three models, and the observation cannot directly obtain  and  h , the difference in the observed estimates of  LOS for the three models may be larger than what we found in our simulations.This should be taken into consideration when interpreting the results and comparing them with observations.
The sample that includes all bright singles and binaries exhibits much larger dispersion values ( LOS,SB ) than the values ( LOS,S ) of the sample containing only singles.By excluding detectable binaries, the values ( LOS,SCB ) are significantly lower than  LOS,SB , roughly 1.5-2 times of  LOS,S .This procedure helps to obtain more accurate estimates of  LOS .

Binaries with BHs
The Bin-BH model at 11.5 Gyr exhibits several binaries which contain one or two BHs (BwBHs), as depicted in Figure17.It is important to investigate whether these BwBHs can be detected, serving as evidence for the existence of BHs.Table 3 provides a summary of the parameters for these binaries, which include three types: BBHs, BH with MS (BH-MS), and BH with WD (BH-WD).Other types of BH-star binaries are not detected.
The presence of BBHs has also been illustrated in Figure 6, with the possibility of some being detected by GW detectors.Three BBHs are inside the clusters and the other three distribute in the tidal stream.
An interacting BwBH that contains an accreting BH primary and a non-BH secondary star is particularly interesting as a potential Xray or radio source that could be detected, providing evidence for the presence of BHs in Pal 5. Unfortunately, there is no BwBH that contains a bright post-main sequence star at 11.5 Gyr, only a few BH-MS and BH-WD exist.
We calculate the Roche lobe radius using Equation 53from Eggleton (1983); Hurley et al. (2002), with the semi-major axis replaced by the peri-center distance : where  =  2 / 1 .The original formula assumes a circular orbit, which misses the eccentric binaries where the accretion may occur at the peri-center separation.To account for this, we use the pericenter distance  instead.When the stellar radius of the secondary star ( 2 ) is greater than or equal to the Roche lobe radius ( RL2 ), the secondary star fills its Roche lobe, and the accretion process might result in observable radiation.
The  2 / RL2 values of BH-MS binaries in our models are below 10 −3 , indicating that no accretion occurs in these cases.The BH-WD binaries have the potential to become ultraluminous X-ray sources (ULXs).Detailed studies of the dynamical formation scenarios for these ULXs in globular cluster environments have been conducted by Ivanova et al. (2010).One BH-WD binary in our simulations has a period of 2.5 days and a peri-center distance () of 2 ⊙ , located ∼ 4.5 pc away from the cluster center.The ratio  2 / RL2 is approximately ∼ 0.04, which does not yet reach the criterion for accretion.
In our investigation of the BH-MS binaries, we have discovered that their formation occurs through a similar dynamical channel.The MS star originates from a primordial binary of two MS stars (MS-MS).The BH originates from a primordial binary of two massive stars, which forms a BBH.The formation process of the BH-MS binaries in the Bin-BH model involves several steps: (i) The BBH undergoes several interactions with other BHs in the cluster.
(ii) After one of the BHs escapes from the cluster following a strong interaction with an intruder, it becomes a single BH.
(iii) This single BH eventually encounters the MS-MS binary and participates in a binary exchange event.
(iv) As a result of the binary exchange, the BH joins the MS-MS binary, forming the BH-MS binary.
The described process is visually illustrated in Figure 14.The dynamical formation of BH-MS binaries in star clusters have been discussed in several works (Kremer et al. 2018;Di Carlo et al. 2023;Rastello et al. 2023;Tanikawa et al. 2023).
Although no observable events from interacting BwBH occur at 11.5 Gyr, we can estimate the frequency of such events by collecting the interacting BwBHs recorded in the evolution of star clusters.The criterion to select interacting BwBHs are  2 / RL2 ≥ 1. Events that occurred in the first 100 Myr are excluded, as they mostly involve primordial binaries that are not significantly affected by stellar dynamics.The results are summarized in Table 4.
The Bin-BH and Bin-BH-Alt models have a dozen of such interacting BwBHs, including both primordial and dynamically formed BwBHs.The dynamically formed BwBHs contribute to approximately half of the interacting BwBHs.The secondary stars involved in these BwBHs include several types, with one being BH-NS, which can trigger a GW merger.
The Bin-noBH model also includes 5 events, all of which consist of primordial binaries.Among these events, four are BH-MS binaries, and one is a BH-NS binary.Despite the high supernovae kick velocities in the Bin-noBH model, these binaries were strongly bound before the supernovae, and the random natal kick did not disrupt the binaries.Instead, the binaries escaped from the cluster after the kick.
In general, the formation rate of an interacting BwBH is estimated to be about one per 2 Gyr.Therefore, the possibility of detecting an interacting BwBH in the present-day Pal5 is practically zero.
The noBin-BH and Bin-noBH-F models do not exhibit any interacting BwBH events, and thus, they are not included in the table.One common feature of these two models is the absence of massive primordial binaries, which is different from all other models that have OB binary properties from Sana et al. (2012).As a result, the possibility of dynamical formation of BwBHs is also low in these models.One important channel for the formation of interacting BwBHs is through the dynamical exchange of binary components after a close encounter between a BH and a binary.The lack of primordial binaries in these models suppresses this formation channel.
Multi-epoch observations of |Δ LOS | can also be used to detect non-interacting BwBHs.For instance, utilizing multi-epoch MUSE spectroscopy, Giesers et al. (2018Giesers et al. ( , 2019) ) discovered three BwBHs in NGC3201.The stellar companions in these BwBHs have mass values of 0.6 − 0.8  ⊙ .The four BH-MS binaries in the Bin-BH model at 11.5 Gyr have comparable companion masses.Therefore, it is possible to detect BHs in Pal 5 via multi-epoch observations of |Δ LOS |.However, due to the long periods of these binaries, a long-term observation plan (several years) is needed to accurately constrain the masses of the BHs.Despite the fact that these binaries are not  LOS variable over a short baseline of a few months, they may still be found: they should appear as member stars according to their position in the CMD, parallax and propor motion, but they have a large  LOS offset.A solar-type star orbiting a 15 M ⊙ BH with a 10 4 d period has an orbital velocity of ∼ 25 km/s.This predicted signal is worth looking for.

Color-magnitude diagram
By utilizing the galevnb code, we can convert our simulation data into mock photometry.As an example, we present the colormagnitude diagram (CMD) of the Bin-BH model at 11.5 Gyr, using HST 555 and 814 filters, and CSST  and  −  filters (Fig- ure15).
In the CSST filters, we observe binary stars distributed between the MS and WD sequence.These binaries consist of a WD and a lowmass main sequence star (LMS).Similar features in the CMD have been seen in -body simulations by Pang et al. (2022) (see figure 5 in Pang et al. 2022) 1 .In these binary systems, the luminosity is mainly dominated by the WD, as both components have very similar masses.They are considered as candidates for cataclysmic variable (CV) stars.
The CSST -band magnitudes of WD and CV are below 26 mag, while the corresponding HST F555W magnitudes are above 26 mag.Therefore, CSST has the advantage of potentially detecting many WD and CV candidates in Pal 5.
We also highlight the BH-MS binaries shown in Table 3.Among them, three have the HST F555W magnitude below 21 mag and the CSST -band magnitude below 16 mag.If the multi-epoch spectroscopy observation can reach this magnitude limit, it is possible to detect these binaries via the observation of |Δ LOS |.

Mass functions
The present-day mass function of a star cluster is influenced by various factors, including the IMF, mass segregation, and tidal evaporation.To investigate the impact of primordial binaries and black holes (BHs) on the mass function, we compare the mass functions of our -body models with the observed ones.In order to make a meaningful comparison with the observed data, we select snapshots from our models that closely match the observed surface number density profile (Σ()), as shown in the lower panel of Figure 7.
It is important to consider the resolution limitations when comparing with observations.The widest binary in our models has a semi-major axis of approximately 1.8 × 10 4 AU.Given the distance to Pal 5, a spatial resolution of less than 1 ′′ is required to resolve this binary.The best resolution achievable by HST is around 0.05 ′′ , which means that only a small fraction of wide binaries with periods above 1.4 × 10 7 days can potentially be resolved.Therefore, we assume that most binaries remain unresolved in observations and calculate their magnitudes by summing the fluxes of their two components.Figure 15 shows the color-magnitude diagram (CMD) of unresolved binaries, which appear redder and brighter compared to the single stars.
To investigate this effect, we compare the (actual) total masses ( tot ) of binaries with the masses converted from their F555W-band magnitudes ( obs ).
For main-sequence binaries, we calculate the absolute F555Wband flux and then determine the mass of a single star that has the closest flux value, which serves as the converted mass  obs .The comparison between  tot and  obs is depicted in Figure 16.
Table 3.The parameters of BwBHs for the Bin-BH model at 11.5 Gyr. 1 and  2 denote the masses of the primary and secondary components, respectively;  represents the peri-center distance;  2 / RL2 indicates the secondary stellar radius relative to the Roche lobe overflow radius; and  represents the distance of the binary from the cluster center.

Type
The difference between  tot and  obs is highly sensitive to the mass ratio  =  1 / 2 and luminosity ratio as well.Here, the mass ratio  is defined as the minimum mass divided by the maximum mass of the two components in a binary.A higher  leads to a larger difference between the  tot and  obs values.Consequently, the  obs of equal-mass unresolved binaries can be significantly lower than their true  tot .
Furthermore, for binaries with the lowest  values, there is a systematic offset between  tot and  obs .As a result, if unresolved main-sequence binaries cannot be distinguished from single stars, the total masses of all these binaries would be underestimated.
The offset between  tot and  obs is determined by the minimum .There is a nonlinear relation between stellar luminosity () and mass ().For MS stars in the mass range of 0.3-0.8 ⊙ ,  ∝  4 , and thus, we can roughly estimate the relation between the total binary mass ( tot ) and the binary mass used in the mass function estimation  3. The left panel corresponds to the HST F555W-F814W and F555W filters, while the right two panels correspond to the CSST g-i and g, and u-y and u filters, respectively.( obs ) as follows: In our model, the minimum  is about 0.12, which corresponds to a maximum  obs / tot ≈ 0.93.
To compute the mass functions, we collect stars within the same observational fields used by the HST observation from the Smith field (Grillmair & Smith 2001) and the Kuepper field (unpublished; reported in Baumgardt et al. 2023), as shown in Figure 17.The center position of the star cluster model is defined as the centre-of-mass of allowing us to obtain mass functions as a function of radial distance.The intersection between the two observational fields and the three radial bins (referred to "Field" regions) are used for selecting samples of stars.
It's important to note that due to the limited observational coverage and stochastic scatter, the comparison between the observed and modeled mass functions may be affected.To improve statistical robustness, we also select stars for measuring the mass functions using only the three radial bins of the -body models (referred to "Ring" regions).
By comparing the mass functions obtained from the -body models and from the observed data, we can investigate the effects of primordial binaries and black holes on the mass function of Pal 5.
We conducted an analysis to assess the impact of unresolved binaries on the determination of the mass function in the Kuepper field, using the Bin-BH model.The results are depicted in Figure 18.We considered two scenarios for the treatment of binaries in the mass function: • RB (Resolved Binaries): All binaries are resolved, meaning that individual masses of binary components are counted in the mass function.
• URB (Unresolved Binaries):  obs is utilized for mass estimation.This scenario represents a real observation where binaries are unresolved.
The mass functions obtained from the RB and URB scenarios display steeper slopes compared to the observational mass function.
In Figure 19, we present a comparison between the mass functions obtained from the -body models using the URB method and the observational data.The upper panel of Figure 19 shows the number counts ().The -body models exhibit a comparable number of stars within the three Field regions when compared to the observed data.The Bin-noBH model shows a slightly higher number of stars, indicating that a longer evolution time of more than 12 Gyr might be necessary for a better match.However, this slight discrepancy does not impact our comparison with the observed normalized counts.
The median and lower panels of Figure 19 display the normalized cumulative distributions,  f (), for the Field regions and  a () for the Ring regions.In the inner radial bin, no significant difference is observed when comparing  f () and  a ().However, for the median and outer radial bins, a noticeable stochastic scatter is present in  f ().This scatter is particularly evident in the  f () of the Bin-BH model in the outer radial bin.These findings suggest that the observational data may also exhibit similar scatter, and it is important to consider this when comparing the -body model with the observational data.
The standard way to characterize a mass function is by using a power-law form given by the equation: where  is a normalisation constant and  is the power-law index used for fitting.We employ the fitting method outlined in Khalaj & Baumgardt (2013) to determine the statistical error accurately.The formula for fitting  is: where  represents the total number of stars,   is the mass of an individual star,  min is the minimum mass of stars, and  is the ratio of the maximum to the minimum masses of stars.Iterative calculations are necessary to solve this fitting equation.The corresponding error can be described as: The power-law indices of the mass functions () obtained from fitting are summarized in Table 5.In the inner radial bin, the  values for the three Bin models are in rough agreement with the observational data, while the noBin-BH model shows a significantly higher .This result remains consistent when comparing the mass functions within the Field and the Ring regions.
In the middle and outer radial bins, all of the -body models exhibit higher  values compared to the observational data.This discrepancy is more pronounced when considering the normalized cumulative distribution in the Ring regions ( a ()).These differences suggest that the -body models exhibit more pronounced mass segregation than what is indicated by the observational data, although we need to take into account the potential stochastic scatter inherent in the observational data.The presence of BHs does not appear to have a clear impact on the mass functions.The models incorporating primordial binaries exhibit better agreement with the observed data, particularly in the inner radial bin.

Uncertainty of initial condition
Due to the computational expense, we are unable to explore the entire parameter space of the initial condition of Pal 5, resulting in several aspects not being addressed in this study.These include assumptions regarding the properties of primordial binaries, the evolution of the Galaxy, the uncertainty associated with stellar evolution, the gravitational wave kicks following mergers of binary black holes (BBHs), and the realistic formation environment of the cluster.
In our study, we have adopted two extreme assumptions for the primordial binaries (Kroupa and FlatLog) with a 100% initial binary fraction.However, these assumptions may not accurately reflect the true properties of primordial binaries in Pal 5. Nonetheless, Fig, 9 suggests that the initial period distribution has no significant impact on the survival fraction of binaries as a function of period, as long as the cluster possesses a similar initial density profile and orbit in the Galaxy.Furthermore, the evolution of the binary fraction ( hb ) can be utilized to derive the period evolution for different assumptions regarding the initial binary populations.By using a 100% initial binary fraction, we also explore the maximum potential dynamical impact of primordial binaries.The wide range of periods considered allows us to investigate the behavior of hard and soft binaries with and without black holes (BHs).
Our model assumes a static Galactic environment, which is consistent with the setup employed in G21 to facilitate proper comparison.Incorporating a realistic time-dependent Galactic potential, which may be important to understand the density profile of the stream (Pearson, Price-Whelan & Johnston 2017), is challenging due to the limited observational constraints on Galactic evolution.It is plausible that Pal 5 was formed in a significantly different Galactic environment, potentially leading to variations in mass loss and density evolution compared to our models.However, we believe that the overall trend driven by the presence of BHs should be similar.Thus, our results offer a general perspective on how the existence of BHs impacts the binary populations.
The retention of BHs in clusters after supernovae remains an open question based on stellar evolution models.Our models do not consider gravitational wave kicks following BBH mergers, which could lead to an overprediction of massive BBHs with masses exceeding 100  ⊙ .Although such BBHs can influence the timescale of cluster disruption as shown in Figure 4 and 6, their impact on the period distribution of binaries is limited since the hard-soft boundary is not determined by a single specific BBH.
The initial conditions of the clusters assume spherically symmetric Plummer models, similar to previous N-body simulations of GCs.However, the initial complexity of GC formation, including irregular cluster structures prior to achieving virial equilibrium and the presence of gas, may affect the binary populations during the gasembedded phase.

Observation of binaries
In Section 3.2.4,we conducted an analysis to assess the feasibility of detecting binaries by measuring the radial velocity difference (|Δ LOS |) through multiple half-year observations.The maximum time-interval reaches half a year.The results indicate that approximately 40 binaries could be identified, covering a period distribution ranging from a few to 10 4 days.The model without BHs tends to exhibit a higher fraction of long-period binaries.While this observation cannot directly constrain the existence of BHs, it can provide insights into the presence of wide (long-period) binaries.Such information may be valuable in constraining the initial period distribution by utilizing the  hb values depicted in Figure 9.
To obtain a stronger constraint on the existence of BHs, it is crucial to obtain additional observations of binaries in the period range around 10 5 days, which has proven to be challenging thus far.Furthermore, it is necessary to observe binaries in different regions of Pal 5, including the inner region and the distant tail.Given the uncertainties associated with the properties of primordial binaries, assuming an initial period distribution becomes essential for constraining the density evolution based on the observed period distribution of present-day binaries.Notably, wide binaries disrupted within the dense cluster can survive along the low-density tidal tail.Therefore, the difference in the fraction of wide binaries inside the cluster and in the distant tail can help constrain both the initial period distribution of binaries and the density evolution of clusters, ultimately shedding light on the existence of BHs.
Another approach to constrain the BH population is by detecting BH-star binaries.We find four BH-MS binaries with relatively high MS masses, as shown in Table 3 and 4 and also illustrated in Figure 8. Figure 15 suggests that the CSST has the potential to detect CVs, thereby providing additional constraints on binaries with WDs.
Multi-epoch spectroscopic observations for |Δ LOS | offer another possibility to detect non-interacting BH-star binaries.By utilizing this data, we can obtain better constraints on  LOS , providing an indirect constraint on the dynamical impact from BHs in the cluster center.

CONCLUSIONS
In this study, we performed -body simulations of the Galactic halo globular cluster Pal 5 with and without the inclusion of BHs, while considering a significant fraction of primordial binaries.Our main objectives were to investigate the influence of binaries and BHs on the cluster's dynamical evolution and to understand how the presence of BHs affects the binary populations within Pal 5. Additionally, we aimed to determine whether the observations of binary populations could provide indirect evidence for the existence of BHs in Pal 5.
Our findings indicate that the presence of primordial binaries has a noticeable but not drastic effect on the cluster's dynamical evolution, consistent with previous work Wang et al. (2022).In models with BHs, the existence of primordial binaries alters the half-mass relaxation time ( rh ) and reduces the number of BBHs that contribute to binary heating.However, the influence on mass loss and radial evolution is more complex.Models with primordial binaries (Bin-BH and Bin-BH-Alt) exhibit shorter initial  rh compared to models without primordial binaries (noBin-BH model).After 1 Gyr, the situation reverses due to larger half-mass radius ( h ) and lower total BH mass ( BH ) in the Bin models.This trend changes again after 8 Gyr when a massive BBH forms in Bin-BH-Alt, accelerating the cluster's dissolution (see Figure 6).Thus, the tidal dissolution time does not exhibit a simple dependence on the presence of primordial binaries.
In models without BHs and a low initial density (Bin-noBH and Bin-noBH-F), the evolution is more sensitive to the presence of primordial binaries compared to the BH models.Achieving a similar cluster at 11.5 Gyr requires a higher initial density in these cases.
Conversely, the assumption of BH existence significantly affects the population of wide binaries.Over long-term evolution, hard binaries are less affected by dynamical disruption.The fraction of hard binaries remains independent of the initial period distribution (Figure 9).The remaining fraction of wide binaries depends on the evolution of the hard-soft boundary.The period distribution of models with BHs peaks at a shorter period compared to models without BHs, consistent with the hard-soft boundary.However, we find that not all wide binaries outside the hard-soft boundary are immediately disrupted.Many wide binaries outside this boundary can persist in the cluster for a long time.This suggests that the observation of wide binaries may not readily constrain the actual hard-soft boundary and be used to determine the cluster's density evolution history.
We have found that multi-epoch spectroscopic observations can detect most binaries with bright stars and periods below 10 4 days.By excluding these binaries, the measurement of  LOS of bright stars can be significantly improved, providing better indirect constraints on the BH population through dynamical analysis.
Additionally, we have identified 4 BH-MS binaries in the Bin-BH model at 11.5 Gyr, which could potentially be detected using the same method, offering an additional possibility to provide evidence for the existence of BHs.
We also investigated how binaries and BHs influence the presentday mass function of Pal 5. Our results suggest that models with primordial binaries have mass function more consistent with the observational data, while the impact of BHs on the mass function is weak.All -body models exhibit mass segregation features that are not observed in the outer region of Pal 5.However, it is important to consider the potential impact of stochastic scatter, which may in-fluence the conclusions drawn from the comparison.This indicates the need for alternative initial mass functions or additional observations of mass functions, with improved statistical precision, to better understand the underlying reasons for this discrepancy.

Figure 1 .
Figure 1.The orbit of the Pal 5 in the Galactocentric frame.The upper and lower panels show the projected trajectory in the  G −  G and  G −  G planes, respectively. G is the projected radial coordinate in the  G −  G plane.The symbols '+' and 'x' represent the zero-age and present-day positions, respectively.

Figure 2 .
Figure 2. Initial periods () v.s.eccentricities () of primordial binaries for the Kroupa binary model and the FlatLog model.The central plot of each panel shows - of individual binaries.The upper and the right histograms show the normalized distribution of  and , respectively.The distribution of massive binaries is shown by blue lines.

5 10Figure 3 .
Figure 3.The evolution of two-component half-mass relaxation time for all models ( rh ; upper two panels) and  factors (lower panel) for BH models.

Figure 4 .
Figure 4.The evolution of half-mass radius of all objects ( h ; dashed curves) and the half-mass radius of BHs ( h,BH ; solid curves).The wBH-1 and noBH-1 models from Gieles et al. (2021) are shown as references.

Figure 5 .
Figure 5.The evolution of the bound mass (), the BH mass fraction (  BH ) and the bound mass of BHs ( BH ) for all models.The data of the wBH-1 and noBH-1 models are shown for comparison.

Figure 6 .
Figure 6.The evolution of the semi-major axes of BBHs within the core radius ( c ) of the three BH models.The colors of the lines indicate the masses of the BBHs.We can observe a reduction of the semi-major axes of individual BBHs, indicating their dynamical hardening over time ( > 10 AU) and inspiral by GW radiation ( < 1 AU).

Figure 7 .
Figure 7.The surface number density (Σ ()) profiles are presented for the  -body models along with observational data from Ibata et al. (2017).The upper panel displays snapshots of the  -body models at the present-day Galactic position and at apporximately 11.5 Gyr.The lower panel shows  -body snapshots that match the observed Σ () profile.The ages of the corresponding snapshots ( mat ) are indicated in the legend.Vertical lines are used to indicate the the 'effective radius' -the radius containing half the number of stars in projection -( hn ) of the clusters.

Figure 8 .
Figure 8.The contour of - b at 11.5 Gyr for the Bin-BH model (upper panel) and the Bin-noBH model (lower panel).Binaries with one or two compact objects are excluded in the contour.Instead, BH-MS and BH-WD binaries are marked as blue and lightblue stars, respectively.Three curves show the hard-soft boundaries  hs ( ) at zero age, 100 Myr and 11.5 Gyr, respectively.The white region outside the color region indicates no binary.

Figure 9 .
Figure9.The period distribution of binaries within  h at two different stages: the initial phase (represented by steps) and at 5 Gyr (shown as filled histograms).The upper panel displays the number of binaries within  h normalized by the bound mass of the cluster ( hb ) .The lower panel shows the ratio between  hb at 5 Gyr and the initial  hb .The vertical dashed and solid lines represent the hard-soft boundary of period within  h at 0 and 5 Gyr, respectively.

Figure 10 .
Figure10.Upper panel: binary fractions of all objects along the 3D radial direction for the Bin-BH and Bin-noBH models; Lower panel: prediction for the observed binary fractions with an I-band magnitude range of 20.5 and 23 mag (corresponding to main sequence stars) and mass ratio > 0.6.

Figure 11 .
Figure 11.The line-of-sight velocity difference of binaries (|Δ LOS |) as a function of period for multi epochs of observation.The initial snapshots of the three models are chosen at  =  mat .Each binary type, classified according to the sse (Single Stellar Evolution) code, is represented by a different color.The stellar types include: MS (Main Sequence), HG (Hertzsprung Gap), GB (First Giant Branch), CHeB (Core Helium Burning), AGB (Asymptotic Giant Branch), and WD (White Dwarf).

Figure 12 .
Figure12.The number counts of bright binaries with post-main-sequence component for three models at .The legend "tot" include all binaries and the "obs" include only detectable binaries with |Δ LOS | > 0.3 km/s.

Figure 13 .
Figure 13.The line-of-sight velocities of individual bright stars and binaries are plotted, and detectable binaries with |Δ LOS | > 0.3 km/s are indicated as green dots.

Figure 14 .
Figure 14.Illustration of the BH-MS formation process.The black and grey circles represent BHs, and the blue circles represent MS stars.

Figure 15 .
Figure 15.The color-magnitude diagram of the Bin-BH model at 11.5 Gyr.Red points are single stars.Other points are unresolved binaries where colors represent mass ratio ().The black crosses are BH-MS binaries shown in Table3.The left panel corresponds to the HST F555W-F814W and F555W filters, while the right two panels correspond to the CSST g-i and g, and u-y and u filters, respectively.

Figure 16 .
Figure 16.The total masses ( tot ) v.s. the F555W-band flux-converted masses ( obs ) for main-sequence binaries of the Bin-BH model at 12 Gyr.The grey line shows the case of  tot =  obs .Colors represent mass ratio ().

Figure 17 .Figure 18 .
Figure 17.The 2-dimensional density map of the noBin-BH model at 11.8Gyr.The color contours with solid lines represent the Smith and Kuepper fields, which have available HST data.The boundaries of the three ring radial bins are indicated by dashed grey circles.Two approaches are employed for selecting samples to measure the mass functions: 1) using the intersection between the Smith/Kuepper fields and the ring regions (referred to as "Field" regions); and 2) using only the ring regions themselves (referred to as "Ring" regions) to enhance statistical accuracy.

Figure 19 .
Figure 19.The mass functions of four  -body models in three radial bins, with the observational data shown as a reference.The upper panel displays the number counts (), the middle panel shows the normalized cumulative distribution  f () for the Field regions, and the lower panel shows the normalized cumulative distribution  a () for the Ring regions.

Table 2 .
The table displays the line-of-sight velocity dispersion ( LOS ) estimated from bright stars and binaries.The last column,  1D , represents the estimation of  LOS based on Equation 6, which serves as the unit for the other four columns.In particular, the column  LOS,S,hn presents the  LOS value derived from single stars within  < 3 arcmin (17 pc, approximately the  hn ).The remaining three columns depict  LOS within  < 10 arcmin (58 pc), where  LOS,S ,  LOS,SB , and  LOS,SCB represent the  LOS values from only single stars, both single stars and binaries, and both single stars and undetectable binaries with |Δ LOS | ≤ 0.3 km/s, respectively.LOS .In our -body model, we simulate the impact of excluding binaries with |Δ LOS | > 0.3 km/s on the determination of  LOS .Figure13displays the individual line-of-sight velocities of bright stars ( LOS ), undetectable bright binaries with |Δ LOS | ≤ 0.3 km/s ( LOS,SB ), and detectable binaries with |Δ LOS | > 0.3 km/s, aligned with the projected distance.Most binaries with  LOS > 1 km/s are detectable, and thus, we can remove them for the calculation of  LOS .Table2demonstrates how removing detectable binaries improves the determination of  LOS .To have a consistent comparison among the three models, we scale the value of  LOS by the estimated 1-

Table 4 .
1 [  ⊙ ]  2 [  ⊙ ]The accretion events of BwBHs after 100 Myr.The "Primordial" column indicates whether the binary is primordial (formed during the initial star cluster formation) or dynamically formed (formed through interactions within the star cluster after its formation).The "Type" column indicates the combination of binary companions.The secondary stellar types involved in the accretion events include: MS, HG , GB, CHeB, AGB, HeHG (Hertzsprung Gap Naked Helium star), WD and NS (Neutron star).

Table 5 .
Fitting result of the power-law indices () of the mass functions in different radial bins.The column labeled "region" distinguishes between the Smith and Kuepper fields (referred to "Field") and the ring regions (referred to "Ring").