Compact Object Modeling in the Globular Cluster 47 Tucanae

The globular cluster 47~Tucanae (47~Tuc) is one of the most massive star clusters in the Milky Way and is exceptionally rich in exotic stellar populations. For several decades it has been a favorite target of observers, and yet it is computationally very challenging to model because of its large number of stars ($N\gtrsim 10^6$) and high density. Here we present detailed and self-consistent 47~Tuc models computed with the \texttt{Cluster Monte Carlo} code (\texttt{CMC}). The models include all relevant dynamical interactions coupled to stellar and binary evolution, and reproduce various observations, including the surface brightness and velocity dispersion profiles, pulsar accelerations, and numbers of compact objects. We show that the present properties of 47~Tuc are best reproduced by adopting an initial stellar mass function that is both bottom-heavy and top-light relative to standard assumptions \citep[as in, e.g.,][]{Kroupa2001}, and an initial Elson profile \citep{Elson1987} that is overfilling the cluster's tidal radius. We include new prescriptions in \texttt{CMC} for the formation of binaries through giant star collisions and tidal captures, and we show that these mechanisms play a crucial role in the formation of neutron star binaries and millisecond pulsars in 47~Tuc; our best-fit model contains $\sim 50$ millisecond pulsars, $70\%$ of which are formed through giant collisions and tidal captures. Our models also suggest that 47~Tuc presently contains up to $\sim 200$ stellar-mass black holes, $\sim 5$ binary black holes, $\sim 15$ low-mass X-ray binaries, and $\sim 300$ cataclysmic variables.


