Tidal Disruptions of Stars by Black Hole Remnants in Dense Star Clusters
Abstract
In a dense stellar environment, such as the core of a globular cluster (GC), dynamical interactions with black holes (BHs) are expected to lead to a variety of astrophysical transients. Here we explore tidal disruption events (TDEs) of stars by stellarmass BHs through collisions and close encounters. Using stateoftheart body simulations, we show that these TDEs occur at significant rates throughout the evolution of typical GCs and we study how their relative rates relate to cluster parameters such mass and size. By incorporating a realistic cosmological model of GC formation, we predict a BH – mainsequencestar TDE rate of approximately in the local universe () and a cosmological rate that peaks at roughly for redshift 3. Furthermore, we show that the ejected mass associated with these TDEs could produce optical transients of luminosity with timescales of about a day to a month. These should be readily detectable by optical transient surveys such as the Zwicky Transient Facility. Finally, we comment briefly on BH – giant encounters and discuss how these events may contribute to the formation of BH – whitedwarf binaries.
0000000240863180]Kyle Kremer \move@AU\move@AF\@affiliation Department of Physics & Astronomy, Northwestern University, Evanston, IL 60202, USA \move@AU\move@AF\@affiliation Center for Interdisciplinary Exploration & Research in Astrophysics (CIERA), Evanston, IL 60202, USA 0000000215687461]Wenbin Lu \move@AU\move@AF\@affiliationTAPIR, Walter Burke Institute for Theoretical Physics, Mail Code 35017, Caltech, Pasadena, CA 91125, USA 0000000341758881]Carl L. Rodriguez \move@AU\move@AF\@affiliationPappalardo Fellow; MITKavli Institute for Astrophysics and Space Research, Cambridge, MA 02139, USA \move@AU\move@AF\@affiliationDepartment of Physics, Allegheny College, Meadville, Pennsylvania 16335, USA \move@AU\move@AF\@affiliation Center for Interdisciplinary Exploration & Research in Astrophysics (CIERA), Evanston, IL 60202, USA 000000027132418X]Frederic A. Rasio \move@AU\move@AF\@affiliation Department of Physics & Astronomy, Northwestern University, Evanston, IL 60202, USA \move@AU\move@AF\@affiliation Center for Interdisciplinary Exploration & Research in Astrophysics (CIERA), Evanston, IL 60202, USA
1 Introduction
The high stellar densities at the centers of globular clusters (GCs) make them hotbeds for exotic astrophysical objects, including blue stragglers (e.g., Sandage1953; Ferraro1995; Piotto2003; Geller2011; Leigh2013; Chatterjee2013b), millisecond pulsars (e.g., Lyne1987; Sigurdsson1995; Ransom2008; Ivanova2008; Fragione2018a; Ye2018), and lowmass Xray binaries (e.g., Clark1975; Verbunt1984; Heinke2005; Ivanova2013; Giesler2018; Kremer2018a). The influence of dynamics is particularly pronounced for stellarmass black holes (BHs), which preferentially reside in the highestdensity central regions of their host clusters, as a consequence of mass segregation. Recent work has shown that GCs efficiently produce merging binary BHs (BBHs) through series of binarymediated dynamical encounters (e.g., Moody2009; Banerjee2010; Bae2014; Ziosi2014; Rodriguez2015a; Rodriguez2016a; Askar2017; Banerjee2017; Giesler2018; Hong2018; Fragione2018b; Samsing2018a; Rodriguez2018b). Thus, as LIGO and Virgo have begun to unveil the gravitationalwave universe through detections of merging BBHs (e.g., Abbott2016a; Abbott2016c; Abbott2016b; Abbott2016d; Abbott2017; Abbott2018b), the topic of binary BH formation in GCs has risen in importance.
In the past two decades, our understanding of stellarmass BH populations in GCs has evolved significantly. From an observational perspective, several BH candidates have been identified in GCs through both Xray/radio (e.g., Maccarone2007; Strader2012) and radialvelocity (Giesers2018) observations. Recent computational work has supported these observational findings. Current models of GCs show that they can retain large populations of BHs, even at late times ( Gyr) and the BHs readily mix with luminous stellar populations forming BH – nonBH binaries. Furthermore, it is now understood that the BHs can play a dominant role in determining the largescale structural properties of their host cluster (e.g., Mackey2007; Mackey2008; Kremer2018b; Kremer2019a).
In a GC, stars will frequently pass through small (most commonly, three or fourbody) resonant encounters that lead to a greatly increased rate of close passages between objects relative to direct single–single interactions (e.g., Bacon1996; Fregeau2004). If these small encounters involve BHs and other stars, sufficiently close passages can lead to tidal disruption events (TDEs) of stars by the BHs. Previous work has shown that these events may serve as a dynamical formation channels for a variety of sources observed in GCs. For instance, Ivanova2005 showed that direct physical collisions between neutron stars (NSs) and giant stars can lead to bright, ultracompact Xray binaries, similar to those observed in Milky Way GCs, in which the NS is accreting from a whitedwarf (WD) companion. Similarly, Ivanova2010 showed that collisions between BHs and giants may result in accreting BH – WD binaries, possibly similar to the Xray/radio source observed in NGC 4472 (Maccarone2007). More recently, Ivanova2017 showed that (grazing) tidal captures of subgiants by stellarmass BHs may lead to BH – subsubgiant binaries, similar to the BH candidate observed in M10 (Shishkovsky2018).
Additionally, TDEs involving stellarmass BHs may produce interesting electromagnetic transient signatures. Perets2016 showed that disruptions of stars by stellarmass BHs and NSs, which they called “microTDEs,” may give rise to bright, energetic, long Xray/gammaray flares, which could possibly resemble ultralong gammaray bursts. More recently, Lopez2018 demonstrated that when a star is disrupted by a BBH, the BH spins may be significantly modified if a significant fraction of stellar material is accreted. Samsing2019 showed that for disruptions by BBHs, observations of breaks in the light curve may be used to determine the BBH orbital period, and they proposed to link this inferred orbital period to the BBH merger rate from clusters, thus connecting tidal disruption and gravitationalwave signals.
In this analysis, we use our Hénontype Monte Carlo code CMC to explore the rates of TDEs by stellarmass BHs in GCs, and we then estimate the electromagnetic signatures expected from these events. Our paper is organized as follows. In Section 2 we summarize the methods used to model the longterm evolution of GCs, including our treatment of BH TDEs. In Section 3 we discuss the total number of BH – mainsequencestar (MS) TDEs identified in our models and discuss the dependence of these events on cluster properties, in particular the initial cluster size and mass. In Section 4 we discuss the expected electromagnetic signatures of such events and in Section 5 we give our predictions for the event rate. We briefly discuss BH–giant encounters in Section 6. Finally, we discuss our results and conclude in Section 7.
2 Globular cluster models
\twocolumngridCMC (for Cluster Monte Carlo) is a rigorously tested, Hénontype Monte Carlo code that computes the longterm evolution of GCs with realistic initial properties and containing up to several million stars (Henon1971a; Henon1971b; Joshi2000; Joshi2001; Fregeau2003; Pattabiraman2013; Chatterjee2010; Chatterjee2013; Rodriguez2015a). CMC incorporates stellar and binary evolution using the evolution codes SSE and BSE (Hurley2000; Hurley2002), modernized to reflect our most updodate understanding of compact object formation, including prescriptions for natal kicks, massdependent fallback, and (pulsational) pairinstability supernovae (Fryer2001; Belczynski2002; Hobbs2005; Morscher2015).
We use the same prescription to treat pulsationalpair instabilities and pairinstability supernovae as Rodriguez2018a, which follows Belczynski2016b. Briefly, we assume that any star with a preexplosion helium core between 45 and 65 will undergo pulsations that eject a significant amount of mass, until the final product is at most 45 . Assuming that of the mass is lost during the conversion from baryonic to gravitational mass at the time of collapse, this leaves a BH of mass . Finally, we assume that stars with helium core masses in excess of are completely destroyed by a pairinstability supernova (no remnant is formed). Of course, there are many uncertainties associated with these numbers (in particular, some studies, (e.g., Woosley2016; Spera2017) have predicted that the lower boundary of this mass gap may be as high as ).
For NS formation and evolution, we have updated the SSE and BSE codes to include changes to the magnetic field and spinperiod evolution, as well as the natal kick prescriptions for NSs formed in ECSNe (Kiel2008; KielHurley2009). See Ye2018 for further detail on the treatment of NSs and radio pulsars in CMC.
In a GC, the stellar dynamics in the highdensity core can frequently lead to close encounters between stars. For encounters of two single stars (single–single interactions) with masses and , the crosssection for an encounter with pericenter distance is given by
(1) 
where is the local velocity dispersion (see, e.g. Fregeau2007). Here, we are interested in the particular case of close encounters between BHs and MS stars. For sufficiently close encounters, the tides raised on the MS star may dissipate enough energy to bind the two objects or disrupt the structure of MS star. Following Fabian1975, we define the characteristic radius for tidal interaction as
(2) 
where is the BH mass, and are the stellar mass and radius, respectively, and is a dimensionless parameter that captures the details of the tidal interaction. In order of increasing , close encounters may take the form of physical collisions (e.g., FryerWoosley1998; Hansen1998; ZhangFryer2001), tidal disruptions (e.g., Perets2016), tidal captures (e.g., Fabian1975; Ivanova2017), or more distant tidal encounters that only weakly perturb the star (e.g., AlexanderKumar2001). In addition to the dependence on the pericenter distance, the details of a particular close encounter and associated tidal interactions also depend upon the internal structure of the star. For instance, tidal interactions can be quite different for young lowmass MS stars, which have relatively uniform density, compared to evolved massive MS stars which are much more centrally concentrated.
Capturing the subtleties of tidal interactions during close encounters requires detailed hydrodynamic calculations of the events, which is beyond the computational scope of CMC. Therefore, we simply consider two limits for our choice of the maximum pericenter distance, , that leads to a TDE. As a conservative lower limit, we take (models 124 in Table 2 adopt this assumption) and as an upper limit, we assume computed by equation 2 with (models 2526 adopt this assumption). Henceforth, we use the term TDE as shorthand to denote all close encounters with , computed in one of these two limits.
Although BH – MSstar TDEs are the main focus of this study, we comment briefly on NS – MSstar TDEs in Section 7.4. These are computed in the same way as described above for BHs. In addition to close encounters between MSstars and BHs (NSs), we also record all close encounters of BHs (NSs) with giants and white dwarfs (WDs). For these events, we consider a maximum pericenter distance of for all models, where is the stellar radius of the giant or WD. We discuss BH–giant and BH–WD TDEs further in Section 3 and we examine BH–giant encounters specifically in Section 6.
For stellar pairs where neither object is a BH or a NS (for example, MS – MS) CMC includes only physical collisions, treated in the “stickysphere” approximation, i.e., the two stars are merged assuming conservation of mass and momentum whenever , where and are the radii of the two stars (Chatterjee2013b, for more detail, see). All stellar radii are computed using the evolutionary tracks of SSE (Hurley2000).
In CMC, for each dynamical timestep, , we determine whether a given neighboring pair of single stars should undergo a TDE by evaluating the quantity
(3) 
where is the local number density of stars. CMC uses a globally shared timestep that considers all relevant physical processes (twobody relaxation, close interactions, and stellar evolution; see Fregeau2007). For the simulations presented here, at each timestep, is evaluated for each pair of stars and compared with a random draw from . If for a particular pair of objects, a TDE is assumed to have happened; otherwise, we perform twobody relaxation. See Fregeau2007; Chatterjee2010; Goswami2012; Pattabiraman2013 for further detail.
The final outcome of TDEs involving stellarmass BHs (and NSs) is highly uncertain. Performing detailed calculations of the 3D hydrodynamics associated with these events is beyond the scope of this study (see Fixelle2019 for a first attempt). Here, for simplicity, we assume that in the event of a BH–star (or NS–star) TDE, the star is instantaneously destroyed, and no mass is accreted by the BH (NS). Qualitatively, this is because the accretion of even a small amount of mass onto the BH (NS) is expected to easily release enough energy to completely unbind the stellar debris. This process is discussed in more detail in Section 4.
In addition to TDEs arising from single–single encounters, we also consider close encounters of stars with BHs that occur through binary–single and binary–binary interactions (higher multiples are not included in CMC). In a manner similar to the procedure for evaluating whether pairs of single stars should undergo a TDE or relaxation, CMC evaluates a quantity for each binary–single or binary–binary neighboring pair. In this case, is given by an expression analogous to equation 1, but with and for binary–single and binary–binary pairs, respectively. Here is the binary semimajor axis and are parameters setting the distance of closest approach for an interaction to be considered a strong interaction. We have found that the values capture almost all the relevant energygenerating binary interactions; see Fregeau2007 for a more detailed discussion. As before, is compared to a random draw from for each pair of objects at each dynamical timestep. If , the strong interaction for these objects is integrated directly using the fewbody package (Fregeau2004; Fregeau2007) which has been updated to incorporate radiation reaction for all encounters involving BHs (for more information, see Rodriguez2018b; Rodriguez2018a). During all fewbody integrations, when a pair of stars pass within , a TDE is assumed to have occured, as in the single–single case.
For this analysis, we use the same set of models described in Rodriguez2018b (models 124), plus two additional models (models 2526) that adopt an upper limit for tidal disruptions, as described above. A number of initial cluster parameters are fixed throughout, including the initial King concentration parameter (), the initial binary fraction (), and the stellar initial mass function (IMF; we assume an IMF as in Kroupa2001, with initial masses sampled in the range ). Four initial parameters are varied in this grid of models: the total particle number (, , , and ), the virial radius ( and 2 pc), the metallicity (, 0.05, and 0.25 ), and the galactocentric distance (, 8, and 20 kpc). We assume that the metallicity and galacocentric distance values are correlated, giving us a model grid. Table 2 lists the initial parameters for all models used in this study, as well as the total number of BHs retained at Gyr, the cumulative number of BBH mergers, and cumulative number of BH–MS, BH–giant, and BH–WD TDEs in each model. All together, this grid of models approximates the full distribution of cluster masses, sizes, and metallicities observed in the Milky Way.
3 BH–MS TDEs in cluster models
As shown in Table 2, we identify 898 total BH–MS TDEs in models 124 of this study. Of these, 433, 363, and 102 occur through single–single, binary–single, and binary–binary encounters, respectively. For models 2425 (for which we increase the cross section for TDEs to obtain an upper limit; see Section 2), we identify 176 total TDEs, including 89, 80, and 7 from single–single, binary–single, and binary–binary encounters, respectively. Figure 3 shows all BH–star TDEs identified in models 124. For the top panel, we show the mass of the star involved in each TDE on the vertical axis and we show the TDE time on the horizontal axis. Black points indicate BH–MS TDEs, orange points indicate BH–giant TDEs, and blue points indicate BH–WD TDEs. In total we identify 246 BH–giant TDEs and 6 BH–WD TDEs in models 124. We discuss BH–giant TDEs in more detail in Section 6. Because BH–WD TDEs occur at such a relatively low rate compared to the other stellar types, we do not consider them further in this analysis.
The gray curve in the top panel of Figure 3 marks the turnoff mass as a function of time. As seen in the figure, a handful of BH–MS TDEs involve MS stars with masses lying above the turnoff curve. These MS stars are themselves collision products of MS stars and would be identified as blue stragglers (e.g., Sandage1953; Ferraro1995; Piotto2003; Chatterjee2013b).
The middle panel of Figure 3 shows the distribution of MS masses for all BH–MS TDEs shown in the top panel (black curve) compared to the initial mass function (IMF; gray curve). As the middle panel shows, the mass distribution of MS stars disrupted by BHs follows the IMF closely. The median MS mass for all BH–MS TDEs is .
Finally, we show in the bottom panel of Figure 3 the distribution of BH masses for all BH–MS TDEs. In black, we show all first generation BHs (BHs that have not already undergone a binary BH merger) and in blue we show all second generation BHs (BHs that were formed from binary BH mergers earlier in the cluster evolution). The peak at in the black curve comes from our treatment of pair instability supernovae (see Section 2). The handful of first generation BHs with masses slightly above 40.5 are the result of either stable mass transfer or stellar collisions prior to BH formation (see Rodriguez2018a for further detail). In total, only 6 of the roughly 900 total BH–MS TDEs in models 124 () occur with a second generation BH. The median BH mass for all BH–MS TDEs is .
In Figure 3 we show how the number of BH–MS TDEs per cluster varies with cluster mass and size (which are specified in our models by initial and , respectively). Here filled and open circles denote cluster models with pc and pc, respectively. As the figure shows, the number of TDEs varies proportionally with cluster mass and inversely with cluster size. This is as anticipated. The TDE rate scales with the cluster number density, . For clusters of similar total mass, decreasing the size () will increase the density leading to a higher TDE rate. Likewise, for clusters of similar physical size, increasing the number of particles (total mass) will increase the density and again lead to more TDEs.
Comparing models 25 and 26 to models 7 and 19, respectively (which have the same initial conditions), we see that using the upper limit for TDEs described in Section 2, the number of BH–MS TDEs is increased by a small factor. This is anticipated: the median BH mass of all disruptions is approximately and the median MS mass is approximately , thus the impact parameter for disruptions is increased, on average, by a factor of compared to our conservative lower limit where we assume (see Section 2 for further detail).
4 Electromagnetic Signatures
In this section we present an analytic prediction for the electromagnetic signature associated with BH–MS TDEs. We start with a broadbrush picture of the disruption process and disk formation from the fallback debris (Section 4.1). Afterwards, we adopt a simple model for the evolution of superEddington accretion disk and its wind (Section 4.2). The wind carries the radiation generated in the disk (by viscous dissipation) and releases it at larger radii where photons can diffuse away (Section 4.3). The final results, the lightcurves and temperature evolution, are shown in Figures 4.3 and 4.3. The possibility of jet formation is discussed in Section 4.4.
4.1 Tidal disruption and disk formation
Consider a MS star of radius and mass , BH mass ; the mass ratio is . Adopting the upper limit case discussed in Section 2 where we assume the maximum pericenter distance that will lead to a TDE, , is equal to as given by equation 2, we can write
(4) 
where is the gravitational radius of the BH. For stellarmass BHs with masses up to and MS stars with masses as low as , equation (4) gives (as compared to for MSdisruptions by supermassive BHs). The initial orbit of the star has pericenter distance , where is the dimensionless impact parameter. Here, we only consider the most frequent cases where is close to unity. Since is comparable to the stellar radius , the star’s selfgravity cannot be ignored in the hydrodynamic disruption process. This suppresses the amount of marginally bound debris and hence the latetime fallback rate evolution is steeper than the canonical powerlaw (e.g. 2013ApJ...767...25G; Perets2016).
Instead of modeling the detailed disruption process (see Perets2016), we focus on the broadbrush picture of disk formation, viscous accretion, and the electromagnetic signals. Roughly speaking, the stellar debris has a spread in specific (kinetic potential) energy
(5) 
and specific angular momentum
(6) 
where the parameter roughly corresponds to fluid elements at different locations inside the star ( for the part of the star closest to the BH at disruption). The bound debris corresponds to negative specific energy () and has eccentricity
(7) 
We define a circularization radius corresponding to a circular Keplerian orbit with the angular momentum in equation (6),
(8) 
For stellar mass BHs of mass (roughly the typical BH mass identifed in our cluster models; see Section 3), we see that is of order unity and hence the orbits of the bound debris are not highly eccentric. This is very different from TDEs by supermassive BHs where we have . Therefore, as seen in smoothedparticle hydrodynamics (SPH) calculations (see, e.g., Fixelle2019), we expect a nearly circular accretion disk to form quickly after the most tightly bound mass falls back to the pericenter, and the radius of the disk is given by (we ignore the weak dependence on ).
The viscous accretion timescale at radius is
(9) 
where is the Keplerian angular frequency, is the dimensionless viscosity parameter (1973A&A....24..337S), and we have taken the disk height ratio to be 0.5 (appropriate for a superEddington thick disk and the uncertainty can be absorbed into ). Since the latetime fallback rate drops rapidly on timescale (Perets2016), nearly all the bound debris will fall back within the viscous time and accumulate near . Thus, we take the disk mass to be , and then the peak accretion rate is given by
(10) 
Note that the maximum luminosity of the accretion system is , which can only be achieved under conditions for extremely efficient energy release (e.g. in the form of relativistic jets, see §4.4). Defining the Eddington accretion rate as ( being the Eddington luminosity), we have . Thus, the viscous heating rate exceeds the Eddington luminosity by about five orders of magnitude and hence the disk is indeed geometrically thick near .
4.2 SuperEddington accretion disk
The electromagnetic signals of MS–BH TDEs depend on the physics of superEddington accretion, and the relevant MHD processes coupled with radiative transport are not understood in detail (2011ApJ...733..110B; 2011ApJ...736....2O; 2014MNRAS.441.3177M; 2016MNRAS.456.3929S; 2017arXiv170902845J), see recent reviews by 2013LRR....16....1A; 2014SSRv..183...21B. To construct a simple model, we approximate the disk mass distribution as a “ring” located at the peak radius of the surface density distribution and calculate the time evolution of the bulk properties of the disk. The disk mass and angular momentum evolve as (2008MNRAS.388.1729K)
(11) 
where the viscous time is taken to be (fixing the height ratio ) and the parameter is the ratio between the specific angular momentum of the disk wind and that of the disk at the wind launching point. The disk angular momentum at any given time is given by , so the disk radius evolves as
(12) 
The analytical solution to equations (11) and (12) is
(13) 
where is the characteristic evolution timescale and the initial conditions for the disk evolution are taken as (cf. Section 4.1)
(14) 
The time evolution of the disk accretion rate and radius for a number of choices of and are shown in Figures 4.2 and 4.2.
It is widely believed that a large fraction of the disk mass will be blown away in the form of a wind (1973A&A....24..337S; 1994ApJ...428L..13N; 1999MNRAS.303L...1B; 2004MNRAS.349...68B; 2009MNRAS.400.2070S). We approximate the radiusdependent accretion rate as a powerlaw (1999MNRAS.303L...1B)
(15) 
where the powerlaw index , and the upper limit of is set by energetic requirement and the lower limit corresponds to no wind mass loss. Numerical simulations of adiabatic accretion flows show that the powerlaw index likely lies in the range – (2012ApJ...761..129Y, and references therein). The radius of the innermost stable circular orbit (ISCO) depends on the BH spin . If the wind launched at radius has specific angular momentum (where is the lever arm), then the rate at which angular momentum carried away by the wind is (2014ApJ...784...87S), which means . The angular momentum evolution only affects the latetime behavior of the electromagnetic emission at time . Since we are mainly concerned with the peak luminosity and the peak duration, we take the leverarm factor for simplicity and hence
(16) 
Under the above prescription, the system evolves selfsimilarly at and we obtain the wellknown solution: , and (see equation 13). The prescription stays valid until the accretion rate drops below (or the sphericalization radius drops below , 1979MNRAS.187..237B), and then the disk is expected to undergo thermal instability and collapse into a thin one near . The accretion rate may drop by many orders of magnitude at the transition to thindisk phase (2014ApJ...784...87S). In this paper, we focus on the thickdisk phase where an optically thick wind can generate bright optical emission. The time evolution of disk properties in the thickdisk phase are shown in Figures 4.2 and 4.2.
4.3 Radiation from the superEddington wind
Having the disk dynamics in hand, we calculate the radiationhydrodynamic evolution of the wind. The local viscous heating rate per unit area is given by (1974MNRAS.168..603L)
(17) 
We assume that a fraction of this heating power escapes^{1}^{1}1Note that is not much less than unity. For instance, if the superEddington wind launched at each radius has asymptotic velocity of order the local escape speed, then the wind kinetic power has efficiency . , shared by the wind kinetic power and radiative power. Then the total luminosity of the disk is given by
(18) 
where and the function is given by
(19) 
Except for extreme values of (0 or 1), most mass is loaded near (lowest escape speed) but most accretion power is released near radius (highest escape speed), so we expect internal collisions to occur and the total power will be thermalized near radius under spherical symmetry (2012MNRAS.420.2912B). Thus, the bulk kinetic energy and the radiationdominated thermal energy are roughly in equipartition at .
As the wind expands in radius, nearly all the heat is converted back into bulk kinetic energy due to adiabatic expansion, so the asymptotic wind speed is given by , which means
(20) 
Note that the asymptotic wind speed is independent of the accretion rate and is subrelativistic in most situations (unless ).
The radiation carried by the wind is advected until the photontrapping radius where photons can diffuse away over a dynamical time. At a given time , a wind shell at radius was launched at the retarded time given by
(21) 
where and are the disk radius (inner boundary of the wind, equation 13) and the wind speed (equation 20) at time , respectively. Thus, the time evolution of the photontrapping radius can be estimated by matching the diffusion time with the expansion time,
(22) 
where is the retarded time of the wind shell currently at radius , is the outer boundary of the wind launched at , and is the Thomson scattering optical depth of the wind outside radius ,
(23) 
The integrated wind mass loss rate is
(24) 
and the wind density profile is taken to be
(25) 
where we have taken into account wind propagation by using the retarded time for each shell at radius . At the relevant densities and temperatures, the Rosselandmean opacity is dominated by Thomson scattering (for cosmic abundance). The time evolution of the photontrapping radius is shown in Figures 4.2 and 4.2 for two TDEs with the same BH mass but different stellar masses and , respectively. The qualitative result is that the photontrapping radius increases with time initially (due to wind expansion), reaches a maximum that is many orders of magnitude larger than the disk radius (due to high optical depth), and then decreases rapidly as the wind massloss rate drops at late time.
Once we find the photontrapping radius, the bolometric luminosity of the escaping photons is given by the rate at which radiation are advected across the trapping surface, i.e.
(26) 
where the radiation energy density is given by adiabatic expansion (valid at )
(27) 
and and are evaluated using the disk luminosity , disk radius , and asymptotic wind speed at the retarded time for this trapping shell. The temperature (or average energy) of the escaping photons^{2}^{2}2The thermalization radius is defined where the effective optical depth (1979rpa..book.....R), where is the absorption optical depth. In the case when , the color temperature is lower than that in equation (28) by a factor of , but the bolometric luminosity in equation (26) is unaffected. Modeling the (boundfree and freefree) absorption opacity requires detailed nonLTE radiative transfer calculations and is left for future works. is given by
(28) 
where is the StefanBoltzmann constant.
Figures 4.3 and 4.3 show the bolometric lightcurve, asymptotic wind velocity, and the temperature of the photons escaping from the trapping surface. We find that BHMS TDEs generate bright optical transients on timescales of a few to 100 days with peak luminosities from to . The photospheric velocity may range from 0.01 to 0.1 and decreases with time. The largest uncertainty in our model is the mass flux powerlaw index , and the luminosity decreases towards larger (as a result of stronger mass loss and hence less accretion power). We note that our simple model for the disk evolution does not include the rising segment of the wind power at very early time, which depends on the details of the tidal disruption and disk formation on a timescale of a few hours. Instead, we add a short segment^{3}^{3}3A steeper powerlaw generally leads to a more rapidly rising lightcurve, but the latter is much shallower in that , because of the lag between wind launching and photon escaping. The lag time is longer for larger mass flux index due to a denser wind. of in the first sec (or sec) for the case of (or ), which roughly captures the rapid onset of the gas fallback as shown in Figure 1 of Perets2016. Input from threedimensional SPH calculations is needed to discuss the intricacies of disk formation and evolution in more detail in the context of the calculations presented here.
4.4 Jet formation and gammaray burst
Given the high peak accretion rate (equation 10), if the accretion efficiency can reach that of typical active galactic nuclei (e.g. 1982MNRAS.200..115S), then the peak luminosity of the accretion system is of order with duration –sec. Thus, it has been proposed (Perets2016) that TDEs by stellarmass BHs may be responsible for the population of ultralong duration gammaray bursts (ulGRBs, 2013ApJ...766...30G; 2014ApJ...781...13L). Numerical simulations of superEddington accretion flows in the magneticarrested disk (MAD, with dynamically important magnetic flux and rapid BH spin) regime show that high accretion efficiency can be achieved (2011MNRAS.418L..79T; 2015MNRAS.454L...6M), and these models have been applied to ultraluminous Xray sources (2017MNRAS.469.2997N) and jetted TDEs by supermassive BHs (2014MNRAS.437.2744T; 2019MNRAS.483..565C).
2014ApJ...781...13L estimated the event rate of ulGRBs to be a factor of a few lower than that of classical GRBs, of order at redshift , after correcting for the selection bias that faint, longlived events tend to fall below the trigger threshold of Swift BAT. However, ulGRBs are most likely strongly beamed, so the true rate may be much (a factor of 10–100) higher than the BHMS rate predicted by our GC simulations (see Figure 5). Another challenge to the BHMS TDE scenario is the detection of a bright supernova (SN2011kl) associated with the ultralong GRB 111209A (2015Natur.523..189G). SN2011kl had a peak time of 14 restframe days and peak bolometric luminosity of , and its overall spectral and lightcurve shapes were similar to that of GRBassociated broadline typeIc supernovae (with no evidence of H or He). This is inconsistent with the expected Hrich gas composition of BHMS TDEs, although the superEddington disk wind may generate sufficiently bright optical emission (see the case in Figure 4.3).
Therefore, we conclude that it is unlikely that BH–MS TDEs are responsible for the majority of ulGRBs. Nevertheless, when appropriate conditions (e.g. BH spin and magnetic flux) are met, the launching of relativistic jets in BH–MS TDEs is still possible. In the following section, we compute the comoving and cumulative rates of BH–MS TDEs as a function of redshift, showing that, in principle if indeed a relativistic jet is launched, these events can potentially be detectable at very high redshifts.
5 Event Rates
To calculate the event rate of BH–MS TDEs, we use an integral equation similar to that used in RodriguezLoeb2018. The BH–MS TDE rate as a function of cosmic time is given by
(29) 
where is the comoving rate of star formation in GCs (in units of ) per halo mass at redshift corresponding to a cosmic time , from ElBadry2018. is the rate of BH–MS TDEs at time for for a cluster with an initial mass and virial radius . As was done in RodriguezLoeb2018, we assume the individual cluster rate can be described as
(30) 
where we fit the 5 parameters, , separately for clusters with and . We truncate the rate to zero below 100 Myr, which we find reproduces both the rate and total number of BH–MS TDEs from our CMC models.
is the cluster initial mass function (CIMF), which we assume to be a powerlaw between and , with a possible exponential truncation
(31) 
We consider exponential truncations of (corresponding to a purely power law) and (as suggested by observations of young massive star clusters in the local universe, e.g. PortegiesZwart2010). is the mean initial mass of GCs given an assumed CIMF (this is used to convert the mass formed in GCs into a number of GCs). is the mean initial mass of GCs given the assumed CIMF.
We have introduced a new parameter, , corresponding to the fraction of clusters that do not form an intermediatemass BH (IMBH) by runaway collisions of BHs. It has recently been shown Antonini2018 that for clusters with central escape speeds the probability of forming an IMBH through repeated collisions of BHs goes to 1. Since the presence of an IMBH can suppress the standard three and fourbody encounters that facilitate collisions between BHs and MS stars, we assume that any cluster which forms an IMBH does not contribute to the rate calculated. This is incorporated into equation (29) by assuming that a certain fraction of clusters with a given escape speed contribute zero mergers. We find that an equation of the form
(32) 
with and , provides a good fit to Figure 6 of Antonini2018. We define as the central escape speed from a Plummer model of mass and virial radius (e.g., HeggieHut2003),
(33) 
In practice, equation (32) only truncates the contribution to the rate from the most massive and compact clusters. We find that the local comoving merger rate is decreased by only in the local universe. For comparison, this same correction decreases the BBH merger rate from RodriguezLoeb2018 from to in the local universe. The IMBH correction does not alter the rate when the CIMF is truncated by .
Figure 5 shows the rate as a function of redshift using equation 29. The top panel shows the comoving rate and the bottom panel shows the cumulative rate. The solid and dashed black lines show the BH–MS TDE rate for both the fiducial CIMF (solid) and assuming an exponential truncation at (dashed). In the top panel, we also show in blue (green) the rate if only those models with pc (pc) are considered in the calculation. As anticipated, the predicted rate is higher for the models with pc, simply because these models have relatively higher stellar densities, leading to more TDEs.
As discussed in Section 2, if we consider our upper limit for the TDE cross section, the rate shown in Figure 5 by a factor of roughly 3, based on our typical values of BH and MS masses. Therefore, we can extrapolate an upper limit for the BH–MS TDE rate by simply multiplying the rates of Figure 5 by this factor. In the local universe (), we predict a comoving BH–MS TDE rate of and an extrapolated upper limit rate of . If we consider a CIMF truncated at , these rates decrease to and , respectively. For comparison, RodriguezLoeb2018 predicted a typical BBH merger rate of in the local universe, roughly comparable to our extrapolated upper limit TDE rate if we assume no truncation for the CIMF.
Assuming the density of Milky Way equivalent galaxies (MWEGs) in the local universe is one per (Abadie2010), we predict a BH–MS TDE rate of with an upper limit of roughly in the local universe. Note that the latter rate is in approximate agreement with the analytic predictions made by Perets2016.
6 BH–giant close encounters
In addition to the BH–MS TDEs which are the main focus of this study, we also report all close encounters of BH–giant pairs identified in our models. The orange points in the top panel of Figure 3 show these events. In total, we find 246 giant TDEs in models 124. Although the cross section is higher for giants compared to MS stars (simply due to their relatively large stellar radii) the TDE rate is limited by the relatively small number of giants, which is a consequence of their shorter lifetimes.
Detailed consideration of the outcome of BH–giant close encounters is outside the scope of this study. We simply cite the predicted outcomes from earlier work, and reserve a more detailed study of these events for a future analysis. In particular, Ivanova2010 showed that close encounters of giants and stellarmass BHs may serve as a formation channel for BH–WD binaries, an important result given that observations suggest that a fraction of the stellarmass BH candidates identified in GCs to date may indeed have WD companions (Strader2012; Bahramian2017). For massive clusters (number densities of and velocity dispersions of ), Ivanova2010 predicted BH–WD binaries will form through BH–giant close encounters at a rate of per BH per Gyr, using a calculation of the form:
(34) 
(see, e.g., Ivanova2005; Ivanova2010) where is the fraction of giants in the stellar population within the GC core, is the number density, is the typical BH mass, is the typical stellar radius of giants, and is the velocity dispersion. describes how close a typical BH–giant encounter must be to result in a disruption. In our models, we take , but as discussed in Ivanova2010, as high as may also be appropriate.
In Figure 6, we show rate of formation of BH–WDs through BH–giant close encounters computed for various late time cluster snapshots (Gyr) in models 124 plotted against the total number of BHs retained in the cluster at each respective time. We compute the rate, using equation 34 with , , , , and computed uniquely for each cluster snapshot. As in Figure 3, filled circles denote models with pc and open circles denote models with pc. We limit ourselves to late times here to reflect the typical ages of presentday GCs. BH–WD binaries formed through close encounters of BHs and giants at earlier times are unlikely to survive to the present as a result of the frequent dynamical encounters which may break apart the binaries.
Figure 6 shows similar trends to those shown in Figure 3 for the BH–MS TDEs. In particular, the formation rate of BH–WDs varies directly with (as seen by comparing the different colors) and inversely with cluster size (; as seen by comparing the open versus filled circles), as expected from equation 34. The dependence of the rate on the total number of retained BHs is more complex. The overall trend is a direct relation between and , as expected because more BHs are formed in more massive clusters. However, for models of fixed mass (i.e., points of a single color in Figure 6), it can be seen that varies inversely with . This results from the way BH populations shape the dynamical evolution of their host cluster. As discussed in Kremer2019a, BH populations provide an internal “heating” source for their host cluster. As a result, clusters with many BHs are relatively diffuse while clusters with few BHs may undergo corecollapse leading to relatively high central densities. Because scales with density, this means clusters with smaller BH populations will actually lead to more BH–giant encounters. Furthermore, as the figure shows, this process is coupled with the initial : clusters with pc (filled circles) retain fewer BHs at late times compared to their pc counterparts (open circles). This is also consistent with Kremer2019a, who showed that the initial cluster size is the key parameter for determining the number of BHs retained in a cluster at present.
As Figure 6 shows, close encounters of BHs and giants may lead to the formation of up to approximately one BH–WD binary per Gyr per cluster, possibly sufficient to explain the handful of accreting BH–WD binary candidates observed in local clusters. In addition to the BH–giant encounter channel considered here, BH–WD binaries compact enough to start mass transfer may also form through series of binarymediated exchange encounters, as discussed in Kremer2018a.
Finally, we note that, in addition to being potentially observed as masstransferring lowmass Xray binaries, BH–WD binaries may also be observed as gravitational wave (GW) sources by lowfrequency GW detectors such as LISA (e.g., Kremer2018c). Future electromagnetic and GW observations of these sources will continue to provide better constraints on the formation of these types of binaries.
7 Conclusions and Discussion
7.1 Summary
We have explored the disruption of stars by stellarmass BHs in GCs through close encounters. We summarize our main findings below.

