Formation of Low-mass Black Holes and Single Millisecond Pulsars in Globular Clusters

Close encounters between neutron stars and main-sequence stars occur in globular clusters and may lead to various outcomes. Here we study encounters resulting in tidal disruption of the star. Using $N$-body models, we predict the typical stellar masses in these disruptions and the dependence of the event rate on host cluster properties. We find that tidal disruption events occur most frequently in core-collapsed globular clusters and that roughly $25\%$ of the disrupted stars are merger products (i.e., blue straggler stars). Using hydrodynamic simulations, we model the tidal disruptions themselves (over timescales of days) to determine the mass bound to the neutron star and the properties of the accretion disks formed. In general, we find that roughly $80-90\%$ of the initial stellar mass becomes bound to the neutron star following disruption. Additionally, we find that neutron stars receive impulsive kicks of up to about $20\,$km/s as a result of the asymmetry of unbound ejecta; these kicks place these neutron stars on elongated orbits within their host cluster, with apocenter distances well outside the cluster core. Finally, we model the evolution of the (hypercritical) accretion disks on longer timescales (days to years after disruption) to estimate the accretion rate onto the neutron stars and accompanying spin-up. As long as $\gtrsim1\%$ of the bound mass accretes onto the neutron star, millisecond spin periods can be attained. We argue the growing numbers of isolated millisecond pulsars observed in globular clusters may have formed, at least in part, through this mechanism. In the case of significant mass growth, some of these neutron stars may collapse to form low-mass ($\lesssim3\,M_{\odot}$) black holes.

1. INTRODUCTION Close stellar encounters have long been understood to be a prominent feature of dense star clusters. Encounters involving neutron stars have drawn particular interest since the 1980s when the first globular cluster millisecond pulsars (MSPs) were discovered (Lyne et al. 1987). In the standard picture, cluster MSPs are thought to form in low-mass X-ray binaries where the neutron star is spun up through accretion of material from its companion (e.g., Phinney & Kulkarni 1994). These binaries are expected to be assembled primarily through dynamical encounters where a neutron star is dynamically exchanged into a binary (e.g., Sigurdsson & Phinney 1995;Camilo & Rasio 2005;Ye et al. 2019). Once formed, the neutron star binary is hardened by In old core-collapsed clusters, the inner regions are dominated by massive main-sequence stars, white dwarfs, and neutron stars (e.g., Kremer et al. 2021b). The high densities within the centers of core-collapsed clusters lead naturally to an increased rate of close stellar encounters of these objects, and therefore, an increased rate of stellar collisions and tidal disruptions (e.g., Heggie & Hut 2003). A number of studies have shown that close tidal interactions involving specifically neutron stars and stars may lead to the formation of compact binaries (e.g., Fabian et al. 1975;Ray et al. 1987;Ivanova et al. 2005;Ye et al. 2021). Even closer encounters inevitably lead to collisions of neutron stars and stars (e.g., Krolik et al. 1984;Rasio & Shapiro 1991;Lombardi et al. 2006;Perets et al. 2016;Kremer et al. 2019b). Previous analyses have shown that such collisions may lead to common envelope-like events that may ultimately result in Thorne-Żytkow objects (Thorne & Zytkow 1977) or, if significant accretion and spin up occurs, MSPs (e.g., Davies et al. 1992;Davies & Benz 1995;Lee et al. 1996;Camilo & Rasio 2005).
Here we explore the possibility of MSP formation through accretion onto neutron stars following the tidal disruption of main-sequence stars. Our analysis proceeds in three steps: In Section 2, we analyze a suite of realistic N -body models that span the full parameter space of the Milky Way globular clusters to predict the rates and demographics of these tidal disruption events (TDEs). Using models tuned to NGC 6752, NGC 6624, and 47 Tuc, we predict the number of TDEs in clusters with known isolated MSP populations. Motivated by the N -body results, in Section 3 we perform hydrodynamic simulations of a few representative encounters to investigate the outcome of the TDEs. Finally, motivated by the hydrodynamics, in Section 4, we present a simple analytic model for the long-term evolution of the accretion disks formed. Incorporating key uncertainties, we predict the total mass and angular momentum accreted by the neutron star and address the key question of whether millisecond spin periods are attainable. We discuss our results and conclude in Section 5.