INTRODUCTION
As laboratories for studying both small-scale stellar physics and gravitational dynamics as well as largerscale galaxy formation and evolution, globular clusters (GCs) have been a main focus of computational N -body modeling for several decades. They are efficient in producing exotic objects such as X-ray binaries (e.g., Heinke et al. 2005;Clark 1975), millisecond pulsars (MSPs; e.g., Ransom 2008), and cataclysmic variables (CVs; e.g., Knigge 2012) given their high stellar densities. GCs may also account for a significant fraction of the overall binary black hole (BH) merger rate in the local universe (Abbott et al. 2021;Kremer et al. 2020; Rodriguez launch of HST, Guhathakurta et al. (1992) exploited its ability to resolve individual stars in crowded fields, constructed an accurate color-magnitude diagram (CMD), and placed an upper limit on the mass of compact objects at the cluster center. Subsequent HST campaigns aimed at obtaining detailed photometry of the cluster core (de Marchi et al. 1993) and deriving stellar luminosity and mass functions (Santiago et al. 1996). CV candidates in the cluster were found through many epochs of observation, long exposure times and almost continuous observing (Edmonds et al. 1996;Shara et al. 1996;Gilliland et al. 1995). Combining X-ray observations from Chandra and 120 orbits of HST, Edmonds et al. (2003a,b) observed the largest number of optical counterparts to the X-ray sources in any GC, and detected many CVs, active binaries, and constrained the total number of MSPs. More recently, to further explore the dynamical structure of 47 Tuc, several studies obtained proper motions and line-of-sight velocities with HST (e.g., McLaughlin et al. 2006;Watkins et al. 2015;Heyl et al. 2017), Gaia (Baumgardt et al. 2019;, MUSE (Kamann et al. 2018), and archival/literature data from a combination of instruments (Baumgardt 2017;Baumgardt & Hilker 2018) to construct velocity distributions and examined the rotation and velocity anisotropy of the cluster (e.g., Bellini et al. 2017).
Extensive previous studies have modeled 47 Tuc using static models, including single-mass King models (e.g., Illingworth & Illingworth 1976), multi-mass generalizations of King models (e.g., Da Costa & Freeman 1985), Michie-King models (e.g., Meylan 1988Meylan , 1989, models with added rotations (e.g., Davoust 1986;Bellini et al. 2017), and models incorporating the effects of the Galactic tidal field on the structure of the cluster (e.g., Hénault-Brunet et al. 2019). However, there have been few previous attempts at modeling the cluster using the N -body method, despite this wealth of observational data on 47 Tuc. This is not surprising since modeling GCs as massive and dense as 47 Tuc requires significant computational time and sophisticated N -body codes. It takes more than a year for the current best direct Nbody code, NBODY6++GPU (Wang et al. 2015), to simulate a massive GC similar to 47 Tuc with ∼ 10 6 stars. Therefore, the more rapid Monte Carlo technique is necessary to explore the parameter space of initial conditions for 47 Tuc (see Section 2). The only two previous studies of 47 Tuc using N -body modeling either rely upon scaling-up of low mass models (Baumgardt & Hilker 2018), or apply artificially large natal kicks for BHs (Giersz & Heggie 2011) so there were only a few BHs in their models. Models without a realistic number of BHs may still show close matches with observations such as the surface brightness profile (Giersz & Heggie 2011), but would potentially miss the BH dynamics that also significantly affects other exotic objects through gravitational encounters. To uncover the cluster's dynamical properties and understand all evolutionary pathways that shape the observable properties of 47 Tuc and other similar clusters, we develop accurate and self-consistent models of 47 Tuc that will match different observed cluster features simultaneously. Significantly, these models include most up-to-date prescriptions for compact object formation, which in particular enables us to examine pulsar and BH dynamics in 47 Tuc for the first time.
The CMC Cluster Catalog models in Kremer et al. (2020) are representative of most typical present-day Milky Way GCs, so we first searched for a best-fit 47 Tuc model in the catalog using the χ 2 statistic method . We found that while some catalog models can match the inner parts of the observed surface brightness and velocity dispersion profiles, there are large discrepancies at the outer parts of the good-fit models (Figure 10). This illustrates the necessity to simulate custom models for 47 Tuc with smaller core radii (less compact objects, especially BHs) and more low-mass stars to puff up the outer part of the cluster. This paper is organized as follows. In Section 2, we summarize the main computational methods used to model 47 Tuc, and describe the implementation of tidal capture and collision with giant stars to our cluster simulations. We demonstrate that our 47 Tuc models fit the observations by comparing to multiple observational properties such as the surface brightness profile and velocity dispersion profile in Section 3. We show the compact object populations in a best-fit model in Section 4. In Section 5, we discuss models constructed in search of the initial conditions that best agree with the observational data, caveats and uncertainties in our models. Finally, in Section 6, we summarize our results.

METHODS
We simulate 47 Tuc using the Cluster Monte Carlo code (CMC; e.g., Rodriguez et al. 2021b, and references therein), a public, fully parallelized Monte Carlo code based on classic work by Hénon (1971a,b). It is capable of modeling the complete dynamical evolution of realistic dense star clusters containing ∼ 10 6 stars over a Hubble time in just a few days. It allows us to compute self-consistent models of even the largest star clusters, covering the full range of masses all the way to the brightest systems. CMC incorporates all relevant physics, including two-body relaxation, tidal mass loss, and strong dynamical interactions of single stars and binaries (e.g., exchange interactions, tidal capture, direct collisions). Stellar evolution is fully coupled to the stellar dynamics and is computed with the publicly available COSMIC software (Breivik et al. 2020) based on Hurley et al. (2000Hurley et al. ( , 2002. The Fewbody package is used to directly integrate all three-and four-body gravitational encounters (Fregeau et al. 2004;Fregeau & Rasio 2007), including post Newtonian dynamics for BHs (Antognini et al. 2014;Amaro-Seoane & Chen 2016;Rodriguez et al. 2018;Rodriguez et al. 2018b). In the following subsections, we describe simulation features adopted here that differ from our past CMC assumptions.

Tidal Capture
Following the early works of Fabian et al. (1975), Press & Teukolsky (1977), and Lee & Ostriker (1986), we implemented tidal capture of MS stars into CMC. When a compact object and a MS star or two MS stars fly close during single-single encounters, tidal force will perturb the stars and excite non-radial oscillation within them. The oscillation energy comes from the two objects' kinetic energy, and the stars are captured if the oscillation energy exceeds that of their initial kinetic energy. The amount of energy deposited in a star from one passage of encounter can be written as where M * and R * are the mass and the radius of the star, and M is the mass of a compact object or another MS star. R p is the distance of closest approach. T l is a dimensionless function measuring energy contribution from different harmonic modes l. η * is a measure of the duration of periastron passage, and is written as We adopt the fitting formulae obtained from Lee & Ostriker (1986) in Kim & Lee (1999) to calculate the cross sections and pericenter distances R p . The dimensionless function T l and the oscillation energy in the MS stars are estimated using the formulae in Portegies Zwart & Meinen (1993), which were calculated by fitted to the hydrodynamical studies (Lee & Ostriker 1986;Ray et al. 1987). We include only the contribution of the quadrupole and octupole terms in the tides, while higher-order corrections do not contribute significantly to the oscillation energy.
Since these prescriptions assumed polytropic stellar models, we apply polytropic index n = 1.5 for low-mass MS stars with M < 0.7 M and naked helium MS stars, which correspond to star type 0 and 7 in COSMIC (Breivik et al. 2020), respectively. We use n = 3 for other MS stars with M > 0.7 M , corresponding to star type 1 (McMillan et al. 1987;Kim & Lee 1999). For a comparison between the simple polytropic stellar models with detailed stellar models, see Appendix B.
The initial total energy of the two stars at infinity can be written as E orb = 1 2 µv 2 ∞ , where µ is the reduced mass and v ∞ is the relative velocity at infinity. We assume that if the total oscillation energy deposited in the stars during the first passage is larger than the initial total energy, then the two stars are bound and form a tidal capture binary. The binary semi-major axis immediately after tidal capture is set to be 2R p , assuming angular momentum conservation during tidal capture and immediate circularization of the binary after tidal capture. We use this simple treatment of final binary separations to maximize the number of binaries formed from tidal capture (Lee & Ostriker 1986), especially NS-MS star binaries that could become redback MSPs in clusters (see also Section 4.1). In the special cases when Roche lobe overflow starts at 2R p , we merge the two interacting objects instead of forming a binary.

Giant Star Collision
We implemented special treatments for single-single collisions involving giant stars into CMC. Here we define "collision" as when the two interacting objects pass within a small factor of the sum of their radii or the radius of a giant star. When giant star collision is activated in CMC, we assume that a giant star collides with another star or compact object when R p ≤ 1.3R g following Lombardi et al. (2006), where R p is the pericenter distance between the two objects and R g is the radius of the giant star. For the collision of two giants, we set the collision pericenter distance limit to R p ≤ 1.3(R g1 + R g2 ), where R g1 and R g2 are the radii of the giants. At the same time, the "sticky sphere" direct collision prescription (R p < R 1 + R 2 , where R 1 and R 2 are the radii) is applied to all other stars, where the mass of the collision product is equal to the total mass of the two colliding objects (Fregeau & Rasio 2007).
These collisions with giant stars can lead to the formation of tight binaries. Qualitatively, we expect a collision involving at least one giant star to proceed similar to the common envelope process where the cores orbit inside an extended envelope. Through drag forces, the binary inspirals and deposits energy into the envelope until eventually the envelope is expelled, leaving behind a compact binary. We assume that there is no mass transfer to a companion star from the envelopes of the giant stars during these collisions. For single-single collisions involving only one giant, we assume that a binary star formed consists of the core of the giant star and the interacting main-sequence star (MS star) or compact object.
To calculate the binary semi-major axis, we adopt equations (3)-(5) below -equations (4)-(6) from Ivanova et al. (2006) -based on smoothed particle hydrodynamics simulations (Lombardi et al. 2006). The semi-major axis is chosen as the minimum of a coll and a ce from these equations.
a coll is calculated as where a ce can be obtained by solving the equation with e ce = 0. M g , M , and M c are the giant's mass, the mass of the non-giant object, and the core mass of the giant, respectively. α ce = 1 and λ = 0.5 are parameters related to common-envelope evolution (e.g., Hurley et al. 2002; see also discussions in Appendix C). For collisions between two giant stars, we instead calculate the binary semi-major axis using the expression with e ce = 0. M g1 , M c1 , M g2 , and M c2 are the masses and the core masses of the first and second giants, respectively. The final binary is composed of the cores of the two giants and both giant envelopes are ejected. In special cases where the MS star or a giant core starts mass transfer after the binary semi-major axis is computed, we assume the collision does not lead to binary formation. Instead we merge the two original interacting objects to avoid Roche Lobe overflow during a common-envelope phase, the outcome of which is highly uncertain. The stellar types of the merger products follow those in Table 2 of Hurley et al. (2002).

47 Tuc Simulations
For 47 Tuc simulations, we fix the metallicity Z = 0.0038 (Harris 1996(Harris , 2010, the Galactocentric distance r gc = 7.4 kpc (Harris 2010;Baumgardt et al. 2019) and the initial binary fraction f b = 0.022 (Giersz & Heggie 2011). Binary companion masses are drawn from a flat distribution in mass ratio q, in the range 0.1 − 1 (e.g., Duquennoy & Mayor 1991). Binary orbital periods are drawn from a distribution flat in log scale from Roche Lobe overflow to the hard/soft boundary, and binary eccentricities are drawn from a thermal distribution (e.g., Heggie 1975). To calculate the tidal radius of a cluster, we assume a circular orbit with radius r gc around the Galactic center and a flat rotation curve of the galaxy with circular velocity v gc = 220 km s −1 . The tidal radius of a cluster is given by where M clu is the total mass of the cluster (Spitzer 1987;Chatterjee et al. 2010). This calculation of tidal radius is valid for 47 Tuc since the cluster's orbit in the galaxy has a small eccentricity (< 0.2) and is almost circular . Following Giersz & Heggie (2011), we adopt a twopart power-law initial mass function (IMF) with masses in the range of 0.08 − 150 M , a break mass at 0.8 M and power-law slopes of α 1 = 0.4 and α 2 = 2.8 for the lower and higher mass part, respectively. Neutron stars (NSs) formed from core-collapse supernovae receive natal kicks drawn from a Maxwellian distribution with a standard velocity deviation σ NS = 265 km s −1 (Hobbs et al. 2005). The natal kicks for NSs from electroncapture supernovae and accretion-induced collapse are also drawn from a Maxwellian distribution but with a smaller standard deviation σ ecsn = 20 km s −1 (Kiel et al. 2008;Ye et al. 2019). We apply mass fallback kicks for BHs and sample their natal kicks from the same Maxwellian distribution for core-collapse NSs but with scaled-down velocities (Fryer et al. 2012). The kicks are reduced in magnitude by a fractional mass of fallback where v NS is the velocity of NS drawn from Hobbs et al. (2005) for core-collapse supernovae and f f b is the fractional parameter of the stellar envelope that falls back upon core collapse (Belczynski et al. 2002;Fryer et al. 2012;Morscher et al. 2015). Compared to Giersz & Heggie (2011), our simulations treat BHs and NSs more realistically by applying fallback BH natal kicks and allowing for NSs to form in electron-capture supernovae and accretion-induced collapses with low kicks, which are essential for explaining the large number of pulsars in GCs (e.g., Pfahl et al. 2002;Podsiadlowski et al. 2004;Ye et al. 2019).
We assume the 47 Tuc models are initially described by Elson density profiles (Elson et al. 1987), where γ is a free parameter of the power-law slope (γ = 4 gives a Plummer sphere; Plummer 1911), and a is a scaling parameter. Unlike previous CMC papers that varied the initial cluster parameter space in a grid format (e.g., Kremer et al. 2020;Weatherford et al. 2021), here we have run 15 sim-ulations (see Table 3 in Appendix D) exploring a specific parameter space for 47 Tuc in the initial number of stars, density profile, binary fraction, virial radius, tidal radius and IMF to search for a best-fit model where the model cluster's late time (9 − 13.8 Gyr) surface brightness profile (SBP) and velocity dispersion profile (VDP) first match the observations. We also vary γ of the Elson profile to analytically fit the outer slope of the observed surface brightness profile. All simulations are evolved to 13.8 Gyr. Table 1 shows the initial conditions of our best-fit model for 47 Tuc. Note-From left to right: initial number of stars, Elson profile γ index, initial virial radius, metallicity, Galactocentric distance, initial tidal radius, minimum mass, break mass, and maximum mass in the IMF, respectively, power-law slopes for the lower mass and higher mass part of the IMF, respectively, and initial binary fraction.

MATCHING OBSERVED CLUSTER PROPERTIES
We found best-fit models by comparing the χ 2 likelihood between different snapshots (time steps) of the simulation and selecting the ones with small χ 2 as in Rui et al. (2021). The model SBP is calculated by making two-dimensional projections at each time step of a simulation, assuming spherical symmetry (also see Kremer et al. 2018b). We exclude very bright stars with luminosity L > 12 L to avoid large fluctuations in the SBP, and also low-mass, faint stars with M ≤ 0.85 M . Changing the cutoff luminosity from 12 L to 15 L has negligible effect on the SBP fit. We use the method in Pryor & Meylan (1993) to find the model VDP, including in the calculation all giant stars and upper MS stars with a cutoff mass M > 0.85 M (the main-sequence turnoff mass of the model is ∼ 0.9 M ). Small variations in the cutoff mass for the VDP calculations do not affect the fit. For more details on the SBP and VDP calculations, see Kremer et al. (2018b) and Rui et al. (2021), and the references therein. Similar to the VDP, we use the same selection of stars to calculate the projected model number density profile (NDP). Note that, unlike the VDP, the model NDP is very sensitive to changes in the cutoff mass. The cluster properties only vary negligibly between ∼ 9 − 12 Gyr, and Figure 1 shows the snapshot at 10.55 Gyr as a representative of the best-fit 47 Tuc models. At this time, the best-fit model is able to closely match the observed SBP, VDP, NDP, and pulsar acceleration distribution simultaneously.
The initial tidal radius of the cluster with a circular orbit at r gc = 7.4 kpc is 182 pc. The best-fit cluster starts out to r max = 800 pc initially, and is overfilling its tidal radius (r max is a free parameter in our simulations). Observationally, young star clusters (potential GC progenitors) appear to not be tidally truncated, and may be surrounded by unbound stars at the outer region that have not yet been stripped away by the Galactic tidal force (e.g., Elson et al. 1987Elson et al. , 1989Mackey & Gilmore 2003). Using r max = 300 pc has negligible effect on the final cluster properties while using r max = 200 pc only changes the final properties slightly. Furthermore, it is shown in the top panel of Figure 1 that the model SBP deviates from the observed SBP at the edge of the cluster where r 1500 arcsec. Attempting to decrease this discrepancy, we tried a simulation with r gc = 5.5 kpc (pericenter of the observed 47 Tuc orbit instead of the apocenter at 7.4 kpc; Baumgardt et al. 2019) and initial tidal radius of 149 pc. Although this did lower the outermost portion of the model SBP, the discrepancy was still present. This could be because our models do not take into account rotation (see Section 5.5) or because the observations at the edge of the cluster are incomplete due to contamination from field stars (the fit is better for the NDPs in the bottom panel where the observations at the edge of the cluster are taken with Gaia; de Boer et al. 2019).
The age of 47 Tuc can be estimated by modeling the properties of the cluster's white dwarf (WD) cooling sequence, comparing its CMD to theoretical isochrones, or fitting stellar evolution tracks to observed eclipsing bi- naries in the cluster. The estimated ages are uncertain and span a wide range from ∼ 9 − 14 Gyr (e.g., Dotter et al. 2010;Hansen et al. 2013;VandenBerg et al. 2013;Brogaard et al. 2017;Thompson et al. 2010Thompson et al. , 2020, and references therein). Our models predict that the age of 47 Tuc is between ∼ 9 Gyr to ∼ 12 Gyr, consistent with these estimates. Figure 2 compares the best-fit model's velocity anisotropy profile to the observations. Some recent studies (e.g., ) have shown that some GCs exhibit velocity anistropy, which may be caused by rotation. In particular, observations of 47 Tuc (Watkins et al. 2015;Heyl et al. 2017; using both HST and Gaia show small velocity anisotropy both in the inner and outer parts of the cluster. Our best-fit model is consistent with being isotropic and, overall, is within the velocity anisotropy uncertainties of the observations.  (Watkins et al. 2015;Heyl et al. 2017; show slight tangential velocity anisotropy in the inner part and radial velocity anisotropy in the outer part. The model is isotropic in general but still within the uncertainties of the observational points. Since 47 Tuc contains the second largest population of MSPs (28 observed to date; Paulo Freire's pulsar catalog 2021), 47 Tuc models must be able to explain its pulsar accelerations in the cluster potential. Accurate pulsar line-of-sight accelerations are obtainable from timing analysis of pulsar binary orbital periods. Upper limits on the acceleration are also calculable from the pulsar spin derivatives, in the case of an isolated pulsar or binary with incomplete timing (Ridolfi et al. 2016;Freire et al. 2017;Freire & Ridolfi 2018;Ridolfi et al. 2021).
The observed line-of-sight accelerations and upper limits for 47 Tuc MSPs are shown in Figure 3 (left panel, orange and black markers 1 ).  . The theoretical core radius is mass-density-weighted and calculated using the method in Casertano & Hut (1985). The model observational half-light radius is defined to be the 2D-projected radius that contains half of the cluster's total light, and the model observational core radius is estimated using the method described in Morscher et al. (2015); Chatterjee et al. (2017). .
We estimate the maximum line-of-sight accelerations of model pulsars using equation 2.5 in Phinney (1992), where M cyl is the mass of the cluster within a cylindrical tube of projected radius R ⊥ from the cluster center along the line-of-sight. This is shown by the green curves in Figure 3 (left panel). We also include the distribution of maximum possible acceleration allowed by the cluster potential in Figure 3 ( of radius r. To 1 σ uncertainty, all but one of the observed pulsars have line-of-sight accelerations within the bounds allowed by the model potential, demonstrating further that our best-fit model is consistent with the observations. The only outlier at ∼ 4 arcsec is still enclosed by the maximum possible total acceleration allowed by the cluster potential. The right panel of Figure 3 shows the cumulative number of all model NSs (light blue curve), model MSPs (orange curve), and observed 47 Tuc MSPs (black curve) as a function of the projected distance from the cluster center. All of the observed 47 Tuc pulsars are MSPs with spin periods < 30 ms. As expected, the observed MSP distribution agrees better with the modeled MSP distribution than the total NS distribution. The total NS population is distributed slightly further out because MSPs are formed through extended periods of mass accretion and are slightly more massive (∼ 0.1 M more) than most other NSs. Mass segregation leaves the MSPs closer to the cluster center than typical NSs. Also NSs closer to the cluster center are more dynamically active, so they more easily acquire companion stars through encounters and become MSPs . Table 2 lists the properties of the best-fit 47 Tuc models at the present day, including the total numbers of compact objects and blue straggler stars (BSSs). In addition, Table 4 in Appendix D shows the more general mean and standard deviation values of the same properties for best-and near-fit models from simulations 1-3 in Table 3. The inclusion of simulation 3 with a high initial binary fraction leads to large standard deviations for the numbers of some compact objects. Figure 4 shows the cumulative number distributions of all stellar populations (colors denoted in legend) at an early and a late time of the 47 Tuc simulation. As expected, the cluster core is dominated by BHs due to mass segregation, mixing with the most numerous MS stars initially (see also Kremer et al. 2020), and also WDs at late times. This is also reflected in the cluster's SBP, where it is flat within ∼ 20 arcsec, consistent with the features of a non-core-collapsed cluster.

EXOTIC OBJECTS IN 47 TUC
Here we report the number of exotic objects in the best-fit 47 Tuc model, which can be used to estimate the total numbers of compact objects, including BHs, MSPs, LMXBs and CVs, and numbers of BSSs, in 47 Tuc. Most compact object binaries discussed below, including the binary BHs in Section 5.3, are formed through dynamical interactions; only ∼ 10% of MSPs and ∼ 35% of CVs are formed in primordial binaries. This high fraction of dynamical compact object binaries is not surprising given the high density of the 47 Tuc models and the low initial primordial binary fraction we adopted, which is motivated by observations of 47 Tuc (Milone et al. 2012; see also Section 5.4).

Pulsars
A total of 28 pulsars (all of them MSPs) have been observed in 47 Tuc and 18 of them are in binaries (Paulo Freire's pulsar catalog 2021). Among the MSPs in binaries, 6 of them are "black widow" MSPs with low-mass, planet-or white/brown dwarf-like stellar remnant companions (  where the giant star's core becomes a NS in subsequent binary and stellar evolution (see also Section 2.2). A few of the giant star collisions are between a WD and a giant star. Figure 5 shows the distribution of formation times of all MSPs formed in a Hubble time in the best-fit 47 Tuc simulation. At early times of the cluster's evolution ( 3 Gyr), collisions with giant stars play a significant role (∼ 80%) in forming the binaries which host MSPs later, or lead to more binary-mediated dynamical encounters and subsequent MSP formation. This is shown by the peak in the formation time distribution in Figure 5. After ∼ 3 Gyr, more than half of the BHs have been ejected and the cluster is denser than earlier. The higher density leads to enhanced dynamical interactions and the formation of subsequent MSP-host binaries through exchange encounters or tidal capture. In total, ∼ 40% of MSPs formed after 3 Gyr.
In Figure 6, we further compare the orbital periods and companion masses of all model MSP binaries and tidal capture NS-MS star binaries from 9 − 12 Gyr to the observations. Roughly speaking, CMC models can 2 We count the 26 MSPs with very low-mass ( 0.01 M ) CO/ONeMg WD companions, where the WDs are likely unstable as isolated MSPs. This is because there may exist minimum masses for WDs with sufficiently high temperature (T 10 4 K) during mass transfer where there is no equilibrium solution of the WDs' equation of states and the WDs may "evaporate" (Bildsten 2002;Deloye & Bildsten 2003    produce all three groups of pulsar binaries that bracket the ones observed (black widows, redbacks and MSPs with WD companions, respectively). It is important to point out that tidal capture can produce redback-like NS-MS star binaries (orange circles in Figure 6, 5 per snapshot) with similar companion masses and orbital periods, but these 5 NSs are not spun up to MSPs in the models. In a realistic scenario during tidal capture, a MS star is likely to expand from the injection of orbital energy and fill its Roche Lobe (e.g., Kochanek 1992). Thus it is reasonable to assume that a captured, puffedup MS star can transfer mass to a NS and spin it up to a MSP. However, CMC currently lacks detailed treatments of MS stars in tidal capture interactions, and COSMIC does not include the latest pulsar physics such as pulsar irradiation (e.g., Chen et al. 2013;Ginzburg & Quataert 2020, and references therein). We will consider these details and study model MSPs and tidal-captured NS-MS star binaries that overlap with the observed pulsar binaries in future works. In total, the best-fit model produced 54 MSPs including the NS-MS star binaries formed in tidal capture interactions. This number is larger than the actual observed number of MSPs in 47 Tuc likely because pulsar searches are limited by selection effects such as the interstellar dispersion and scattering, and survey thresholds (e.g., Lorimer 2008, and references therein). Our estimate is consistent with the upper-limit estimates of the potential total number of MSPs in 47 Tuc from previous studies (Heinke et al. 2005;Abdo et al. 2009).

Low-mass X-ray Binaries and Cataclysmic Variables
Observations have detected 370 X-ray sources within 2 79 of 47 Tuc using the Chandra Observatory We define LMXBs as BHs or NSs in mass-transferring binaries with WD, giant star, or MS star companions in our models. There are ∼ 14 LMXBs on average per snapshot in the time scale of 3 Gyr in our best-fit models, where ∼ 3 are NS-MS star binaries, ∼ 5 are BH-WD binaries, ∼ 1 is a BH-giant star binary, and ∼ 5 are BH-MS star binaries. There are a similar number of LMXBs at 10.5 Gyr. Similar to LMXBs, we define CVs as masstransferring WD binaries with MS star companions (e.g., Kremer et al. 2021). There are on average ∼ 334 CVs at ∼ 9 − 12 Gyr, and a similar number at 10.5 Gyr. Most of the CVs (∼ 260) and almost all of the LMXBs identified in our models are within the theoretical half-mass radius, consistent with the detected number of X-ray sources. These are only general comparisons between the number of LMXBs and CVs in models and from observations, since we do not consider properties such as duty cycles and luminosities of these mass-transferring binaries. It is not surprising that our models predict more LMXBs and CVs than are observed. The observa-tions of LMXBs and CVs are also limited by selection effects, including the stellar crowdedness in GCs and the faintness of some of these binary systems (Heinke et al. 2003 estimated that there are 7 times more quiescent LMXBs than bright LMXBs in GCs). The number of CVs from the best-fit 47 Tuc models is consistent with the estimates in Ivanova et al. (2006) and Rivera Sandoval et al. (2018), which showed that clusters similar in mass to 47 Tuc can have ∼ 300 CVs. This provides more support for the number estimates of compact objects in our 47 Tuc model. −37 (2σ uncertainties) stellar-mass BHs in 47 Tuc, a few percent of which may be found in mass-transferring binaries similar to the X9 source (e.g., Kremer et al. 2018a). It has also been suggested that 47 Tuc may host an intermediate-mass BH (Kızıltan et al. 2017) with mass of 2300 1500 −850 M . Initially, the model cluster produced ∼ 1270 BHs. Most of them are ejected later through dynamical interactions, as the escape velocities of the cluster decrease significantly over a Hubble time (Figure 7). On average there are ∼ 186 BHs and ∼ 3 binary BHs over a 3 Gyr timescale from ∼ 9 Gyr to ∼ 12 Gyr, while there are 182 BHs and 2 binary BHs at the best-fit model ( Table 2). The number of retained BHs at present is consistent with the 2σ upper limit, and the total mass of the BHs per snapshot (∼ 2350 M ) is consistent with the 1σ upper limit in Weatherford et al. (2020). However, the BH mass in our best-fit models is twice as large as in Hénault-Brunet et al. (2019), indicating different IMFs may vary the final retained BH mass significantly (although they did not specify the maximum and minimum masses used in their IMF). Furthermore, the total mass of BHs in the 47 Tuc model is about the same as the mass of a potential intermediate-mass BH in the cluster (Kızıltan et al. 2017). However, since the pulsar accelerations in 47 Tuc can be reproduced with a group of stellar-mass BHs, the existence of an intermediatemass BH is not necessary (also see Mann et al. 2019). We will study how an intermediate-mass BH affects the properties of a cluster in future works.

Blue Stragglers
We define model BSSs based on the MS star turnoff (TO) luminosity, where the model MS stars reach the highest temperature. Temperatures of isolated stars are calculated using the stars' luminosities and radii (COSMIC outputs) assuming a black body radiation model, and temperatures of binary stars are calculated from the combined luminosity-weighted temperature of the component stars, where the total luminosity of a binary is the sum of the luminosities of its component stars (Weatherford et al. 2018, equation 1). To avoid contamination from normal MS stars, we specify BSSs as upper-MS stars with temperature T > T T O and luminosity L > 2L T O , either isolated or in binaries (e.g., Kremer et al. 2020). Parada et al. (2016) detected 114 BSSs using UV data taken from HST of 47 Tuc at 160 arcsec. In our models, we find a total of ∼ 101 BSSs per unit time, and ∼ 100 of them are within the theoretical half-mass radius (see Figure 8 and Table 2), less than the observed number of BSSs. This may be because COSMIC under-produces BSSs through the mass transfer channel (Leiner & Geller 2021). Most BSSs ( 95%) in the 47 Tuc models are formed in direct collisions of MS stars. In addition, the 2σ upper limit of the number of BSSs (Table 2) is very close to the observed number of BSSs, suggesting that the potential under-production of BSSs through the mass transfer channel is only a small effect in 47 Tuc.
The lifetime of BSSs in the best-fit model is also shown in Figure 9  2 Gyr. With the short lifetimes and the continuing formation of BSSs through stellar collisions, the average number of BSSs at late times of the cluster's evolution (Figure 9; bottom panel) is quite stable.

King Profile versus Elson Profile
Many previous studies have shown that a King profile (King 1966) may not be the best match for the observed SBPs of the majority of GCs in the Milky Way because of its sharp cutoff of the density distribution close to the tidal radius (e.g., McLaughlin & van der Marel 2005;Miocchi et al. 2013;de Boer et al. 2019;Hénault-Brunet et al. 2019). Instead, the Wilson profile (Wilson 1975), which is more spatially extended than the King profile, is shown to fit equally well or better where the observed SBPs exhibit a smoother decrease in the outermost radial density (e.g., GCs better than the King or the Wilson profile, including 47 Tuc (de Boer et al. 2019;Hénault-Brunet et al. 2019). LIMEPY models are isothermal in the central regions and are described by polytropes in the outer regions. They can reproduce King or Wilson profiles with certain choices of truncation parameter. Furthermore, it has been shown that young star clusters in the Large Magellanic Cloud are better described by an Elson profile with no apparent tidal truncation (Elson et al. 1987;Mackey & Gilmore 2003). The Elson profile is a generalization of the Plummer profile (Plummer 1911) that arises when γ = 4 in equation (8). The Elson profile is similar to the King profile in the inner regions of a cluster, but exhibits larger density or surface brightness at radii approaching the tidal radius (see Figure 4 in Mackey & Gilmore 2003). A typical range of γ for young star clusters in the Large Magellanic Cloud is ∼ 2.2−3.2 (Mackey & Gilmore 2003;Portegies Zwart et al. 2010), close to what we adopt for 47 Tuc.  Figure 10. Top 10 closest-fit models for 47 Tuc (gray curves) from the CMC Cluster Catalog ). Nine of them are different time steps from the same simulation, where initially it has a canonical mass function from Kroupa (2001), number of stars N = 1.6 × 10 6 , King concentration parameter W0 = 5, virial radius rv = 1 pc, and binary fraction f b = 5%. The metallicity is Z = 0.002 and the Galacticocentric distance is 8 kpc. The last model is from a simulation with initial rv = 2 pc and the other initial conditions the same as the above simulation.
in Figure 1). The top panels of Figure 1 and 11 show that models using both profiles can fit the observed SBP. However, the King profile model starts to deviate from observations at around 300 arcsec, which encloses ∼ 20% of the model's total mass. We also explored whether varying the tidal radius of the King models could improve their fit to the outer part of the SBP. Changing the Galactocentric distance from 47 Tuc's apocenter 7.4 kpc to its pericenter 5.5 kpc (Baumgardt et al. 2019), thereby shrinking the tidal radius of the cluster, does not affect the SBP. This indicates that models with a King profile are underfilling. Furthermore, we also ran a set of small models (initial masses ∼ 2 × 10 5 M ) with different initial King concentration parameters W 0 = 3, 5, 7, 9 in order to alter the cluster potential's shape. We find that unless W 0 is very large (  9), there is no apparent change in the final SBPs. More importantly, the model VDP does not match the observed VDP of  Figure 11. A massive CMC model with a King profile (the black curve and black markers). This model has initial number of stars N = 2×10 6 , initial King concentration parameter W0 = 7.5, and initial virial radius rv = 1.4 pc. Its initial binary fraction, IMF, metallicity and Galactocentric distance are the same as the best-fit 47 Tuc models (Table 1). These initial conditions are chosen such that the model SBP can match the observed 47 Tuc SBP.
47 Tuc (Figure 11, bottom panel). Figure 10 further demonstrates that a flatter density profile than the King profile is needed to explain the observed VDP at radii 50 arcsec. Taking all these into account, we adopt the Elson profile for better fits to the outer part of the cluster in this study.

Initial Mass Function
Whether or not the IMF is universal for GCs is still an ongoing debate. For example, Marks et al. (2012) found that stellar IMFs are more top-heavy (more high-mass stars) as a cluster's metallicity becomes lower and the gas density of the pre-GC cloud becomes larger. Later, Haghi et al. (2017) showed that top-heavy IMFs can better explain the relation between the observed massto-light ratios and the metallicities of GCs in the Andromeda galaxy. However, in contrast to the suggested initial power-law slope for the mass function of massive stars ( 1 M ) in 47 Tuc (Marks et al. 2012), both Giersz &Heggie (2011) andHénault-Brunet et al. (2019) found that steeper power-law slopes (comparing to Kroupa 2001) better match the 47 Tuc observations. The power-law slope we adopt for high-mass stars (α 2 = 2.8 for 0.8 ≤ m/M in Table 1)  We also ran simulations with α 2 = 2.3, 1.8 (the initial mass and other initial conditions are the same as in Table 1) to further check if a power-law slope for massive stars more consistent with the Kroupa IMF or a topheavy IMF is allowed. Figure 12 shows the SBP and VDP of a model at 10.55 Gyr with an initial α 2 = 2.3. Because the shallower IMF increases the number of massive stars and BHs, the model is puffier (having larger core radius from BH burning; e.g., Kremer et al. 2019) and does not fit either the observed SBP or VDP. The simulation with initial α 2 = 1.8 dissolved before 5 Gyr due to even more extreme BH burning (e.g., Chatterjee et al. 2010;Weatherford et al. 2021).
It is significant to point out that a top-light IMF may not be the only way to produce a good fit for 47 Tuc. There may be a degeneracy between having a top-light IMF or instead, having a smaller initial virial radius r v , both can reduce the number of BHs and the cluster core radius. A smaller initial r v corresponds to a shorter two-body relaxation time (Spitzer 1987) and a faster dynamical clock for the cluster ). Thus a model with smaller r v will process and eject its BHs faster through dynamical encounters in the same time (e.g., 13 Gyr). However, for 47 Tuc and similar simulations with large initial masses (∼ 2.5 × 10 6 M ), a simulation with smaller initial r v of 1 pc takes a much longer time than simulations with r v > 1 pc (> 2 months compared to 2 weeks). Therefore we do not explore small r v models in this paper, and leave them for future studies (e.g., Weatherford et al. in preparation).
Furthermore, Hénault-Brunet et al. (2019) found that a power-law with α 1 = 0.52 for m/M < 0.5 and α 2 = 1.35 for 0.5 < m/M < 1 matched their 47 Tuc models. This together with the IMF used for 47 Tuc modeling in Giersz & Heggie (2011), which we also adopted here, all have shallower power-law slopes for lower-mass stars than the Kroupa IMF (α 1 = 1.3 ± 0.5 for 0.08 ≤ m/M < 0.5 and α 2 = 2.3 ± 0.3 for 0.5 ≤ m/M < 1; Kroupa 2001). In addition, Baumgardt & Hilker (2018) derived a global power-law slope of 0.53 for 0.2 < m/M < 0.8 for 47 Tuc. This is slightly steeper than the 0.4 present-day mass function slope for ∼ 0.2 < m/M < 0.7 in our 47 Tuc model, which could be due to the selection bias against faint low-mass stars in the crowded central area of 47 Tuc. In Section 3 we show that, by adopting an IMF that is both bottomheavy and top-light, we can produce models that match the 47 Tuc observations. Although we have only explored a small set of initial parameters and did not vary the lower-mass power-law slope of the IMF because the 47 Tuc-size simulations are very computationally expensive. It is possible a power-law slope for lower-mass stars more consistent with that in the Kroupa IMF may be ale to produce a 47 Tuc model.

Implications for Black Hole Merger Rate
The top-light IMF is likely to reduce the total number of binary BH mergers (e.g., Weatherford et al. 2021). Over a Hubble time, the 47 Tuc simulation produced 220 binary BH mergers, including both those that merged in the cluster and those that merged after being ejected. To check how this steep power-law slope for massive stars affects binary BH mergers, we also calculate the number of mergers in the simulation with high-mass power-law slope α 2 = 2.3 (also see Figure 12). Surprisingly, this simulation produces ∼ 2700 BHs at late times, about 15 times more than in the best-fitting 47 Tuc models, but only 210 merging binary BHs over a Hubble time. This occurs because, although it produces more BHs, the model with α 2 = 2.3 is much less dense in the core than what we see in the best-fitting 47 Tuc simulation ( Figure 12). Thus its BHs undergo fewer dynamical encounters and there are fewer mergers per BH (e.g., Kremer et al. 2020, Figure 13). Hypothetically, how-ever, if we were to generate 47 Tuc models with a standard Kroupa IMF, but with smaller initial virial radius and a high central density at present, more binary BH mergers should be produced ). In Weatherford et al. (2021), models with a high-mass power-law slope α = 2.3 produce more binary BH mergers than those with α = 3. The difference may be attributed to the very different IMF, initial virial radius and/or initial density profile used in this study than in Weatherford et al. (2021).

Binary Fraction
Binary stars play an important role in the dynamical evolution of GCs (e.g., Heggie & Hut 2003). For example, BH binaries can act as energy sources in a cluster core, delaying the core collapse of clusters (e.g., Kremer et al. 2019). Learning the primordial binary fractions of GCs can also help reveal details of the formation of GCs and the properties of their pre-GC molecular clouds. Furthermore, since the majority of stars in the Milky Way galaxy are expected to be born in clusters, the binary fractions in GCs are also connected to that in the field (Lada & Lada 2003). However, the GC binary fractions, especially the relation between primordial and present-day binary fractions, are not well constrained (e.g., Fregeau et al. 2009;Milone et al. 2012, and references therein). Milone et al. (2012) found the total fraction of binaries in 47 Tuc to be ∼ 2%, and the binary fractions at mass ratio q > 0.5, 0.6, 0.7 to be 0.009 ± 0.003, 0.007 ± 0.003, 0.005 ± 0.003, respectively at the radial region between the cluster's observed core and half-light radii (Harris 1996(Harris , 2010. Ji & Bregman (2015) obtained a binary fraction in 47 Tuc to be 3.01 ± 0.13% within about the half-light radius using their direct couting method for binaries with high mass ratios q > 0.5. More recently, Mann et al. (2019) found the binary fractions in 47 Tuc for q in the ranges of 0.5 − 0.7, 0.7 − 0.9, and > 0.9 to be 0.0265, 0.0098 and 0.0039, respectively. These observations are all consistent with low binary fractions.
We "observed" the binary fractions in the best-fit model by identifying binary stars in its CMD using the technique in Rui et al. (2021), which mimics the observational cuts in Milone et al. (2012). We found that at radii between the observed core and half-light radii (Harris 1996(Harris , 2010, the best-fit model has binary fractions 0.0107 ± 0.00023, 0.0085 ± 0.00021 and 0.0063 ± 0.00018 at mass ratio q > 0.5, 0.6 and 0.7, respectively. This is consistent with the observation in Milone et al. (2012) at all three mass ratios. It has initial N = 2.95e6 and initial rv = 3.5 pc. See text for more details.
To explore how changing the primordial binary fraction affects the model, we also experimented with a cluster with initial binary fraction f b = 8%, initial number of stars N = 2.95e6 and initial virial radius r v = 3.5 pc. The small variation in the initial N and r v is to compensate for the increased number of binaries so that the cluster has similar initial mass to the best-fit 47 Tuc model and that its SBP and VDP at late times will fit the observations. Figure 13 compares the SBP and VDP of this model to the observations, and shows the model properties match well the observations. However, the model has binary fractions 0.0352 ± 0.00043, 0.0274 ± 0.00038, and 0.0199 ± 0.0032 at q > 0.5, 0.6, and 0.7, respectively, which are about three times higher than observed in Milone et al. (2012). Furthermore, models with 8% binary fraction may be overproducing MSPs and CVs (there are 132 MSPs and 787 CVs at 10.38 Gyr). Therefore, 47 Tuc is more likely to have a low primordial binary fraction of < 8% and closer to ∼ 2%.

Rotation in 47 Tuc
Many studies have shown that 47 Tuc is rotating and has velocity anisotropy using observations from HST, Gaia, and MUSE (e.g., Watkins et al. 2015;Bellini et al. 2017;Kamann et al. 2018;. These additional dimensions of the cluster's dynamical structure can improve our understanding of how the cluster forms and evolves, and may alter the cluster properties derived from the traditional fitting of SBPs with spherically symmetric profiles such as the King profile (King 1966) or the Elson profile (Elson et al. 1987). However, CMC assumes spherical symmetry (Rodriguez et al. 2021b) and is not able to treat rotation. Our model cluster's velocity dispersion is also in general isotropic 3 (see Figure 2). Compared to the structural parameters derived in Bellini et al. (2017) from their best-fit model that took into account rotation and asymmetry, our model has a larger mass but smaller observational core and half-light radii, which are more consistent with the radii in Harris (1996Harris ( , 2010. Although note that their model SBP (Bellini et al. 2017, Figure 9) is slightly lower in the core compared to Trager et al. (1995), which may contribute to the differences. Given our best-fit model matches the observations using a spherically symmetric Elson profile, the differences in the properties of 47 Tuc when taking into account rotation and velocity anisotropy are likely small.

CONCLUSIONS
In this study, we simulate 47 Tuc using the Monte Carlo N -body code CMC. As one of the most massive, densest and closest GCs in the Milky Way, there is a large data set of different 47 Tuc observations, but few N -body simulations because simulations with the cluster's mass and density are computationally time consuming. We produce 47 Tuc models that can simultaneously match the observed surface brightness, velocity dispersion and number density profiles, the pulsar acceleration profile, and the pulsar distribution in the cluster. These models are self-consistent and take into account all the relevant physics (e.g., BH dynamics) for the first time. The best-fit model has an IMF that is both bottom-heavy and top-light with power-law slopes α 1 = 0.4 and α 2 = 2.8 (Table 1). The models were also overfilling the tidal radius initially.
We report the number of exotic objects in the models that most closely match the observations; the bestfit models result in 54 MSPs, 14 LMXBs, 334 CVs, 186 BHs, 1368 NSs and 101 BSSs in the cluster at the present day. The numbers of MSPs, LMXBs and CVs are likely upper limits given that our simple treatments of tidal capture and giant star collision to maximize binary formation, and that we do not consider duty cycles and luminosities of LMXBs and CVs. We also demonstrate that tidal capture in GCs may be important to the formation of redback MSPs and many MSPs may 3 CMC can model clusters with anisotropic velocities (Rodriguez et al. 2021b), but we do not consider anisotropy in this study.
be formed from giant star collisions. Furthermore, an intermediate-mass black hole is not necessary to explain the pulsar accelerations observed in the cluster. We also discuss how different density distribution profiles (King vs. Elson profiles), IMFs, and binary fractions affect the observed properties of 47 Tuc, such as the binary BH merger rates. We point out that there may be a degeneracy in 47 Tuc models given a top-light IMF or an initial small virial radius, which have similar effects on reducing the number of BHs and thus the properties of the cluster. Due to the significant computational time taken to run a 47 Tuc-size simulation, we did not explore systematically the vast parameter space of initial conditions for star clusters, or the full range of uncertainties in our various assumptions and input parameters (e.g., characterizing the BH and NS natal kicks, or the treatment of binary mass transfer and common envelope phases). Therefore, we do not claim that our best-fit model presented here is unique. Software: CMC (Joshi et al. 2000(Joshi et al. , 2001Fregeau et al. 2003;Fregeau & Rasio 2007;Chatterjee et al. 2010Chatterjee et al. , 2013Umbreit et al. 2012;Morscher et al. 2015;Rodriguez et al. 2016Rodriguez et al. , 2021b, Fewbody (Fregeau et al. 2004 Figure 14. Flowchart illustrating single-single dynamical encounters in CMC. Note that if both stars are giants, the collision is performed when rp < 1.3(Rg1 + Rg2). Default f coll = 1 in CMC. "MS" is main-sequence stars, "CO" is compact objects, "RLOF" is Roche lobe overflow and "TDE" is tidal disruption events.

B. COMPARISON WITH DETAILED STELLAR MODELS
To verify that our simple application of polytropic stellar models agrees with more sophisticated calculations, we compare the oscillation energy and the critical capture pericenter distance between the tidal capture of two MS stars in this study with those in McMillan et al. (1987). Figure 15 shows the comparison. For η 10, the oscillation energy calculated with polytropic stellar models agrees well with the more detailed stellar models. At v ∞ ∼ 10 km s −1 , the velocity dispersion of typical GCs and of 47 Tuc, there can be as big as ∼ 20% discrepancy in the critical pericenter distances between different stellar models for low-mass MS stars, which is not likely to affect the outcomes of tidal capture (see Section 5). For simplicity, we assume that the structure of the stars does not change during tidal capture.

C. TIDAL CAPTURE AND GIANT STAR COLLISION UNCERTAINTIES
We are only considering tidal capture and giant star collision during single-single encounters in this study for simplicity. To estimate how the tidal capture and giant star collision rates will be affected if we include these interactions during binary-single and binary-binary encounters, we calculated the direct collision ("sticky-sphere" collision) rates for all types of stars in the CMC Cluster Catalog models . We found that for MS stars and  (1987). The left panel shows the tidal energy dissipation in terms of the dimensionless transit time. The right panel shows the ratio of the critical capture periastron separation to the radius of the star as a function of the relative velocity at infinity. "MS" is for MS stars and "evol" is for stars just evolve off the main-sequence.
compact objects, direct collisions occur ∼ 2.5 times more frequently during binary-single and binary-binary encounters than single-single encounters. However, for giant stars, direct collisions are ∼ 5 times more frequent in single-single encounters. Thus including giant star collisions in binary interactions is not likely to have a large influence in the overall rates, but the tidal capture rates is likely to increase by at least twice if the interaction is occuring in binarybinary and binary-single encounters. The implementation of tidal capture into Fewbody is currently undergoing (Vick et al. in preparation). We assume, during common-envelope evolution, α = 1 (parameter for scaling the efficiency of transferring orbital energy to the envelope) and λ = 0.5 (binding energy factor for stellar envelope) in binary evolution in COSMIC (Hurley et al. 2002;Breivik et al. 2020). These simple values give us a sense of how giant star collisions vary with different values for α and λ. For example, for a smaller binding energy factor where λ is close to zero, the release of recombination energy that can enhance envelope ejection is turned off (Claeys et al. 2014) and all giant star collisions lead to mergers. We point out that the values of α and λ are highly uncertain (e.g., Ivanova et al. 2013, and references therein). Previous studies suggested an α from as small as ∼ 0.2 (Zorotovic et al. 2010;Toonen & Nelemans 2013;Camacho et al. 2014), to as large as ∼ 5 (Fragos et al. 2019). In addition, while most treatments adopt a constant λ, it is more likely that λ varies according to the types of stars and the stellar mass fraction in their envelopes (Ivanova 2011). By varying α and λ and matching the number and properties of compact objects in GCs, especially pulsars, we could put better constraints on these two parameters but is out of the scale of this study.
Previous studies have shown that tidal-force induced oscillation during tidal capture may cause large fluctuation and instability in MS stars, and may lead to their disruption (e.g., Kochanek 1992). We compare the binding energy of the MS stars 0.15 GM 2 R with their induced oscillation energy during tidal capture (Kochanek 1992). Only four stars are disrupted for all R p , and none is disrupted when only considering Rp R < 1.3( M T 2M ) 1 3 (Kochanek 1992), where M , R is the mass and radius of the MS stars, respectively, and M T is the total mass of the two interacting objects. Figure 15 shows some discrepancies in the oscillation energy and critical radius calculated with polytropic stellar models and detailed stellar models (McMillan et al. 1987). For the oscillation energy induced by tidal force (left panel), the most pronounced difference is for n = 1.5 and 0.4 MS when η 6. In our 47 Tuc simulation, only about 4% of tidal capture encounters have η > 6, so the difference for η > 6 has negligible effect on the final tidal capture binaries. For the critical radius (right panel), the largest deviation is ∼ 20% again for n = 1.5 and 0.4 MS at v inf ∼ 10 km s −1 . However, these parameters mostly involve interactions between low-mass MS stars with another low-mass MS stars or WDs, and the discrepancy will not affect the number of tidal capture NS-MS star binaries. Note-Columns from left to right: number of objects, virial radius, power-law slope of an Elson profile or King concentration parameter, maximum radius (in virial radius) to sample the cluster, metallicity, tidal radius, binary fraction, E-Elson or K-King profile, masses and slopes of the IMF, and total cluster mass. *-the best-fit 47 Tuc simulation No. 2 and 3 also fit the SBPs and VDPs, so are used for calculating the more general uncertainties of the properties and numbers shown in Table 2.
No. 4 and 5 disrupted at early times of the simulation. The run of No. 8 was stopped at early times because it is very computationally expensive. The IMF for No. 7 and 8 are adopted from Hénault-Brunet et al. (2019). The tidal radius 171 pc is adopted from Lane et al. (2012).