Using our Monte Carlo code CMC to model the evolution of GCs, we show that stellarmass BHs disrupt MS stars, giants, and WDs throughout the lifetime of the cluster. These TDEs can occur through both single–single and binarymediated encounters (binary–single and binary–binary).

The number of TDEs per cluster is determined primarily by the cluster’s initial size (parameterized in our models by the initial virial radius) and initial mass. For our most massive and compact models ( and pc), we get up to 200 BH–MS TDEs over the cluster lifetime. For our lower mass cluster models (), the total number of TDEs can be as low as zero to a few.

By incorporating a realistic cosmological model for GC formation, we derive a rate of BH–MS TDEs of approximately in the local universe and a cosmological rate that peaks at roughly at a redshift of about 3.

We show that the wind mass loss associated with these BH–MS TDEs produces optical transients of luminosity to on timescales of about a day to a month. In Section 4 we show lightcurve predictions expected from these events.

BH–giant close encounters occur at rates of up to per Gyr per cluster. These events may serve as a dynamical formation channel for BH–WD binaries in GCs, which may be observed as Xray or gravitational wave sources.
7.2 Detectability
The Zwicky Transient Facility (ZTF) reaches = 20.8 mag (5) during a single exposure of s, surveying of the sky. Assuming a blackbody temperature of K (see Figures 4.3 and 4.3), the detection horizon in luminosity distance is roughly , where is the bolometric luminosity near the optical peak. The allsky rate is about in the optimistic case where and peak luminosity . In the pessimistic case where and peak luminosity , the allsky rate is about .
7.3 Possible effects on the BH population
When a BH disrupts a star, the disruption may have various effects upon the BH itself, which in turn may alter the overall BH dynamics in the cluster. For instance, if a significant fraction of the mass of the disrupted star were accreted by the BH following the TDE, the BH may be spun up through accretion. Merging BBHs that are highly spinning can get gravitational wave recoil kicks as high as (e.g., Campanelli2007; Lousto2012), significantly larger than the escape speed of a typical GC or even a galactic nucleus. Thus, if a significant number of BHs in a given cluster attain high spins through TDEs of MS stars, all BBH merger products will be ejected from the cluster promptly after merger. This has important implications for the production of secondgeneration mergers in clusters (see, e.g., Rodriguez2018a).
Additionally, as discussed in Section 4, in a typical BH–MS TDE, some fraction of stellar material will become unbound from the system. This unbound mass is ejected in an asymmetric manner, so the BH receives an impulsive kick in response. The unbound debris has positive specific energies in the range , where is given by
(35) 
where we have used . The escape velocity of the star is given by . The maximum speed of the unbound debris is then . The total linear momentum, , carried away by the unbound debris depends on the mass distribution over specific energy, , but a rough estimate is that (this expression is exact for a flat distribution of ). Then the BH receives a kick in the direction of velocity
(36) 
For a typical TDE with a MS star of mass of (with escape velocity of km/s) and a BH of mass of , we obtain km/s. These kicks are low compared to both typical cluster escape velocities () as well as typical dynamical recoil kicks attained from small BH resonant encounters. Thus, these TDE kicks are unlikely to affect the overall BH dynamics in a significant way.
Accounting for these kicks as well as potential BH spinup as a result of these TDEs is beyond the present study. Followup analyses may explore the potential impact of these effects upon the evolution of the BH populations in clusters in more detail.
7.4 Neutron star TDEs
As discussed briefly in Section 2, CMC records close encounters of neutron stars (NSs) and stars in a manner similar to the BH–star encounters. These NS–star interactions may also lead to disruption events, where the NS acts as the disrupting object. Such events were considered in, e.g., Hansen1998 and Perets2016. In particular, Hansen1998 argued that accretion onto the NS during such events may trigger collapse to a BH. In total, we identify 28 NS–MS TDEs events in our models 124. Roughly, this translates to a galactic rate of per Milky Waylike galaxy, two orders of magnitude lower than the rate predicted for BH–MS TDEs. Thus we conclude that NS–MS TDEs do not occur at an astrophysically interesting rate in this set of cluster models.
As discussed in Ye2018, the dynamical interaction rate for NSs in a GC is closely related to the cluster’s BH population. For clusters with BH populations sufficiently large to dynamically affect the cluster through “BH heating,” the NS population will be relegated to the outer parts of the cluster where densities are relatively low. Only when the BH population has been sufficiently depleted will the NSs see a significant boost in their encounter rate, but even then the encounter rate for NSs is limited by the fact that, on average, NSs have masses of the same order as other populations, in particular white dwarfs and MS binaries. As a result, dynamical interactions involving NSs are never as frequent as those involving BHs in typical GCs, explaining the relatively low rate of NS–MS collisions compared to BH–MS TDEs. However, in a corecollapsed cluster with relatively higher stellar densities and BHs (see, e.g., Kremer2019a), NS–MS TDEs may become more frequent.
7.5 Link to process enrichment?
Observational evidence of process enrichment exists in several Milky Way GCs, including M5, M15, M92, and NGC 3201 (e.g., Roederer2011; Bekki2018). The LIGO/Virgo detection of the binary NS merger GW170817 (Abbott2017) and the followup electromagnetic observations showed that NS mergers produce large amounts of process elements (Kasen2017). However, explaining the observed process abundance in GCs with NS merger events proves to be difficult. Recent observational evidence shows that Milky Way GCs exhibit multiple stellar populations which are formed over a series of star formation episodes that can span tens to even hundreds of Myrs (for a recent review on the formation of multiple populations in GCs see, e.g., Gratton2012). For a particular binary NS merger event (or series of events) to enrich a GC’s stars with process material, the event(s) must occur while star formation is still taking place (see, e.g., Bekki2018; Zevin2019 for further discussion).
Merging binary NSs are expected to form in GCs through two mechanisms: binary evolution of primordial binaries (e.g., Ivanova2003; Dominik2012; Tauris2017) and dynamical formation of NSs at late times (e.g., Ye2018). In the first scenario, where the NS components typically form through iron corecollapse supernovae (which are expected to lead to large natal kicks; Hobbs2005), it is not straightforward to produce binary NSs that remain bound to the cluster or have merger times sufficiently low such that the merger takes place within the cluster environment during GC star formation (Zevin2019; also see Safa2018 for related discussion in the context of ultrafaint dwarf galaxies). In the second scenario, where binary NSs are formed dynamically through exchange encounters, one must wait for the NSs to masssegregate to the cluster core, which takes Gyrs, due to the relatively low NS masses (see, e.g., Ye2018; Zevin2019). Thus, by the time binary NSs have begun to form dynamically and subsequently merge, any star formation episodes will have almost certainly ceased. In this case, any process material produced in these latetime NS merger events would be unable to enrich the stellar population.
We propose that if accretion onto a BH arising from the BH–MS disruptions of this study leads to the production of neutronrich material, these events may provide a way to enrich GCs with process elements at early times. As shown in Kremer2018b, the cluster NGC 3201, one of the MW clusters in which process enrichment is observed, is consistent with hosting a large population of BHs at present. In particular, our models 7 and 19 have final masses, metallicities, and BH numbers consistent with both observed and theoretical predictions for NGC 3201. In these models, we identify 28 and 35 BH–MS collisions, respectively. Of these, 6 and 2, occur within the first 100 Myr of cluster evolution. (We adopt 100 Myr as an approximate time duration for star formation episodes. The exact duration may be much shorter or longer. Again, see Gratton2012 for further detail). Of course, more careful consideration must be taken to determine whether a single BH–MS disruption could produce sufficient (or indeed, any) process material to explain the observed enrichment. In particular, it is unclear whether the disk expected to form during such an event will reach high enough densities to produce process elements. We simply note that on the basis of event rates alone the BH–MS disruptions could in principle serve as a possible mechanism. We also note that other mechanisms such as the collapsar model discussed in (Siegel2018) may serve as viable alternatives for process production in GCs. A complete exploration of process production in GCs is beyond the scope of this paper; see Zevin2019 for a more detailed discussion on the topic.
7.6 Directions for future work
We note that there are several complexities associated with BH–star TDEs not captured in this analysis. For example, how do effects of metallicity and/or stellar age (and therefore density profile of the star) alter the outcome of these disruption events? Also, are there differences expected if the disrupted star is itself a collision product (i.e., a blue straggler, as discussed briefly in Section 3)? Furthermore, under what circumstances is a relativistic jet expected to be launched? Additionally, recent work (e.g., Samsing2017; Samsing2018e) has shown that the inclusion of tidal coupling within small dynamical encounters may have a significant effect on tidal disruptions. Many of these complexities will require detailed 3D (magneto)hydrodynamic calculations of these encounters. Our first exploration of some of these complexities will be presented in a forthcoming paper (Fixelle2019).
The distinction between direct, physical collisions where a BH passes within the radius of a star and more distant encounters near the tidal disruption boundary like those considered in Section 4 is likely important when considering the electromagnetic signature of these events. For instance, a direct collision may lead to prompt accretion onto the BH which may make the effects of feedback critical. Here, we have remained agnostic to these differences and simply included the physical collision and more distant tidal disruption limits to bracket the expected range of close encounters between BHs and stars leading to TDEs. More careful treatment of the tidal interactions, especially in the context of binarymediated resonant encounters, is necessary to explore outcomes of these events in greater detail. In a future study, we hope to explore binarymediated close encounters between BHs and stars using a small integrator that selfconsistently incorporates relevant hydrodynamic effects to examine some of these complexities in more detail (for previous work on the topic, see, for example, Goodman1991).
Finally, we note that the number of single–single TDEs identified in this analysis is roughly comparable to the number of TDEs through binary encounters (binary–single or binary–binary; compare columns 8, 9, and 10 in Table 2). This is perhaps surprising, given earlier work showing that the number of close encounters can increase substantially for binarymediated encounters relative to single–single encounters alone (e.g., Bacon1996; Fregeau2004; Chatterjee2013). For example, Chatterjee2013 showed that binary interactions lead to a marked increase in MS–MS collisions, which may lead to the formation of blue stragglers. Several reasons may explain why we do not see a similar result here in the case of BH–MS TDEs. First, as discussed in Kremer2018a, interactions involving BHs and stars are sensitive to the details of the BH–star “mixing zone,” the properties of which are determined by a complex interaction between the BH population and the rest of the cluster. In a manner similar to its effect on the formation of accreting BH–star binaries, this interaction may limit the number of binary–mediated BH–star TDEs. Second, as seen in Table 1 in Chatterjee2013, the models with the most marked increase in binarymediated close encounters in that analysis are those with initial global binary fractions as high as (and as high as roughly in the core). Here we consider models with lower binary fraction (). A more expansive set of cluster models covering a broader range in initial binary fraction are necessary to evaluate the relative contribution of single–single encounters and binary encounters to the overall TDE rate.
We thank Sasha Tchekhovskoy, Pablo Marchant, Johan Samsing, Michael Zevin, and Enrico RamirezRuiz for useful discussions. We also acknowledge discussions with Tony Piro and Eliot Quataert on wind reprocessing. This work was supported by NASA ATP Grant NNX14AP92G and NSF Grant AST1716762. KK acknowledges support by the National Science Foundation Graduate Research Fellowship Program under Grant No. DGE1324585. WL is supported by the David and Ellen Lee Fellowship at Caltech.
mybib