RATES AND DEMOGRAPHICS OF TDES FROM N-BODY MODELS
To compute the numbers of TDEs in typical clusters, we use the N -body models from our CMC Cluster Catalog (Kremer et al. 2020) which were computed using CMC (Rodriguez et al. 2022), a Hénon-type Monte Carlo code that includes prescriptions for various physical processes relevant to dense star clusters including two-body relaxation, stellar and binary evolution (computed using COSMIC; Breivik et al. 2020), and direct integration of small-N resonant encounters (Fregeau & Rasio 2007;Rodriguez et al. 2018). A number of parameters are varied within the CMC Catalog namely the initial cluster mass, initial virial radius, metallicity 1 , and radial position within the Galactic potential. Altogether, this catalog effectively spans the full parameter space of the Milky Way globular clusters and captures the formation of a variety of astrophysical objects such as gravitational-wave, X-ray binaries, (millisecond) pulsars, cataclysmic variables, and blue straggler stars.
In CMC, we record a TDE whenever a neutron star passes within the tidal disruption radius of a nearby main-sequence star: r T = (M NS /M ) 1/3 R , where M NS is the neutron star mass and M and R are the mass and radius of the star. In the case of M > M NS , r T < R and the relevant minimum pericenter distance for disruption is simply the stellar radius.
1 As we are specifically interested here in only those clusters that are sufficiently old to have undergone core collapse, we exclude the solar metallicity models published as part of the CMC Cluster Catalog and only examine TDEs occurring in models with Z = 0.01Z and 0.1Z , most typical of the Milky Way globular clusters. Note-Total number of TDEs identified in various models from the CMC Cluster Catalog. In column 2 we list the initial virial radius for each model (Kremer et al. 2019a;Ye et al. 2021). In columns 3-6 we list the present-day cluster mass, core radius, half-light radius, and central density (within core radius) of each model. In column 7, we denote whether the cluster model is core-collapsed at present. In column 10, we list the total number of TDEs in each model that lead to sufficient spin-up to form isolated MSPs (assuming s = 0.2 as discussed in Section 4) and, in brackets, we list the observed number of isolated MSPs, where relevant. In the rows 1-4, we list TDE counts for four "typical" cluster models with varying virial radius, rv. In rows 5-7, we show models that match well three specific clusters. In the last row, we list the total number of TDEs identified in the full cluster catalog and the inferred total number of TDEs expected in the full Milky Way globular cluster system (shown in parentheses).
In Figure 1 we show the stellar mass and disruption time for all TDEs identified in our suite of models. Black scatter points denote "normal" main-sequence stars with mass below the turn-off mass (indicated by the solid gray curve). Roughly 75% of identified TDEs fall into this category. Open blue circles denote stars with mass above the turn-off mass. These stars (which constitute roughly 25% of all identified TDEs) were formed from previous stellar collisions/mergers and would observationally be identified as blue straggler stars (BSSs; e.g., Sandage 1953). We also show as a dashed gray curve the boundary marking twice the turn-off mass. Stars to the right of this dashed curve were formed through two or more stellar collisions. As shown the disrupted stellar masses range from roughly 0.1 M (the assumed lower limit of the stellar mass function) to roughly 5 M . At very early times, there are also a handful of TDEs with M > 5 M , which are not plotted here since they are rare. The median mass of all disrupted stars is roughly 0.9 M (averaged over all time). For BSS (main-sequence star) TDEs, the median disrupted mass is roughly 1.4 M (0.7 M ). In Figure 1, we show as solid blue (black) curves the median mass of BSS (mainsequence star) TDEs versus time. Finally, the majority of TDEs ( 80%) occur at late times (t > 8 Gyr), after their host clusters' stellar-mass black hole populations have been mostly depleted (for discussion, see Kremer et al. 2020).
In Table 1 we list the total number of TDEs as well as BSS TDEs in four representative models from the CMC Cluster Catalog with initial stellar number of N = 8 × 10 5 , Galactocentric distance of 8 kpc, metallicity of Z = 0.1Z and four different initial virial radii from 0.5 − 4 pc. 2 As described in Kremer et al. (2019a), the initial virial radius (which sets the initial density of the model) determines whether or not the cluster undergoes core collapse by the present-day age (t ∼ 12 Gyr). For N = 8 × 10 5 , only models with r v = 0.5 pc have reached core collapse by this time. As a result, the r v = 0.5 pc model yields by far the most TDEs of the four.
In addition to the four representative models, we also list the CMC Catalog models that most effectively match the surface brightness and velocity dispersion profiles (following the method of Rui et al. 2021) of two corecollapsed clusters with large numbers of known isolated MSPs: NGC 6752 and NGC 6624 (Paulo Freire's pulsar catalog 2022). We also provide this same information for our CMC model for 47 Tuc (Ye et al. 2021). Although 47 Tuc is not core collapsed, it is sufficiently massive and dense to still yield a large number of TDEs (and observed isolated MSPs). In the last row of Table 1, we list the total number of TDEs in the full set of models, and the inferred total (in parentheses) for the full Milky Way cluster population determined by scaling up the set of models to match the total Milky Way cluster mass.
In typical core-collapsed clusters like NGC 6752, we predict an event rate of neutron star+star TDEs of up to roughly 10 Gyr −1 , averaged over the full (∼ 12 Gyr) lifetime. At late times (t > 8 Gyr) after core collapse has occurred, we find event rates of up to roughly 30 Gyr −1 per cluster. For massive non-core-collapsed clusters like 47 Tuc (and Terzan 5), we predict a TDE rate of roughly 2 Gyr −1 over the lifetime of the cluster. For lower-mass non-core-collapsed clusters (e.g., simulations 2-4 in Ta Table 2) and on bottom, we show a M = 1.2 M (a "blue straggler star" modeled as an Eddington standard model; simulation 6 in Table 2). In both simulations, rp = r T is adopted. In the M = 0.5 M (M = 1.2 M ) case, roughly 0.45 M (1.1 M ) is bound to the neutron star at the end of the simulation. Links to animations for these simulations (and others) are included in Table  2. The animation corresponding to the top (bottom) panel covers the simulation from t = −0.3 − 7.5 d (t = −0.5 − 6.3 d). In these two animations, the star is partially disrupted during the first pericenter passage. The partially stripped star becomes bound to the neutron star and returns for subsequent passages before ultimately being destroyed completely.
ble 1), we predict an event rate of at most 0.1 Gyr −1 . For our inferred Milky Way population (in which roughly 20% of globular clusters reach core collapse by present day; Harris 1996), we estimate a total TDE event rate of roughly 100 Gyr −1 . These rates are comparable to those found in previous studies. For instance, Sigurdsson & Phinney (1995) showed neutron star+star collisions occur most frequently in the densest clusters (n > 10 5 pc −3 ; comparable to that expected for core-collapsed clusters). In such clusters, this study estimated rates of roughly 100 per cluster lifetime, similar to our findings. Furthermore, this study predicted that roughly 75% of MSPs in these dense clusters should be single. Similarly, Davies & Hansen (1998) estimated neutron star collision rates of roughly 100 Gyr −1 for clusters with n > 10 5 pc −3 , corresponding to a rate of roughly 1000 Gyr −1 in the full Milky Way cluster population, consistent with our estimate given that they did not take into account the fact that core-collapsed clusters likely only spend a fraction of their lives in a core-collapsed state with extreme central densities. Finally, Hansen & Murali (1998) used the inferred MSP birth rate in clusters to estimate a neutron star+star collision rate of roughly 10 − 10 4 Gyr −1 in the Milky Way, also consistent with our estimate.
3. HYDRODYNAMIC EVOLUTION Motivated by the overall TDE demographics from our N -body models, we now explore the hydrodynamic evolution of these TDEs for a few representative cases. This work extends upon previous studies on this topic (e.g., Davies et al. 1992;Rasio 1993;Ivanova et al. 2005;Lombardi et al. 2006). In order to explore the hydrodynamic outcomes, we use the smoothed particle hydrodynamics (SPH) code StarSmasher (Rasio 1991;Gaburov et al. 2018). To model close encounters of neutron stars and main-sequence stars, we follow the method described in Kremer et al. (2022) where the neutron star is treated as a "point particle" and interacts with the 200k SPH particles of the star only via softened gravity. We perform eight SPH simulations, which are summarized in Table 2. In all cases, we adopt v ∞ = 0 km/s, representative of nearly parabolic encounters expected in typical globular clusters. Additionally, we assume a fixed neutron star mass of 1.2 M , typical for neutron stars formed through electron capture supernovae, expected to be the most common formation channel for neutron stars retained in globular clusters (e.g., Ye et al. 2019).
For the first three simulations, we model the encounter of a 0.5 M M-dwarf modeled as an n = 1.5 polytrope, governed by a polytropic equation of state with adia-batic index Γ = 5/3, interacting with a neutron star at three pericenter distances: r p = [0.5, 0.75, 1]×r T . These M = 0.5 M cases are representative of the low mass TDEs expected to occur (see Figure 1). Roughly 40% of the TDEs in our CMC models feature a stellar mass less than 0.8 M (i.e., M-and K-dwarfs). For the second set of three simulations, we model the case of a 1.2 M main-sequence star modeled as an Eddington standard stellar model (n = 3 polytrope index with an equation of state incorporating both ideal gas and radiation pressure, as in Lai et al. 1993), representative of the BSS TDEs shown in Figure 1. We also include two additional simulations (simulations 7 and 8 in the table) adopting M = 0.8 M and M = 2 M . The M = 0.8 M simulation explores the case of the disruption of a near turnoff mass main-sequence star representative of the most common disrupted star (see Section 2). The M = 2 M simulation explores the rarer case of even more massive stellar disruptions (less than 8% of the TDEs identified in Section 2 have stellar masses of 2 M or more). For each of these two extra simulations, we adopt r p = r T and the Eddington standard stellar model. Finally, the initial radii of the 0.5 M , 0.8 M , 1.2 M , and 2 M stellar models are 0.6 R , 0.8 R , 1.2 R , and 1.5 R , respectively.
The three pericenter distances chosen span a reasonable range of encounter types: for r p /r T = 0.5, the neutron star penetrates well into the stellar radius, while the r p /r T = 1 case is a more classic tidal disruption. We do not simulate even more distant encounters where little mass is stripped and stable neutron star+mainsequence star binaries may form through tidal capture (e.g., Fabian et al. 1975;Press & Teukolsky 1977), although these more distant encounters may well play an important role in the formation of binary MSPs (Ye et al. 2021). We discuss this further in Section 5.
In nearly all simulations performed, the star is partially disrupted during the first pericenter passage, becomes bound to the neutron star (i.e., is tidally captured), and is ultimately disrupted fully after one or more additional passages. The only exception is simulation 1, where the star is disrupted fully on the first pericenter passage due to the relatively high penetration depth of the encounter and relatively high mass ratio. As expected, the time between the initial pericenter passage and full disruption of the star (shown as ∆t in column 11 of Table 2) increases with increasing r p /r T .
In columns 4 and 5 of Table 2, we list the total mass bound to the neutron star and total mass ejected from the system after final disruption. In all cases, we find that roughly 80−90% of the initial stellar mass becomes bound to the neutron star, while the remaining roughly 10 − 20% of mass becomes unbound from the system. In columns 6, 7, and 8 we report, respectively, the total angular momentum of the material bound to the neutron star, characteristic disk radius, and characteristic disk angular frequency. These quantities are computed as in Equations (15), (16), and (17) of Kremer et al. (2022). In column 9, we report the characteristic viscous accretion timescale t v ≈ (αh 2 Ω disk ) −1 of the material bound to the neutron star (where h ≈ 1 is the disk scale height ratio and α ≈ 0.1 is the assumed disk viscosity). All quantities in columns 4-9 are reported at a time 2t orb (where t orb is the orbital period of material bound to the neutron star) after the final pericenter passage (when the star is completely disrupted). As discussed in Kremer et al. (2022), 2t orb is chosen to ensure a disk has had sufficient time to form. On longer timescales (t days), the hydrodynamic evolution is governed by the accretion process of the disk, which is not modeled in our SPH set up. Importantly however, the viscous accretion times estimated here are longer than the typical time elapsed between the first pericenter passage and the final disruption of the star (column 11 of Table 2). This indicates that significant accretion onto the neutron star is unlikely to occur on the timescales of our SPH simulations (t days), justifying our assumption to ignore accretion in the SPH modelling. We discuss the possible long-term outcome of the evolving disks in Section 4.
In column 10 of Table 2, we show the characteristic peak mass inflow rate of the disks, defined aṡ M = M disk /t v where M disk and t v are obtained directly from the SPH simulations (column 4 and 9, respectively). We discuss the accretion process further in Section 4. Finally, in column 12 of the table we report the final velocity of the neutron star. As described in Kremer et al. (2022) in the context of black hole TDEs, the compact object is expected to receive a dynamical "kick" as a result of the impulse from material ejected to infinity. For more penetrating encounters, the geometry of the unbound ejecta becomes increasingly asymmetric and as a result, the velocity of the kick increases for encounters that are more nearly head-on. For the TDEs modeled here, we find kick velocities ranging from roughly 2 − 20 km/s. While unlikely to be sufficient to eject the neutron stars from their host cluster, these velocities are likely sufficient to kick the neutron stars onto elongated orbits within their host. We return to this point in Section 5.
To illustrate the key effects, we show in Figure 2 the hydrodynamic evolution of models 3 and 6. We also provides links to videos of all simulations in the final column of Table 2. 3

ACCRETION, SPIN UP, AND FORMATION OF MILLISECOND PULSARS
Motivated by the disk features (e.g., mass, radius, and angular momentum) predicted from the SPH simulations, we now build a simple model for the disk evolution   Note-Outcomes of all SPH simulations. Columns 1 and 2 show initial stellar mass and radius. Column 3 shows pericenter distance of encounter in units of tidal disruption radius. We report all properties in columns 4-9 after the final pericenter passage.Ṁ = M disk /tv (column 10) shows the characteristic peak mass inflow rate. ∆t (column 11) denotes the time elapsed between the initial pericenter passage and the final pericenter passage (when the star is fully disrupted). In column 12 we list the final velocity (reported at the end of the simulation) of the neutron star (attained from the impulsive kick imparted by the asymmetric ejecta mass).
and discuss possible accretion rates onto the neutron star. For the material bound to the neutron star, we find typical masses of roughly 0.5 − 1 M and viscous accretion times of roughly 1−5 d (columns 4 and 9 of Table 2). This implies a characteristic mass inflow rate of roughly 30−300 M yr −1 (column 10 of Table 2), orders of magnitude above the standard Eddington accretion limit of roughly 10 −8 M yr −1 for typical neutron star assumptions (e.g., Phinney & Kulkarni 1994). The possibility of so-called "hypercritical" accretion has been considered in a variety of contexts including accretion during supernovae (e.g., Colgate 1971;Zel'dovich et al. 1972), as a feature of gamma-ray burst models (e.g., Popham et al. 1999;Lee et al. 2005), and for accretion in a common envelope (e.g., Chevalier 1993;Brown 1995;Terman et al. 1995;Fryer et al. 1996;Bethe & Brown 1998;Hansen & Murali 1998;MacLeod & Ramirez-Ruiz 2015a). The key is there must exist a way for the accretion energy to be radiated away without impeding the flow and limiting the mass accretion rate onto the neutron star. In principle, this may occur in several ways. In the case of disk-like accretion flows where the inward mass flow predominately occurs in the plane of the disk (similar to what is expected for off-axis tidal disruptions like those modeled with our SPH simulations described in the previous section), the accretion energy may flow relatively freely out low-density polar regions (e.g., Frank et al. 2002). Even for relatively spherical Bondi-Hoyle type (Bondi & Hoyle 1944) accretion flows, 4 where the geometry doesn't necessarily facilitate a low-density region through which accretion energy can escape unabsorbed, sufficiently high mass inflow rates render the Eddington limit irrelevant because photons can be trapped and advected inward with the accretion flow (e.g., Rees 1978;Begelman 1979). In this case, an accretion shock is expected to form near the neutron star surface within which the temperature and density are sufficiently high for neutrinos (produced through pair annihilation) to become the dominant cooling mechanism, carrying away the accretion energy without impeding the flow onto the neutron star (e.g., Colgate 1971). In the spherical limit, previous studies (e.g., Chevalier 1993Chevalier , 1996Brown 1995;Fryer et al. 1996) have shown that for mass transfer rates above a critical valueṀ cr ≈ 10 −4 M yr −1 , neutrino cooling allows for hypercritical accretion. For accretion flows with some rotational support,Ṁ cr may be slightly larger than the spherical case (e.g., Chevalier 1993), but hypercritical accretion is still expected to be possible (e.g., Chevalier 1996; Armitage & Livio 2000;Brown et al. 2000;Zhang & Dai 2009). The physics of hypercritical accretion for neutrino-dominated accretion disks has also been explored at length in the context of black hole accretion, especially in the context of gamma-4 Quasi-spherical accretion flows may be relevant for relatively head-on stellar collisions where the neutron star becomes embedded fully within the disrupted stellar envelope or even for off-center tidal disruptions if the viscous evolution of the bound material ultimately causes the initially disk-like structure to "puff up" and become quasi-spherical (e.g., Abramowicz et al. 1988;King & Begelman 1999  (1 M ) disk. Top panels show the accretion rate onto the neutron star, middle panels show the total (cumulative) mass accreted, and bottom panels show the total spin angular momentum accreted. Solid curves indicate where the mass transfer rate is sufficiently high to "smother" the magnetic field (in this case the field -assumed here to be 10 11 G -does not play an important role in the accretion flow) while dashed lines indicate the point in evolution where magnetic energy density dominates the accretion flow. Finally, as horizontal dotted lines in the bottom panels, we show neutron star spin angular momentum values for a few representative spin periods. ray burst models (e.g., Popham et al. 1999;Narayan et al. 2001;Di Matteo et al. 2002;Kohri & Mineshige 2002;Lei et al. 2009). In the black hole case, the critical mass transfer rate for neutrino cooling is estimated to be roughly 0.01 M s −1 , much larger than in the neutron star case since in the former case accretion energy can disappear into the black hole.
With the above considerations in mind, we assume in what follows that the accretion flow is disk-like (motivated by the outcomes of our hydrodynamic simulations) and that hypercritical (e.g., super-Eddington) accretion onto the neutron star is possible at all times.
In order to take into account potential mass outflows expected in these hypercritical accretion disks (e.g., Narayan & Yi 1995;Blandford & Begelman 1999;Metzger et al. 2008), we parameterize the mass inflow rate as a power-law in radiuṡ where R acc is the accretion radius (i.e., the radius of the inner edge of the disk) and t v is the viscous accretion time. The exponent s ∈ [0, 1] parameterizes the uncertain amount of material transported from the outer edge of the disk (near the tidal disruption radius) to the accretion radius (Blandford & Begelman 1999). For example, for s ≈ 0.5 (e.g., Yuan et al. 2012) and for R acc = 10 6 cm and R disk ≈ 10 11 cm (Table 2), (R acc /R disk ) 0.5 ≈ a few × 10 −3 . In this case, only a fraction of the disk mass is actually accreted and the remainder is ejected via a disk wind. For the high magnetic field strengths expected in neutron stars, magnetic stresses can in principle dominate the flow in the accretion disk, especially near the neutron star surface where the field is strongest. The Alfvén radius, r A , is the characteristic distance from the neutron star at which the magnetic energy density is equal to the kinetic energy density of the material in the disk (e.g., Shapiro & Teukolsky 1983): (2) where µ = BR 3 NS is the assumed magnetic moment and where we have adopted M NS = 1.2 M and R NS = 10 6 cm. For high mass inflow rates representative of s 0.2, the Alfvén radius lies within the neutron star radius for field strengths expected for old neutron stars found in typical globular clusters (e.g., B 10 11 G; Ye et al. 2019). In this limit, the magnetic field is "smothered" by the large inflow of mass and the magnetic stresses do not play an important role in the accretion process. For lower mass transfer rates (or higher field strengths), the disk is truncated at r A , and the accretion flow onto the neutron star surface is dominated by the magnetic field (i.e., the accretion flows onto the neutron star along the field lines in the magnetosphere).
Following Metzger et al. (2008), the time dependence of the accretion rate for thick disks parameterized by Equation (1) can be expressed aṡ where C = 2s/(2s + 1) and M d,i , R d,i , and t v,i are the initial disk mass, disk radius, and viscous accretion time, respectively. 5 R acc is the accretion radius, which we define as the maximum of R NS and r A at a given time. The total mass accreted by the neutron star after time t can be computed by integrating Equation (3). The time derivative of the total accreted angular momentum can be expressed aṡ where M acc is the total mass accreted after time t and M NS is the (evolving) neutron star mass. We have assumed thatṀ NS =Ṁ acc (sinceṀ acc already incorporates implicitly outflows associated with the s parameter, this assumption is equivalent to stating simply that all material successfully transported to the neutron star surface is accreted) and that the accretion radius R acc is roughly constant in time (especially appropriate for r A < r NS when most of the mass is accreted). Integration of Equation (4) gives the total angular momentum supplied to the neutron star after time t.
In Figure 3, we show the accretion rate, total accreted mass, and total accreted angular momentum for the neutron star versus time (with t = 0 defined as the time of disk formation, i.e., roughly the end of the SPH simulations discussed in Section 3). As different colors, we show different assumed values of the uncertain s parameter, ranging from s = 0 (the case of highest mass inflow rate) to s = 1 (a much lower mass inflow rate case, in which the majority of disk mass is blown away in a wind). In the left column, we show the evolution for an initial disk mass M d,i = 0.4 M , initial disk radius R d,i = 9 R , and viscous accretion time t v,i = 5 d, representative of the tidal disruption of a 0.5 M Mdwarf (e.g., simulation 3 in Table 2). For the plots in the right column, we assume M d,i = 1 M , R d,i = 3 R , and t v,i = 2 d, typical of the 1.2 M BSS TDE case (e.g., simulation 6). In all cases, we assume a magnetic field strength of 10 11 G (e.g., Ye et al. 2019). Solid curves indicate evolution where the mass transfer rate is sufficiently high for the Alfvén radius to lie within the neutron star radius. In this case the magnetic field is "smothered" and does not play an important role in the accretion flow. Dashed lines indicate the point in evolution where magnetic energy density dominates the accretion flow so r A > R NS . In reality, the neutron star magnetic field may decrease (be "buried") through the accretion process (e.g., Bhattacharya & van den Heuvel 1991), in which case the magnetic field may not influence the accretion flow until even later times still (at even lower accretion rates).
For a neutron star with moment of inertia 2M NS R 2 NS /5 and spin period P s , the spin angular momentum is Assuming the majority of accretion occurs at the neutron star radius (appropriate for low s cases as shown in Figure 3) and M NS M acc , the total accreted angular momentum can be expressed simply as J acc ∼ M acc √ GM NS R NS . In this case (assuming the initial spin period is roughly zero), the spin period attained through accretion is given roughly by (6) Thus, the available disk mass is sufficient to spin up the neutron star to millisecond spin periods. The exact spin period attained depends on the efficiency of the accretion. We show as horizontal dotted lines in the bottom panels of Figure 3 the spin angular momentum values corresponding to a few characteristic neutron star spin periods (for M NS = 1.2 M ). As shown, in the high inflow case of s 0.2, a MSP can be produced for both disk masses assumed here. 6 For s 0.5, the neutron star is unlikely to be spun up significantly, regardless of the disk mass.
Of specific relevance is the work by MacLeod & Ramirez-Ruiz (2015a) who evaluated the mass growth of a neutron star embedded within a common envelope (taking into account the asymmetric structure of the structure of the envelope; e.g., MacLeod & Ramirez-Ruiz 2015b). Although a common envelope is not exactly identical to the accretion geometry expected for the tidal disruptions considered here, there are key qualitative similarities. This study found that, in general, a modest amount ( 0.1 M ) of the envelope material is 6 In the most extreme case of s = 0 for M d,i = 1 M (yellow curve in right hand panel of Figure 3), the total angular momentum accreted corresponds to a neutron star spin period of roughly 0.4 ms, near the break-up angular velocity of the neutron star (corresponding to the Keplerian velocity at the neutron star surface). In reality, once the break-up velocity is reached the neutron star is unlikely to be spun up further (although it may continue to accept mass). the fraction of isolated MSPs observed versus core radius over half-light radius for various Milky Way globular clusters. Open stars (circles) denote clusters that have (have not) undergone core collapse. The size of the points is scaled by the total number of MSPs (single or binary) observed in that cluster. On the righthand vertical axis (black points), we show the total number of neutron star+star TDEs occurring in our various CMC Catalog cluster models (again black stars and circles denote core-collapsed and non-core-collapsed clusters, respectively). We argue the overabundance of isolated MSPs observed in the most centrally concentrated clusters can be explained in part by the increased rate of TDEs in these systems.
accreted by the neutron star. They pointed out that such mass growth is likely sufficient to spin up the neutron star (consistent with our result here), but is likely insufficient to lead to collapse to a black hole (a point we return to in Section 5). Figure 4 summarizes the key result of this study. On the left-hand vertical axis (blue color) we plot the fraction of observed isolated MSPs relative to the total number of observed MSPs (defined here as having spin periods less than 30 ms) versus core radius, r c , over halflight radius, r h , for all relevant Milky Way globular clusters. Blue stars indicate clusters that have undergone core collapse and blue circles indicate clusters that have not. The sizes of all blue points are scaled by the total number of MSPs observed in the cluster. On the righthand vertical axis (black), we plot the total number of neutron star+main-sequence star TDEs versus r c /r h for all relevant models from the CMC Cluster Catalog (again stars versus circles denote core-collapsed versus non-core-collapsed). As demonstrated in the figure, the most centrally dense clusters, especially those that have undergone core collapse, feature the highest fraction of isolated MSPs and the highest rate of TDEs.

DISCUSSION & CONCLUSIONS
The choice of 30 ms to define a MSP is consistent with previous studies (e.g., Lorimer 2008;Ye et al. 2019), but is admittedly somewhat arbitrary. Some core-collapsed clusters (e.g., M15 and NGC 6624) contain a handful of mildly recycled isolated pulsars with relatively large spin periods ( 100 ms) that in principle may also have formed through TDEs but with less mass accreted by the neutron star (e.g., Camilo & Rasio 2005). Overall, the median spin period for all isolated pulsars in clusters is roughly 5.3 ms (Paulo Freire's pulsar catalog 2022). Less than 10% of the observed isolated pulsars have periods in excess of 30 ms, thus this definition reasonably captures the bulk of the observed distribution. Additionally, we note that the spin period distribution for the observed binary pulsars in clusters peaks at a slightly lower value than the isolated pulsars (the median spin period of the binary pulsars is roughly 3.7 ms) and a Kolmogorov-Smirnov test reveals these two distributions may in fact be distinct (KS statistic of roughly 0.3). Tentatively, this may hint at different formation channels for the binary versus isolated pulsars (e.g., classic binary mass transfer versus TDEs).
As shown by Equation (6), the final spin period of the neutron star is determined by the mass accreted, which in turn is determined by the mass bound to the neutron star and the accretion efficiency of the disk (the s parameter in Equation 1). For simplicity, we can assume that, as shown in our SPH models, roughly 90% of the disrupted star becomes bound to the neutron star. Under this assumption, from the distribution of stellar masses that undergo TDEs in our CMC models, we estimate that, for s = 0.2 (assuming R disk ≈ 10 11 cm so that [R NS /R disk ] 0.2 ≈ 0.1), in all TDEs the neutron star would accrete sufficient mass to attain MSP spin periods (P s < 30 ms), corresponding to an average of roughly 10 isolated MSPs per cluster. For s = 0.4 ([R NS /R disk ] 0.4 ≈ 0.01), roughly 56% of TDEs (558 out of 988) would create MSPs, corresponding to roughly 5 isolated MSPs per cluster. For s = 0.5 (roughly 0.3% of the bound mass is accreted), only 41 (roughly 4%) of all identified TDEs would lead to MSPs.
Current observations have revealed 98 isolated MSPs in the Milky Way globular clusters. Of these, 31 are observed in six core-collapsed clusters, 29 in the massive non-core-collapsed clusters Terzan 5 and 47 Tuc, and the remaining 38 are found in eleven lower-mass noncore-collapsed clusters. Because of various observational biases, this sample is likely to remain highly incomplete and the true number of isolated (as well as binary 7 ) MSPs could be much higher. In this case, the formation of roughly 5-10 isolated MSPs per cluster overall (including roughly 20-100 per typical core-collapsed cluster, roughly 20 per massive non-core-collapsed cluster like 47 Tuc or Terzan 5, and roughly 1 per typical low density cluster; see column 9 of Table 1) suggested by efficient disk accretion models is quite possibly consistent with observations. As discussed in Section 3, regardless of the accretion and spin-up process, the neutron stars are expected to receive impulsive kicks of up to roughly 20 km/s from the asymmetric ejection of material stripped from the star during disruption. As a consequence of these kicks, we predict isolated MSPs formed through these TDEs should, on average, be found at large offsets from their host cluster's center. The average radial position (in units of their host cluster's core radius) of all the observed isolated MSPs in clusters with known radial positions is roughly 2.4 (Paulo Freire's pulsar catalog 2022) -these objects are clearly offset from their hosts' centers, which may be indicative of the velocity kicks proposed here. For reference, this value for binary MSPs with known cluster positions is roughly 2, marginally lower than the isolated MSP value. We reserve for future study a detailed comparison of the observed offset distribution and the predicted offset distribution expected for the velocity kicks predicted by our models.
Although we argue formation of isolated MSPs is a plausible outcome of neutron star TDEs, it is certainly not the only possibility. For inefficient accretion disks (s 0.2) where only a small fraction of mass is accreted, the spin angular momentum of the neutron star will only increase slightly (Equation 6). In this case, the TDEs would have a negligible effect upon the properties of the disrupting neutron stars even for the most massive disrupted stars. On the other hand, for highly efficient accretion disks (s ∼ 0) that are sufficiently massive (M disk 1.5M ), the neutron star may accrete sufficient material to exceed the (uncertain) maximum allowable neutron star mass and, in this case, may collapse to form a low-mass black hole. The possibility of neutron stars being driven to collapse through accretion in a stellar envelope has been explored in the context of common envelope evolution of binaries (e.g., Chevalier 1993;Bethe & Brown 1998;Armitage & Livio 2000;Bethe et al. 2007;MacLeod & Ramirez-Ruiz 2015a). In the specific case of a neutron star colliding with a massive main-sequence star, where the final collision product qualitatively resembles a collapsar (Woosley 1993), Hansen & Murali (1998) showed that the collapse to a black hole may be accompanied by a (long) gamma-ray burst. For the case s = 0 where the entire disk mass is accreted, we find that 39% (16%) of the TDEs in our CMC catalog models would lead to collapse to a low-mass black hole, assuming a maximum neutron star mass of 2 M (2.5 M ). This translates to roughly 1 − 4 lowmass black holes formed per typical core-collapsed clus-ter (in this case, the number of MSPs quoted previously would be reduced slightly since 16 − 39% of the MSPs would instead become black holes). For s 0.2, where less than 10% of the disk mass is accreted, only one of the TDEs identified in our CMC models would lead to low-mass black hole formation. Therefore this outcome appears significant only if nearly the entire mass bound to the neutron star can be accreted.
A key process not considered here is the potential role of feedback energy in unbinding material initially bound to the neutron star (with binding energy E bind ≈ GM NS M disk /R disk ≈ 10 48 erg). In principle, material may be unbound before the roughly 10 −2 M necessary to attain millisecond spin periods can be accreted by the neutron star, thus inhibiting the viability of these TDEs as a MSP formation mechanism. Feedback may arise through accretion energy (e.g., Armitage & Livio 2000;Papish et al. 2013) or nuclear energy generated through burning of hydrogen (and possibly heavier elements) near the neutron star surface (e.g., Hansen & van Horn 1975;Taam 1985;Bildsten 1998). Energy generated through accretion is expected to be of order E acc ∼ ηM acc c 2 , where η is the uncertain accretion efficiency. For a typical η ≈ 0.01, M acc ≈ 10 −4 M is sufficient to unbind the envelope, assuming the accretion energy can very efficiently couple mechanically with the envelope. In reality, for disks similar to those considered here, a fraction of the accretion energy can likely be released relatively unabsorbed through the polar regions (e.g., a jet-like geometry; Livio 1999).
In the nuclear energy case, previous studies have demonstrated in the context of common envelope episodes, nuclear energy may be sufficient to eject remaining bound material (e.g., Podsiadlowski et al. 2010;Ivanova et al. 2015;Grichener et al. 2018). However, the efficiency of mechanical coupling is also key here; if the coupling is inefficient and most of the nuclear energy can be released as radiation (e.g., Grichener et al. 2018;Soker et al. 2018) then ejection of the envelope may be difficult. We reserve for future study treatment of the possible role of feedback from both accretion and nuclear energy on the long-term outcome of these TDEs. Given that the cross section for close encounters scales linearly with stellar radius in the parabolic regime, the occurrence of TDEs, with r p < r T , implies a comparable number of more distant encounters with r p in the range from r T to a few × r T that may form long-lived binaries through tidal capture (Fabian et al. 1975). In Ye et al. (2021), we argued these tidal captures may eventually lead to the formation of "redback" MSPs (e.g., Strader et al. 2019), provided the companion star fills its Roche lobe and transfers mass onto (and spins up) the neutron star. This would imply an overabundance of redback MSPs in core-collapsed clusters, for the same reasons we argue tidal disruptions lead to an overabundance of isolated MSPs. There are 16 redbacks currently known in Milky Way globular clusters, four of which are in core-collapsed clusters. 8 Given that only 20% of Milky Way clusters are core-collapsed (Harris 1996), this perhaps suggests a marginal overabundance. From a hydrodynamic perspective, it remains unclear whether the ultimate fate of tidal captures is indeed the formation of a detached binary (e.g., Camilo & Rasio 2005) and, if so, whether the amount of mass transferred, is sufficient to spin up the neutron star to millisecond periods. Alternatively, depending on how quickly the tidally distorted star can radiate away the dissipated tidal energy, tidal captures may lead ultimately to mergers as we see clearly for closer encounters. In that case, the numbers of isolated MSPs predicted here may increase by a small factor. We reserve for future work more careful consideration of the distinction between tidal disruptions and captures and the implementation of a self-consistent treatment of the formation and fate of MSPs through tidal disruptions and captures within CMC.
As summarized in Kremer et al. (2021a), the neutron star+main-sequence star TDEs considered here are just one of several possible processes expected in corecollapsed clusters that could in principle create rapidly spinning neutron stars. Accretion and spin-up may similarly occur for tidal disruptions of white dwarfs by neutron stars. As discussed in Metzger (2012), these TDEs may lead to orders-of-magnitude larger mass transfer rates. Concerning the event rates, on the one hand, the cross section for white dwarf tidal disruptions is a factor 100 times smaller than for main-sequence stars (accounting for the relatively tiny radii but relatively high masses of white dwarfs compared to typical cluster main-sequence stars); on the other hand, white dwarfs are expected to be far more abundant in the inner regions of core-collapsed clusters (e.g., Rui et al. 2021). With this in mind, Kremer et al. (2021a) showed that the total rate of white dwarf + neutron star TDEs is roughly a factor of 10 lower than the main-sequence star TDE rate. Alternatively, MSPs may be produced by mergers of pairs of white dwarfs (e.g., Schwab 2021; Kremer et al. 2021b). Depending on various physical processes, white dwarf mergers may alternatively lead to Type Ia supernovae (e.g., Webbink 1984), slowly spinning pulsars possibly connected to the "young pulsars" observed in several Milky Way globular clusters (e.g., Tauris et al. 2013), or magnetars (e.g., King et al. 2001), which may also connect with fast radio bursts similar to FRB 20200120E (e.g., Bhardwaj et al. 2021;Kremer et al. 2021b;Lu et al. 2022). In an upcoming paper (Ye et al., in preparation) we will implement the for-mation of MSPs through neutron star + main-sequence star TDEs and other aforementioned mechanisms within CMC, enabling us to track self-consistently the formation and subsequent dynamical evolution of these objects.

ACKNOWLEDGMENTS
We thank the anonymous referee for their careful review of the paper. We also thank Tony Piro, Phil Hopkins, Nick Kaaz, and Ariadna Murguia-Berthier for helpful discussions. KK is supported by an NSF Astronomy and Astrophysics Postdoctoral Fellowship under award AST-2001751. FK acknowledges support from a CIERA Board of Visitors Graduate Fellowship. SMR is a CIFAR Fellow and is supported by the NSF Physics Frontiers Center awards 1430284 and 2020265. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. This work was supported by NSF Grant AST-2108624 and NASA ATP Grant 80NSSC22K0722 at Northwestern University.