Tidal Disruptions of Stars by Black Hole Remnants in Dense Star Clusters

TL;DR

This paper investigates tidal disruption events caused by stellar-mass black holes in dense star clusters, predicting their rates, observable signatures, and potential contributions to binary formation through advanced simulations and cosmological models.

Contribution

It provides the first detailed rate predictions and observational signatures for TDEs by stellar-mass black holes in globular clusters using state-of-the-art N-body simulations and cosmological models.

Findings

TDE rates in globular clusters are significant and depend on cluster parameters.

Predicted optical transients from these TDEs are detectable by current surveys.

The cosmological rate peaks at roughly 25 Gpc^{-3} yr^{-1} at redshift 3.

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 stellar-mass BHs through collisions and close encounters. Using state-of-the-art $N$-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 - main-sequence-star TDE rate of approximately $3\,\rm{Gpc}^{-3}\,\rm{yr}^{-1}$ in the local universe ($z<0.1$) and a cosmological rate that peaks at roughly $25\,\rm{Gpc}^{-3}\,\rm{yr}^{-1}$ for redshift 3. Furthermore, we show that the ejected mass associated with these TDEs could produce optical transients of luminosity $\sim 10^{41} - 10^{44}\rm\,erg\,s^{-1}$ 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 - white-dwarf binaries.

Figures

Tables

Table 1. Table 1: Initial and final cluster properties for all models

Model	$N$	Z	$r_v$	$M_{\mathrm{tot}}$	$N_{\mathrm{BH}}$	$N_{\mathrm{BBH},\mathrm{merger}}$	BH–MS TDEs				BH–Giant TDEs	BH–WD TDEs
	($\times {10}^5$)	($Z_{\odot }$)	($\mathrm{pc}$)	($\times {10}^5{M}_{\odot }$)			s–s	b–s	b–b	total
Total
1	2	0.01	1	0.5	1	7	0	2	0	2	1	0
2	5	0.01	1	1.46	41	42	5	6	2	13	3	0
3	10	0.01	1	3.09	165	114	23	12	4	39	9	0
4	20	0.01	1	6.3	634	341	57	49	7	113	29	2
5	2	0.25	1	0.43	1	9	1	2	3	6	6	1
6	5	0.25	1	1.35	36	48	4	9	3	16	6	0
7	10	0.25	1	2.84	198	121	17	9	2	28	15	0
8	20	0.25	1	5.94	608	344	59	53	13	125	29	0
9	2	0.05	1	-	-	7	3	0	1	4	3	0
10	5	0.05	1	0.87	9	32	7	6	6	19	12	0
11	10	0.05	1	2.45	89	102	18	23	8	49	18	1
12	20	0.05	1	5.29	423	331	93	71	10	174	57	1
13	2	0.01	2	0.54	32	4	0	1	0	1	0	0
14	5	0.01	2	1.54	114	27	3	2	0	5	3	0
15	10	0.01	2	3.19	416	66	9	8	0	17	2	0
16	20	0.01	2	6.61	1192	196	26	21	5	52	9	0
17	2	0.25	2	0.47	23	10	0	0	0	0	0	0
18	5	0.25	2	1.4	136	20	3	2	1	6	1	0
19	10	0.25	2	2.98	415	75	18	5	2	25	3	0
20	20	0.25	2	6.3	1133	203	30	11	1	42	8	0
21	2	0.05	2	-	-	8	1	1	0	2	1	0
22	5	0.05	2	0.8	55	24	3	7	0	10	0	0
23	10	0.05	2	2.48	260	59	8	10	5	23	2	0
24	20	0.05	2	5.78	871	177	29	23	2	54	12	0
${25}^{\star }$	10	0.05	1	2.82	171	111	67	56	6	129	18	0
${26}^{\star }$	10	0.05	2	2.97	375	78	22	24	1	47	6	0

Equations

$$\Sigma_{\rm{enc}}(r_{p})=\pi r_{p}^{2}\Big{[}1+\frac{2G(M_{1}+M_{2})}{r_{p}\sigma_{v}^{2}}\Big{]}$$

$$r_{\rm{T}}=f_{p}\Big{(}\frac{M_{\rm{BH}}}{M_{\star}}\Big{)}^{1/3}\,R_{\star},$$

$$r_{\mathrm{TDE}}=R_{\ast }{q}^{-1/3}\simeq (1.0\times 10^5{r}_{\mathrm{g}})M_1^{-2/3}{m}_{\ast }^{-1/3}{r}_{\ast },$$

$$\mathcal{E}\simeq \frac{GM}{r_{\mathrm{TDE}}}-\frac{GM}{r_{\mathrm{TDE}}-\eta R_{\ast }}$$

$$\ell \simeq (r_{\mathrm{TDE}}-\eta R_{\ast })\sqrt{2GM/r_{\mathrm{TDE}}},$$

$$e=\sqrt{1+2\mathcal{E}{\ell }^2/r_{\mathrm{g}}^2{c}^4}\simeq \sqrt{1-4\eta q^{1/3}(1-\eta q^{1/3})}.$$

$$r_{\mathrm{c}}={\ell }^2/r_{\mathrm{g}}{c}^2\simeq 2r_{\mathrm{TDE}}(1-\eta q^{1/3}{)}^2.$$

$$t_{\mathrm{vis}}\simeq \frac{(r_{\mathrm{c}}/H{)}^2}{\alpha {\Omega }_{\mathrm{K}}}\sim (0.74\mathrm{d})({\alpha }_{-1}{)}^{-1}{m}_{\ast }^{-1/2}{r}_{\ast }^{3/2},$$

$${\dot{M}}_{\mathrm{peak}}\sim \frac{M_{\ast }/2}{t_{\mathrm{vis}}}\sim (2.5\times 10^2{\mathrm{M}}_{\odot }\mathrm{y}{\mathrm{r}}^{-1})\frac{{\alpha }_{-1}{m}_{\ast }^{3/2}}{r_{\ast }^{3/2}}.$$

$${\dot{M}}_{\mathrm{d}}(r_{\mathrm{d}})=-M_{\mathrm{d}}/t_{\mathrm{vis}},{\dot{J}}_{\mathrm{d}}=-f_{\mathrm{j}}{J}_{\mathrm{d}}/t_{\mathrm{vis}},$$

$${\dot{r}}_{\mathrm{d}}=2(1-f_{\mathrm{j}})r_{\mathrm{d}}/t_{\mathrm{vis}}.$$

$$\begin{array}{rl}{r}_{\mathrm{d}}& =r_{\mathrm{d},0}{\left(1+\frac{t}{t_0}\right)}^{\frac{2}{3}},\\ {\dot{M}}_{\mathrm{d}}& =-\frac{M_{\mathrm{d},0}}{3(1-f_{\mathrm{j}})t_0}{\left(1+\frac{t}{t_0}\right)}^{-\frac{4-3f_{\mathrm{j}}}{3(1-f_{\mathrm{j}})}},\end{array}$$

$$r_{\mathrm{d},0}=r_{\mathrm{TDE}},M_{\mathrm{d},0}=M_{\ast }/2.$$

$$\dot{M}(r)=\dot{M}(r_{\mathrm{d}})(r/r_{\mathrm{d}}{)}^s,r_{\mathrm{isco}}<r<r_{\mathrm{d}},$$

$$f_{\mathrm{j}}\sim s/(s+1/2).$$

$$Q(r)=\frac{3GM\dot{M}(r)}{4\pi r^3}\left[1-{\left(\frac{r_{\mathrm{isco}}}{r}\right)}^{1/2}\right].$$

$$L_{\mathrm{d}}=f_{ϵ}{\int }_{r_{\mathrm{isco}}}^{r_{\mathrm{d}}}Q(r)2\pi r\mathrm{d}r=\frac{3f_{ϵ}{r}_{\mathrm{g}}}{2r_{\mathrm{isco}}}\dot{M}(r_{\mathrm{d}})c^{2}g(x,s),$$

$$g(x,s)=x^{-s}\left(\frac{1-x^{s-1}}{1-s}-\frac{1-x^{s-1.5}}{1.5-s}\right).$$

$$v_{\mathrm{w}}/c={\left[\frac{3f_{ϵ}{r}_{\mathrm{g}}g(x,s)}{r_{\mathrm{isco}}(1-x^{-s})}\right]}^{1/2}.$$

$$r=r_{\mathrm{d}}(t_0)+v_{\mathrm{w}}(t_0)(t-t_0),$$

$$\frac{\tau (r_{\mathrm{tr}})}{c}\frac{r_{\mathrm{tr}}(r_{\mathrm{out}}-r_{\mathrm{tr}})}{r_{\mathrm{out}}}\simeq t-t_0(r_{\mathrm{tr}}),$$

$$\tau (r)={\int }_r^{r_{\mathrm{out}}}{\kappa }_{\mathrm{s}}\rho (r)\mathrm{d}r.$$

$${\dot{M}}_{\mathrm{w}}=\dot{M}(r_{\mathrm{d}})-\dot{M}(r_{\mathrm{isco}}),$$

$$\rho (r)={\dot{M}}_{\mathrm{w}}(t_0)/[4\pi r^2{v}_{\mathrm{w}}(t_0)],$$

$$L_{\mathrm{bol}}(t)=4\pi r_{\mathrm{tr}}^{2}U(r_{\mathrm{tr}})\left[v_{\mathrm{w}}(t_0)-\frac{\mathrm{d}{r}_{\mathrm{tr}}}{\mathrm{d}t}\right],$$

$$U(r_{\mathrm{tr}})=U(r_{\mathrm{d}}){\left[\frac{\rho (r_{\mathrm{tr}})}{\rho (r_{\mathrm{d}})}\right]}^{4/3},$$

$$T(r_{\mathrm{tr}})={\left[U(r_{\mathrm{tr}})c/4{\sigma }_{\mathrm{SB}}\right]}^{1/4},$$

$$\mathcal{R}(t)=\iiint\frac{\dot{M}_{\rm{GC}}}{d\log M_{\rm{Halo}}}\bigg{\rvert}_{z(\tau)}\frac{1}{\langle M_{\rm{GC}}\rangle}P(M_{\rm{GC}})\\ \times R(r_{v},M_{\rm{GC}},\tau-t)f_{\rm noIMBH}(r_{v},M_{\rm GC})dM_{\rm{Halo}}dM_{\rm{GC}}d\tau$$

$$R(M_{\mathrm{GC}},t)\equiv (A{M}_{\mathrm{GC}}^2+B{M}_{\mathrm{GC}}+C)\times t^{-(\gamma +{\gamma }_M{\log }_{10}{M}_{\mathrm{GC}})}$$

$$\varphi (M_{\mathrm{GC}})d{M}_{\mathrm{GC}}\propto M_{\mathrm{GC}}^{-2}\exp (-M_{\mathrm{GC}}/M_{\mathrm{GC}}^{★})d{M}_{\mathrm{GC}}.$$

Full text

globular clusters: general–stars: black holes–stars: kinematics and dynamics–methods: numerical

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., Sandage, 1953; Ferraro et al., 1995; Piotto et al., 2002; Geller & Mathieu, 2011; Leigh et al., 2013; Chatterjee et al., 2013), millisecond pulsars (e.g., Lyne et al., 1987; Sigurdsson & Phinney, 1995; Ransom, 2008; Ivanova et al., 2008; Fragione et al., 2018; Ye et al., 2019), and low-mass X-ray binaries (e.g., Clark, 1975; Verbunt et al., 1984; Heinke et al., 2005; Ivanova, 2013; Giesler et al., 2018; Kremer et al., 2018b). The influence of dynamics is particularly pronounced for stellar-mass black holes (BHs), which preferentially reside in the highest-density 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 binary-mediated dynamical encounters (e.g., Moody & Sigurdsson, 2009; Banerjee et al., 2010; Bae et al., 2014; Ziosi et al., 2014; Rodriguez et al., 2015, 2016; Askar et al., 2017; Banerjee, 2017; Giesler et al., 2018; Hong et al., 2018; Fragione & Kocsis, 2018; Samsing & D’Orazio, 2018; Rodriguez et al., 2018a). Thus, as LIGO and Virgo have begun to unveil the gravitational-wave universe through detections of merging BBHs (e.g., Abbott et al., 2016a, b, c, d, 2017; The LIGO Scientific Collaboration et al., 2018), the topic of binary BH formation in GCs has risen in importance.

In the past two decades, our understanding of stellar-mass BH populations in GCs has evolved significantly. From an observational perspective, several BH candidates have been identified in GCs through both X-ray/radio (e.g., Maccarone et al., 2007; Strader et al., 2012) and radial-velocity (Giesers et al., 2018) 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 ($t\gtrsim 10$ Gyr) and the BHs readily mix with luminous stellar populations forming BH – non-BH binaries. Furthermore, it is now understood that the BHs can play a dominant role in determining the large-scale structural properties of their host cluster (e.g., Merritt et al., 2004; Mackey et al., 2007, 2008; Peuten et al., 2016; Wang et al., 2016; Arca Sedda et al., 2018; Kremer et al., 2018c; Zocchi et al., 2019; Kremer et al., 2019).

In a GC, stars will frequently pass through small-$N$ (most commonly, three- or four-body) resonant encounters that lead to a greatly increased rate of close passages between objects relative to direct single–single interactions (e.g., Bacon et al., 1996; Fregeau et al., 2004). If these small-$N$ 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, Ivanova et al. (2005) showed that direct physical collisions between neutron stars (NSs) and giant stars can lead to bright, ultracompact X-ray binaries, similar to those observed in Milky Way GCs, in which the NS is accreting from a white-dwarf (WD) companion. Similarly, Ivanova et al. (2010) showed that collisions between BHs and giants may result in accreting BH – WD binaries, possibly similar to the X-ray/radio source observed in NGC 4472 (Maccarone et al., 2007). More recently, Ivanova et al. (2017) showed that (grazing) tidal captures of subgiants by stellar-mass BHs may lead to BH – sub-subgiant binaries, similar to the BH candidate observed in M10 (Shishkovsky et al., 2018).

Additionally, TDEs involving stellar-mass BHs may produce interesting electromagnetic transient signatures. Perets et al. (2016) showed that disruptions of stars by stellar-mass BHs and NSs, which they called “micro-TDEs,” may give rise to bright, energetic, long X-ray/gamma-ray flares, which could possibly resemble ultra-long gamma-ray bursts. Recently, Lopez et al. (2018) 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. Samsing et al. (2019) 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 gravitational-wave signals. Additionally, Fragione et al. (2019) recently showed that secular effects relevant in hierarchical triple systems may also lead to BH–star TDEs that occur at significant rates.

In this analysis, we use our Hénon-type Monte Carlo code CMC to explore the rates of TDEs by stellar-mass 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 long-term evolution of GCs, including our treatment of BH TDEs. In Section 3 we discuss the total number of BH – main-sequence-star (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

CMC (for Cluster Monte Carlo) is a rigorously tested, Hénon-type Monte Carlo code that computes the long-term evolution of GCs with realistic initial properties and containing up to several million stars (Henon, 1971; Hénon, 1971; Joshi et al., 2000, 2001; Fregeau et al., 2003; Pattabiraman et al., 2013; Chatterjee et al., 2010, 2013; Rodriguez et al., 2015). CMC incorporates stellar and binary evolution using the evolution codes SSE and BSE (Hurley et al., 2000, 2002), modernized to reflect our most up-do-date understanding of compact object formation, including prescriptions for natal kicks, mass-dependent fallback, and (pulsational) pair-instability supernovae (Fryer & Kalogera, 2001; Belczynski et al., 2002; Hobbs et al., 2005; Morscher et al., 2015). We also integrate all small-$N$ resonant dynamical encounters (including binary–single and binary–binary encounters) using the fewbody package (Fregeau et al., 2004; Fregeau & Rasio, 2007) which has been updated to incorporate gravitational radiation reaction for all encounters involving BHs (for more information, see Rodriguez et al., 2018a, b).

We use the same prescription to treat pulsational-pair instabilities and pair-instability supernovae as Rodriguez et al. (2018b), which follows Belczynski et al. (2016). Briefly, we assume that any star with a pre-explosion helium core between 45 and 65 $M_{\odot}$ will undergo pulsations that eject a significant amount of mass, until the final product is at most 45 $M_{\odot}$. Assuming that $10\%$ of the mass is lost during the conversion from baryonic to gravitational mass at the time of collapse, this leaves a BH of mass $40.5\,M_{\odot}$. Finally, we assume that stars with helium core masses in excess of $65\,M_{\odot}$ are completely destroyed by a pair-instability supernova (no remnant is formed). Of course, there are many uncertainties associated with these numbers. In particular, some studies, (e.g., Woosley, 2016; Spera & Mapelli, 2017) have predicted that the lower boundary of this mass gap may be as high as $50\,M_{\odot}$.

For NS formation and evolution, we have updated the SSE and BSE codes to include changes to the magnetic field and spin-period evolution, as well as the natal kick prescriptions for NSs formed in ECSNe (Kiel et al., 2008; Kiel & Hurley, 2009). See Ye et al. (2019) for further detail on the treatment of NSs and radio pulsars in CMC.

In a GC, the stellar dynamics in the high-density core can frequently lead to close encounters between stars. For encounters of two single stars with masses $M_{1}$ and $M_{2}$, the cross-section for an encounter with pericenter distance $r_{p}$ is given by

[TABLE]

where $\sigma_{v}$ is the local velocity dispersion (see, e.g. Fregeau & Rasio, 2007). 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 Fabian et al. (1975), we define the characteristic radius for tidal interaction as

[TABLE]

where $M_{\rm{BH}}$ is the BH mass, $M_{\star}$ and $R_{\star}$ are the stellar mass and radius, respectively, and $f_{p}$ is a dimensionless parameter that captures the details of the tidal interaction. In order of increasing $r_{p}$, close encounters may take the form of physical collisions (e.g., Fryer & Woosley, 1998; Hansen & Murali, 1998; Zhang & Fryer, 2001), tidal disruptions (e.g., Perets et al., 2016), tidal captures (e.g., Fabian et al., 1975; Ivanova et al., 2017), or more distant tidal encounters that only weakly perturb the star (e.g., Alexander & Kumar, 2001). 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 low-mass 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, and is beyond the computational scope of CMC. Therefore, we simply consider two limits for our choice of the maximum pericenter distance, $r_{\rm{TDE}}$, that leads to a TDE. As a conservative lower limit, we take $r_{\rm{TDE}}=R_{\star}$ (models 1-24 in Table 1 adopt this assumption) and as an upper limit, we assume $r_{\rm{TDE}}=r_{\rm{T}}$ computed by equation 2 with $f_{p}=1$ (models 25-26 adopt this assumption).111The choice of $f_{p}=1$ as our upper limit is meant to simply reflect the classical tidal disruption limit assuming a star of uniform density. While this approximation is reasonable for low-mass MS stars, a different choice for $f_{p}$ may be more appropriate for evolved massive stars. As this is a first attempt at studying these TDEs in our cluster models, we adopt $f_{p}=1$ for all stars for simplicity and reserve a more detailed exploration of the dependence of $f_{p}$ on stellar type for a later study. We do emphasize however, that a stellar-type-dependent $f_{p}$ is unlikely to change the results presented here significantly. Henceforth, we use the term TDE as shorthand to denote all close encounters with $r_{p}<r_{\rm{TDE}}$, computed in one of these two limits.

Although BH – MS-star TDEs are the main focus of this study, we comment briefly on NS – MS-star TDEs in Section 7.4. These are computed in the same way as described above for BHs. In addition to close encounters between MS-stars 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 $r_{\rm{TDE}}=R_{\star}$ for all models, where $R_{\star}$ 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 “sticky-sphere” approximation, i.e., the two stars are merged assuming conservation of mass and momentum whenever $r_{p}<R_{1}+R_{2}$, where $R_{1}$ and $R_{2}$ are the radii of the two stars (Chatterjee et al., 2013, for more detail, see). All stellar radii are computed using the evolutionary tracks of SSE (Hurley et al., 2000).

Close encounters leading to TDEs can occur through both single–single dynamical encounters as well as binary-mediated encounters (binary–single and binary–binary encounters; higher multiples are not included in CMC). For a discussion of how TDEs are computed in CMC, we direct the reader to Fregeau & Rasio (2007), which describes the elements of the calculation for both the single–single case and the binary-mediated case in detail. For a more general discussion of the calculation of TDEs and collisions within a Monte Carlo-type dynamics code, we direct the reader to Chatterjee et al. (2010); Goswami et al. (2012); Pattabiraman et al. (2013) as well as earlier work by Freitag & Benz (2002) and Giersz (2001).

The final outcome of TDEs involving stellar-mass BHs (and NSs) is highly uncertain. Performing detailed calculations of the 3-D hydrodynamics associated with these events is beyond the scope of this study (see Fixelle et al. (2019) 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. However, the possibility that a significant fraction of the disrupted star could be accreted by the BH cannot be entirely ruled out. If indeed a significant portion is accreted, it could affect the cluster’s BH population, as discussed further in Section 7.3.

For this analysis, we use the same set of models described in Rodriguez et al. (2018a) (models 1-24), plus two additional models (models 25-26) 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 ($W_{0}=5$) and the stellar initial mass function (IMF; we assume an IMF as in Kroupa (2001) with initial masses sampled in the range $0.1-150\,M_{\odot}$). We assume an initial binary fraction of $10\%$ in all models with binary mass ratios drawn from a flat distribution, binary semi-major axes are drawn from a distribution flat in log from near contact to the hard/soft boundary, and initial eccentricities drawn from a thermal distribution. Four initial parameters are varied in this grid of models: the total particle number ($N=2\times 10^{5}$, $N=5\times 10^{5}$, $N=10^{6}$, and $N=2\times 10^{6}$), the virial radius ($r_{v}=1$ and 2 pc), the metallicity ($\rm{Z}=0.01$, 0.05, and 0.25 $\rm{Z}_{\odot}$), and the galactocentric distance ($r_{\rm{gc}}=2$, 8, and 20 kpc). We assume that the metallicity and galacocentric distance values are correlated, giving us a $4\times 3\times 2$ model grid. Table 1 lists the initial parameters for all models used in this study, as well as the total number of BHs retained at $t=12\,$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 1, we identify 898 total BH–MS TDEs in models 1-24 of this study. Of these, 433, 363, and 102 occur through single–single, binary–single, and binary–binary encounters, respectively. For models 24-25 (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 1 shows all BH–star TDEs identified in models 1-24. 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 1-24. 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 1 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., Sandage, 1953; Ferraro et al., 1995; Piotto et al., 2002; Chatterjee et al., 2013).

The middle panel of Figure 1 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 $0.5\,M_{\odot}$.

We show in the bottom panel of Figure 1 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 $M_{\rm{BH}}=40.5\,M_{\odot}$ 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 $M_{\odot}$ are the result of either stable mass transfer or stellar collisions prior to BH formation (see Rodriguez et al. (2018b) for further detail). In total, only 6 of the roughly 900 total BH–MS TDEs in models 1-24 ($\approx 1\%$) occur with a second generation BH. The median BH mass for all BH–MS TDEs is $18\,M_{\odot}$.

In Figure 2 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 $N$ and $r_{v}$, respectively). Here, filled and open circles denote cluster models with $r_{v}=1\,$pc and $r_{v}=2\,$ 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, $n\sim N/r_{v}^{3}$. For clusters of similar total $N$, decreasing the size ($r_{v}$) 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.

For all initial models in this study, we adopt a King concentration parameter of $W_{0}=5$ and vary the initial density simply using the initial virial radius. However, in principle, higher values of $W_{0}$ may also yield an increase in the TDE rate and stellar collision rate. Along these lines, several recent papers (e.g., Portegies Zwart et al., 2004; Freitag et al., 2006; Giersz et al., 2015; Mapelli, 2016) have shown that large initial concentrations may even lead to runaway mergers of MS stars in clusters, a possible channel for forming intermediate-mass BHs (IMBHs). An exploration of the effect of initial cluster concentration on the TDE rate is beyond the scope of the present study but is a topic worth consideration in future analyses. We return to this topic of IMBHs briefly in Section 7.6.

Finally, 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 factor of a few (roughly 5 for model 25 versus model 7 and roughly 2 for model 26 versus model 19). This is anticipated: the median BH mass of all disruptions is approximately $18\,M_{\odot}$ and the median MS mass is approximately $0.5\,M_{\odot}$, thus the impact parameter for disruptions is increased, on average, by a factor of $(M_{\rm{BH}}/M_{\star})^{1/3}\approx 3$ compared to our conservative lower limit where we assume $r_{\rm TDE}=R_{\star}$ (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 broad-brush 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 super-Eddington 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 5 and 6. The possibility of jet formation is discussed in Section 4.4.

4.1 Tidal disruption and disk formation

Consider a MS star of radius $R_{*}=r_{*}R_{\odot}$ and mass $M_{*}=m_{*}M_{\rm\odot}$, BH mass $M=10M_{1}M_{\rm\odot}$; the mass ratio is $q=M_{*}/M$. Adopting the upper limit case discussed in Section 2 where we assume the maximum pericenter distance that will lead to a TDE, $r_{\rm{TDE}}$, is equal to $r_{\rm{T}}$ as given by equation 2, we can write

[TABLE]

where $r_{\rm g}=GM/c^{2}=15M_{1}\rm\,km$ is the gravitational radius of the BH. For stellar-mass BHs with masses up to $\sim 50M_{\rm\odot}$ and MS stars with masses as low as $\sim 0.1M_{\rm\odot}$, equation (3) gives $r_{\rm{TDE}}/R_{*}\sim q^{-1/3}\lesssim 8$ (as compared to $r_{\rm{TDE}}/R_{*}\gg 1$ for MS-disruptions by supermassive BHs). The initial orbit of the star has pericenter distance $r_{\rm p}=r_{\rm{TDE}}/\beta$, where $\beta\gtrsim 1$ is the dimensionless impact parameter. Here, we only consider the most frequent cases where $\beta$ is close to unity. Since $r_{\rm{TDE}}$ is comparable to the stellar radius , the star’s self-gravity cannot be ignored in the hydrodynamic disruption process. This suppresses the amount of marginally bound debris and hence the late-time fallback rate evolution is steeper than the canonical power-law $t^{-5/3}$ (e.g. Guillochon & Ramirez-Ruiz, 2013; Perets et al., 2016).

Instead of modeling the detailed disruption process (see Perets et al., 2016), we focus on the broad-brush picture of disk formation, viscous accretion, and the electromagnetic signals. Roughly speaking, the stellar debris has a spread in specific (kinetic $+$ potential) energy

[TABLE]

and specific angular momentum

[TABLE]

where the parameter $\eta\in[-1,1]$ roughly corresponds to fluid elements at different locations inside the star ($\eta=1$ for the part of the star closest to the BH at disruption). The bound debris corresponds to negative specific energy ($0<\eta<1$) and has eccentricity

[TABLE]

We define a circularization radius $r_{\rm c}$ corresponding to a circular Keplerian orbit with the angular momentum in equation (5),

[TABLE]

For stellar mass BHs of mass $M\sim 20M_{\rm\odot}$ (roughly the typical BH mass identifed in our cluster models; see Section 3), we see that $1-e$ 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 $1-e\sim(M_{*}/M)^{1/3}\ll 1$. Therefore, as seen in smoothed-particle hydrodynamics (SPH) calculations (see, e.g., Fixelle et al., 2019), 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 $r_{\rm c}\sim r_{\rm TDE}$ (we ignore the weak dependence on $\eta$).

The viscous accretion timescale at radius $r\sim r_{\rm c}$ is

[TABLE]

where $\Omega_{\rm K}(r_{\rm c})=\sqrt{GM/r_{\rm c}^{3}}$ is the Keplerian angular frequency, $\alpha\sim 0.1\alpha_{-1}$ is the dimensionless viscosity parameter (Shakura & Sunyaev, 1973), and we have taken the disk height ratio $H/r_{\rm c}$ to be 0.5 (appropriate for a super-Eddington thick disk and the uncertainty can be absorbed into $\alpha$). Since the late-time fallback rate drops rapidly on timescale $\Omega_{\rm K}^{-1}$ (Perets et al., 2016), nearly all the bound debris will fall back within the viscous time $t_{\rm vis}(\gg\Omega_{\rm K}^{-1})$ and accumulate near $r\sim r_{\rm c}$. Thus, we take the disk mass to be $\sim M_{*}/2$, and then the peak accretion rate is given by

[TABLE]

Note that the maximum luminosity of the accretion system is $\sim 0.1\dot{M}_{\rm peak}c^{2}\sim 1.4\times 10^{48}(\alpha/0.1)\mathrm{\,erg\,s^{-1}}$, 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 $\dot{M}_{\rm Edd}\equiv L_{\rm Edd}/c^{2}=2.6\times 10^{-8}M_{1}\mathrm{\,M_{\rm\odot}\,yr^{-1}}$ ($L_{\rm Edd}$ being the Eddington luminosity), we have $\dot{M}_{\rm peak}r_{\rm g}/\dot{M}_{\rm Edd}r_{\rm c}\sim 10^{5}\alpha_{-1}\gg 1$. Thus, the viscous heating rate $\sim GM\dot{M}_{\rm peak}/r_{\rm c}$ exceeds the Eddington luminosity by about five orders of magnitude and hence the disk is indeed geometrically thick near $r\sim r_{\rm c}$.

4.2 Super-Eddington accretion disk

The electromagnetic signals of MS–BH TDEs depend on the physics of super-Eddington accretion, and the relevant MHD processes coupled with radiative transport are not understood in detail (Blaes et al., 2011; Ohsuga & Mineshige, 2011; McKinney et al., 2014; Sa̧dowski & Narayan, 2016; Jiang et al., 2017), see recent reviews by Abramowicz & Fragile (2013); Blaes (2014). To construct a simple model, we approximate the disk mass distribution as a “ring” located at the peak radius of the surface density distribution $r_{\rm d}$ and calculate the time evolution of the bulk properties of the disk. The disk mass $M_{\rm d}$ and angular momentum $J_{\rm d}$ evolve as (Kumar et al., 2008)

[TABLE]

where the viscous time is taken to be $t_{\rm vis}(r_{\rm d})=4\alpha^{-1}\Omega_{\rm K}^{-1}(r_{\rm d})$ (fixing the height ratio $H/r_{\rm d}=0.5$) and the parameter $f_{\rm j}$ 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 $J_{\rm d}=(GMr_{\rm d})^{1/2}M_{\rm d}$, so the disk radius evolves as

[TABLE]

The analytical solution to equations (10) and (11) is

[TABLE]

where $t_{0}\equiv[4/3(1-f_{\rm j})]\alpha^{-1}\Omega_{\rm k}^{-1}(r_{\rm d,0})$ is the characteristic evolution timescale and the initial conditions for the disk evolution are taken as (cf. Section 4.1)

[TABLE]

The time evolution of the disk accretion rate and radius for a number of choices of $f_{\rm j}$ and $\alpha$ are shown in Figures 3 and 4.

It is widely believed that a large fraction of the disk mass will be blown away in the form of a wind (Shakura & Sunyaev, 1973; Narayan & Yi, 1994; Blandford & Begelman, 1999, 2004; Strubbe & Quataert, 2009). We approximate the radius-dependent accretion rate as a power-law (Blandford & Begelman, 1999)

[TABLE]

where the power-law index $0\leq s\leq 1$, and the upper limit of $s$ is set by energetic requirement and the lower limit corresponds to no wind mass loss. Numerical simulations of adiabatic accretion flows show that the power-law index likely lies in the range $0.4$–$0.8$ (Yuan et al., 2012, and references therein). The radius of the inner-most stable circular orbit (ISCO) depends on the BH spin $r_{\rm isco}\in[1,9]r_{\rm g}$. If the wind launched at radius $r$ has specific angular momentum $\widetilde{f}_{\rm j}(GMr)^{1/2}$ (where $\widetilde{f}_{\rm j}$ is the lever arm), then the rate at which angular momentum carried away by the wind is $\dot{J}_{\rm d}=-[\widetilde{f}_{\rm j}s/(s+1/2)]J_{\rm d}/t_{\rm vis}$ (Shen & Matzner, 2014), which means $f_{\rm j}=\widetilde{f}_{\rm j}s/(s+1/2)$. The angular momentum evolution only affects the late-time behavior of the electromagnetic emission at time $t\gg t_{\rm vis}(r_{\rm d,0})$. Since we are mainly concerned with the peak luminosity and the peak duration, we take the lever-arm factor $\widetilde{f}_{\rm j}\sim 1$ for simplicity and hence

[TABLE]

Under the above prescription, the system evolves self-similarly at $t\gg t_{\rm vis}(r_{\rm d,0})$ and we obtain the well-known solution: $r_{\rm d}\propto t^{2/3}$, $\dot{M}_{\rm d}\propto t^{-(2s+4)/3}$ and $\dot{M}_{\rm isco}\propto t^{-4(s+1)/3}$ (see equation 12). The prescription stays valid until the accretion rate $\dot{M}_{\rm d}$ drops below $\sim(r_{\rm d}/r_{\rm g})L_{\rm Edd}/c^{2}$ (or the sphericalization radius $\dot{M}_{\rm d}c^{2}r_{\rm g}/L_{\rm Edd}$ drops below $r_{\rm d}$, Begelman, 1979), and then the disk is expected to undergo thermal instability and collapse into a thin one near $r_{\rm d}$. The accretion rate may drop by many orders of magnitude at the transition to thin-disk phase (Shen & Matzner, 2014). In this paper, we focus on the thick-disk phase where an optically thick wind can generate bright optical emission. The time evolution of disk properties in the thick-disk phase are shown in Figures 3 and 4.

4.3 Radiation from the super-Eddington wind

Having the disk dynamics in hand, we calculate the radiation-hydrodynamic evolution of the wind. The local viscous heating rate per unit area is given by (Lynden-Bell & Pringle, 1974)

[TABLE]

We assume that a fraction $f_{\epsilon}$ of this heating power escapes222Note that $f_{\epsilon}$ is not much less than unity. For instance, if the super-Eddington wind launched at each radius has asymptotic velocity of order the local escape speed, then the wind kinetic power has efficiency $\sim[s\dot{M}(r)/r][GM/r]/[2\pi rQ(r)]=(2s/3)[1-\sqrt{r_{\rm isco}/r}]^{-1}$. , shared by the wind kinetic power and radiative power. Then the total luminosity of the disk is given by

[TABLE]

where $x\equiv r_{\rm d}/r_{\rm isco}$ and the function $g$ is given by

[TABLE]

Except for extreme values of $s$ (0 or 1), most mass is loaded near $r_{\rm d}$ (lowest escape speed) but most accretion power is released near radius $r_{\rm isco}$ (highest escape speed), so we expect internal collisions to occur and the total power $L_{\rm d}$ will be thermalized near radius $r_{\rm d}$ under spherical symmetry (Begelman, 2012). Thus, the bulk kinetic energy and the radiation-dominated thermal energy are roughly in equipartition at $r\sim r_{\rm d}$.

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 $v_{\rm w}$ is given by $\dot{M}_{\rm w}v_{\rm w}^{2}/2=L_{\rm d}$, which means

[TABLE]

Note that the asymptotic wind speed is independent of the accretion rate and is sub-relativistic in most situations (unless $s\approx 0$).

The radiation carried by the wind is advected until the photon-trapping radius $r_{\rm tr}$ where photons can diffuse away over a dynamical time. At a given time $t$, a wind shell at radius $r$ was launched at the retarded time $t_{0}(r)$ given by

[TABLE]

where $r_{\rm d}(t_{0})$ and $v_{\rm w}(t_{0})$ are the disk radius (inner boundary of the wind, equation 12) and the wind speed (equation 19) at time $t_{0}$, respectively. Thus, the time evolution of the photon-trapping radius $r_{\rm tr}(t)$ can be estimated by matching the diffusion time with the expansion time,

[TABLE]

where $t_{0}(r)$ is the retarded time of the wind shell currently at radius $r$, $r_{\rm out}=v_{\rm w}(t_{0}=0)t$ is the outer boundary of the wind launched at $t_{0}=0$, and $\tau(r)$ is the Thomson scattering optical depth of the wind outside radius $r$,

[TABLE]

The integrated wind mass loss rate is

[TABLE]

and the wind density profile is taken to be

[TABLE]

where we have taken into account wind propagation by using the retarded time $t_{0}(r)$ for each shell at radius $r$. At the relevant densities and temperatures, the Rosseland-mean opacity is dominated by Thomson scattering $\kappa_{\rm s}=0.34\rm\,cm^{2}\,g^{-1}$ (for cosmic abundance). The time evolution of the photon-trapping radius is shown in Figures 3 and 4 for two TDEs with the same BH mass $M=20M_{\rm\odot}$ but different stellar masses $M_{*}=0.5$ and $5M_{\rm\odot}$, respectively. The qualitative result is that the photon-trapping 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 mass-loss rate drops at late time.

Once we find the photon-trapping radius, the bolometric luminosity of the escaping photons is given by the rate at which radiation are advected across the trapping surface, i.e.

[TABLE]

where the radiation energy density is given by adiabatic expansion (valid at $r\leq r_{\rm tr}$)

[TABLE]

and $U(r_{\rm d})=L_{\rm d}/(8\pi r_{\rm d}^{2}v_{\rm w})$ and $\rho(r_{\rm d})=\dot{M}_{\rm w}/(4\pi r_{\rm d}^{2}v_{\rm w})$ are evaluated using the disk luminosity $L_{\rm d}(t_{0})$, disk radius $r_{\rm d}(t_{0})$, and asymptotic wind speed $v_{\rm w}(t_{0})$ at the retarded time for this trapping shell. The temperature (or average energy) of the escaping photons333The thermalization radius is defined where the effective optical depth $\tau_{*}(r_{\rm th})\simeq\sqrt{\tau_{\rm a}\tau}\simeq 1$ (Rybicki & Lightman, 1979), where $\tau_{\rm a}=\tau_{\rm s}\kappa_{\rm a}/\kappa_{\rm s}$ is the absorption optical depth. In the case when $r_{\rm tr}<r_{\rm th}$, the color temperature is lower than that in equation (27) by a factor of $(r_{\rm th}/r_{\rm tr})^{-3/4}$, but the bolometric luminosity in equation (25) is unaffected. Modeling the (bound-free and free-free) absorption opacity $\kappa_{\rm a}$ requires detailed non-LTE radiative transfer calculations and is left for future works. is given by

[TABLE]

where $\sigma_{\rm SB}$ is the Stefan-Boltzmann constant.

Figures 5 and 6 show the bolometric lightcurve, asymptotic wind velocity, and the temperature of the photons escaping from the trapping surface. We find that BH-MS TDEs generate bright optical transients on timescales of a few to 100 days with peak luminosities from $10^{41}$ to $10^{43}\rm\,erg\,s^{-1}$. The photospheric velocity may range from 0.01$c$ to 0.1$c$ and decreases with time. The largest uncertainty in our model is the mass flux power-law index $s$, and the luminosity decreases towards larger $s$ (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 segment444A steeper power-law generally leads to a more rapidly rising lightcurve, but the latter is much shallower in that $\mathrm{d}\mathrm{log}L/\mathrm{d}\mathrm{log}t<\mathrm{d}\mathrm{log}\dot{M}(r_{\rm d})/\mathrm{d}\mathrm{log}t$, because of the lag between wind launching and photon escaping. The lag time is longer for larger mass flux index $s$ due to a denser wind. of $\dot{M}(r_{\rm d})\propto t^{5}$ in the first $10^{3}\,$sec (or $10^{4}\,$sec) for the case of $M_{*}=0.5M_{\rm\odot}$ (or $M_{*}=5M_{\rm\odot}$), which roughly captures the rapid onset of the gas fallback as shown in Figure 1 of Perets et al. (2016). Input from three-dimensional 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 gamma-ray burst

Given the high peak accretion rate (equation 9), if the accretion efficiency can reach that of typical active galactic nuclei $\sim 0.1$ (e.g. Soltan, 1982), then the peak luminosity of the accretion system is of order $10^{48}\rm\,erg\,s^{-1}$ with duration $\sim 10^{4}$–$10^{5}\,$sec. Thus, it has been proposed (Perets et al., 2016) that TDEs by stellar-mass BHs may be responsible for the population of ultra-long duration gamma-ray bursts (ulGRBs, Gendre et al., 2013; Levan et al., 2014). Numerical simulations of super-Eddington accretion flows in the magnetic-arrested disk (MAD, with dynamically important magnetic flux and rapid BH spin) regime show that high accretion efficiency $>0.1$ can be achieved (Tchekhovskoy et al., 2011; McKinney et al., 2015), and these models have been applied to ultra-luminous X-ray sources (Narayan et al., 2017) and jetted TDEs by super-massive BHs (Tchekhovskoy et al., 2014; Dai et al., 2018; Curd & Narayan, 2019).

Levan et al. (2014) estimated the event rate of ulGRBs to be a factor of a few lower than that of classical GRBs, of order $1\rm\,Gpc^{-3}\,yr^{-1}$ at redshift $z\sim 1$, after correcting for the selection bias that faint, long-lived 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 BH-MS rate predicted by our GC simulations (see Figure 7). Another challenge to the BH-MS TDE scenario is the detection of a bright supernova (SN2011kl) associated with the ultra-long GRB 111209A (Greiner et al., 2015). SN2011kl had a peak time of 14 rest-frame days and peak bolometric luminosity of $\sim 3\times 10^{43}\rm\,erg\,s^{-1}$, and its overall spectral and lightcurve shapes were similar to that of GRB-associated broad-line type-Ic supernovae (with no evidence of H or He). This is inconsistent with the expected H-rich gas composition of BH-MS TDEs, although the super-Eddington disk wind may generate sufficiently bright optical emission (see the $s=0.2$ case in Figure 5).

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.

An alternative possibility mentioned briefly in Perets et al. (2016) is that, on longer timescales, debris from the TDE may slowly fall back, potentially forming a long-lived accretion disk that may produce an X-ray source very similar to an X-ray binary (with the notable difference that the accretion occurs through the disk alone, not fed by a stellar companion). We refer the reader to previous work by Hills (1976) and Krolik (1984) for a more detailed discussion of this scenario.

5 Event Rates

To calculate the event rate of BH–MS TDEs, we use an integral equation similar to that used in Rodriguez & Loeb (2018). The BH–MS TDE rate as a function of cosmic time $t$ is given by

[TABLE]

where $\frac{\dot{M}_{\rm GC}}{d\log_{10}M_{\rm Halo}}$ is the comoving rate of star formation in GCs (in units of $M_{\odot}{\rm yr}^{-1}{\rm Mpc}^{-3}$) per halo mass $M_{\rm Halo}$ at redshift $z(\tau)$ corresponding to a cosmic time $\tau$, from El-Badry et al. (2018). $R(r_{v},M_{\rm GC},t)$ is the rate of BH–MS TDEs at time for for a cluster with an initial mass $M_{\rm GC}$ and virial radius $r_{v}$. As was done in Rodriguez & Loeb (2018), we assume the individual cluster rate can be described as

[TABLE]

where we fit the 5 parameters, $\theta=(A,B,C,\gamma,\gamma_{M})$, separately for clusters with $r_{v}=1\rm{pc}$ and $r_{v}=2\rm{pc}$. 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.

$P(M_{\rm GC})$ is the cluster initial mass function (CIMF), which we assume to be a power-law between $10^{5}M_{\odot}$ and $10^{7}M_{\odot}$, with a possible exponential truncation

[TABLE]

We consider exponential truncations of $M_{\rm GC}^{\bigstar}=\infty$ (corresponding to a purely $\propto M^{-2}$ power law) and $M_{\rm GC}^{\bigstar}=10^{6}M_{\odot}$ (as suggested by observations of young massive star clusters in the local universe, e.g. Portegies Zwart et al. (2010)). $\left<M_{\rm GC}\right>$ 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). $\left<M_{\rm GC}\right>$ is the mean initial mass of GCs given the assumed CIMF.

We have introduced a new parameter, $f_{\rm noIMBH}(r_{v},M_{\rm GC})$, corresponding to the fraction of clusters that do not form an IMBH by runaway collisions of BHs. It has recently been shown Antonini et al. (2018) that for clusters with central escape speeds $\gtrsim 300\rm{km}/\rm{s}$ 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 four-body 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 (28) by assuming that a certain fraction of clusters with a given escape speed $v_{\rm esc}$ contribute zero mergers. We find that an equation of the form

[TABLE]

with $k=0.05$ and $v_{\rm esc}^{\bigstar}=295\rm{km}/\rm{s}$, provides a good fit to Figure 6 of Antonini et al. (2018). We define $v_{\rm esc}$ as the central escape speed from a Plummer model of mass $M_{\rm GC}$ and virial radius $r_{v}$ (e.g., Heggie & Hut, 2003),

[TABLE]

In practice, equation (31) 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 $\sim 10\%$ in the local universe. For comparison, this same correction decreases the BBH merger rate from Rodriguez & Loeb (2018) from $\sim 14\,{\rm Gpc}^{-3}{\rm yr}^{-1}$ to $\sim 12\,{\rm Gpc}^{-3}{\rm yr}^{-1}$ in the local universe. The IMBH correction does not alter the rate when the CIMF is truncated by $M^{\bigstar}_{\rm GC}=10^{6}M_{\odot}$. We discuss the topic of IMBHs further in Section 7.6.

Figure 7 shows the rate as a function of redshift using equation 28. 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 $M^{\star}=10^{6}\,M_{\odot}$ (dashed). In the top panel, we also show in blue (green) the rate if only those models with $r_{v}=1\,$pc ($2\,$pc) are considered in the calculation. As anticipated, the predicted rate is higher for the models with $r_{v}=1\,$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 7 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 7 by this factor. In the local universe ($z<0.1$), we predict a comoving BH–MS TDE rate of $4\,\rm{Gpc}^{-3}\rm{yr}^{-1}$ and an extrapolated upper limit rate of $12\,\rm{Gpc}^{-3}\rm{yr}^{-1}$. If we consider a CIMF truncated at $10^{6}\,M_{\odot}$, these rates decrease to $2.5\,\rm{Gpc}^{-3}\rm{yr}^{-1}$ and $7.5\,\rm{Gpc}^{-3}\rm{yr}^{-1}$, respectively. For comparison, Rodriguez & Loeb (2018) predicted a typical BBH merger rate of $12\,\rm{Gpc}^{-3}\rm{yr}^{-1}$ 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 $(4.4\,\rm{Mpc})^{3}$ (Abadie et al., 2010), we predict a BH–MS TDE rate of $3.4\times 10^{-7}\,\rm{MWEG}^{-1}\rm{yr}^{-1}$ with an upper limit of roughly $10^{-6}\,\rm{MWEG}^{-1}\rm{yr}^{-1}$ in the local universe. Note that the latter rate is in approximate agreement with the analytic predictions made by Perets et al. (2016).

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 1 show these events. In total, we find 246 giant TDEs in models 1-24. 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 here the predicted outcomes from earlier work, and reserve a more detailed study of these events for a future analysis. In particular, Ivanova et al. (2010) showed that close encounters of giants and stellar-mass BHs may serve as a formation channel for BH–WD binaries, an important result given that observations suggest that a fraction of the stellar-mass BH candidates identified in GCs to date may indeed have WD companions (Strader et al., 2012; Bahramian et al., 2017). Estimates from Ivanova et al. (2010) of the hydrodynamics suggest that BH–giant collisions will result in a BH–WD binary with orbital separation $a\gtrsim 1.3\,R_{\rm{G}}$, where $R_{G}$ is the giant’s radius. For giants of mass $\approx 1\,M_{\odot}$ and radii $\approx 3\,R_{\odot}$ this corresponds to $a\gtrsim 10^{-2}\,$AU. Of course, the details of these BH–giant close encounters are sensitive to the highly uncertain common envelope physics and more detailed (hydrodynamic) calculations are necessary to explore the outcomes of these events in detail.

For typical clusters (number densities of $\sim 10^{5}\,\rm{pc}^{-3}$ and velocity dispersions of $\sim 10\rm{km\,s}^{-1}$), Ivanova et al. (2010) predicted BH–WD binaries will form through BH–giant close encounters at a rate of $\sim 1.5\times 10^{-2}$ per BH per Gyr, using a calculation of the form:

[TABLE]

(see, e.g., Ivanova et al., 2005, 2010) where $f_{\rm{G}}$ is the fraction of giants in the stellar population within the GC core, $n$ is the number density, $M_{\rm{BH}}$ is the typical BH mass, $R_{\rm{G}}$ is the typical stellar radius of giants, and $\sigma_{v}$ is the velocity dispersion. $f_{p}=r_{p}/R_{\rm{RG}}$ describes how close a typical BH–giant encounter must be to result in a disruption. In our models, we take $f_{p}=1$, but as discussed in Ivanova et al. (2010), $f_{p}$ as high as $\approx 5$ may also be appropriate.

In Figure 8, we show rate of formation of BH–WDs through BH–giant close encounters computed for various late time cluster snapshots ($t>8\,$Gyr) in models 1-24 plotted against the total number of BHs retained in the cluster at each respective time. We compute the rate, $\Gamma_{\rm{BH-WD}}$ using equation 33 with $M_{\rm{BH}}$, $R_{\rm{G}}$, $n$, $\sigma_{v}$, and $f_{\rm{G}}$ computed uniquely for each cluster snapshot. As in Figure 2, filled circles denote models with $r_{v}=1\,$pc and open circles denote models with $r_{v}=2\,$pc. We limit ourselves to late times here to reflect the typical ages of present-day GCs. BH–WD binaries formed through close encounters of BHs and giants at earlier times may be unlikely to survive to the present as a result of the frequent dynamical encounters which may break apart the binaries,555The typical timescale for a BH–WD binary to undergo a strong encounter that could potentially disrupt the binary can be approximated as $t\sim 3\mathrm{Gyr}\,(n/10^{5}\,\mathrm{pc}^{-3})^{-1}(a/0.1\,\mathrm{AU})^{-1}(M/30\,M_{\odot})^{-1}(\sigma_{v}/10\,\mathrm{km/s})$. Thus we expect a typical BH–WD binary to survive for no longer than a few Gyr in a typical cluster. See, e.g., Kremer et al. (2018b) for further discussion of the rate at which BH–non-BH binaries are expected to be broken through dynamical encounters. thus to remain conservative we include only those binaries that would have formed relatively recently in their host cluster.

Figure 8 shows similar trends to those shown in Figure 2 for the BH–MS TDEs. In particular, the formation rate of BH–WDs varies directly with $N$ (as seen by comparing the different colors) and inversely with cluster size ($r_{v}$; as seen by comparing the open versus filled circles), as expected from equation 33. The dependence of the rate on the total number of retained BHs is more complex. The overall trend is a direct relation between $\Gamma_{\rm{BH-WD}}$ and $N_{\rm{BH}}$, 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 8), it can be seen that $\Gamma_{\rm{BH-WD}}$ varies inversely with $N_{\rm{BH}}$. This results from the way BH populations shape the dynamical evolution of their host cluster. As discussed in Kremer et al. (2019), 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 core-collapse leading to relatively high central densities. Because $\Gamma_{\rm{BH-WD}}$ 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 $r_{v}$: clusters with $r_{v}=1\,$pc (filled circles) retain fewer BHs at late times compared to their $r_{v}=2\,$pc counterparts (open circles). This is also consistent with Kremer et al. (2019), 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 8 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 binary-mediated exchange encounters, as discussed in Kremer et al. (2018b).

Finally, we note that, in addition to being potentially observed as mass-transferring low-mass X-ray binaries, BH–WD binaries may also be observed as gravitational wave (GW) sources by low-frequency GW detectors such as LISA (e.g., Kremer et al., 2018a). 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 stellar-mass 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 stellar-mass BHs disrupt MS stars, giants, and WDs throughout the lifetime of the cluster. These TDEs can occur through both single–single and binary-mediated encounters (binary–single and binary–binary). 2. 2.

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 ($N=2\times 10^{6}$ and $r_{v}=1\,$pc), we get up to 200 BH–MS TDEs over the cluster lifetime. For our lower mass cluster models ($N=2\times 10^{5}$), the total number of TDEs can be as low as zero to a few. 3. 3.

By incorporating a realistic cosmological model for GC formation, we derive a rate of BH–MS TDEs of approximately $4-12\,\rm{Gpc}^{-3}\,\rm{yr}^{-1}$ in the local universe and a cosmological rate that peaks at roughly $25-75\,\rm{Gpc}^{-3}\,\rm{yr}^{-1}$ at a redshift of about 3. 4. 4.

We show that the wind mass loss associated with these BH–MS TDEs produces optical transients of luminosity $\sim 10^{41}$ to $10^{44}\rm\,erg\,s^{-1}$ on timescales of about a day to a month. In Section 4 we show lightcurve predictions expected from these events. 5. 5.

BH–giant close encounters occur at rates of up to $\sim 1$ per Gyr per cluster. These events may serve as a dynamical formation channel for BH–WD binaries in GCs, which may be observed as X-ray or gravitational wave sources.

7.2 Detectability

The Zwicky Transient Facility (ZTF) reaches $g$ = 20.8 mag (5$\sigma$) during a single exposure of $30\,$s, surveying $3\pi$ of the sky. Assuming a blackbody temperature of $2\times 10^{4}\,$K (see Figures 5 and 6), the detection horizon in luminosity distance is roughly $D_{\rm L,max}\simeq 150L_{\rm bol,42}\rm\,Mpc$, where $L_{\rm bol}=10^{42}L_{\rm bol,42}^{1/2}\rm\,erg\,s^{-1}$ is the bolometric luminosity near the optical peak. The all-sky rate is about $10\rm\,yr^{-1}$ in the optimistic case where $s=0.2$ and peak luminosity $L_{\rm bol}\sim 10^{44}\rm\,erg\,s^{-1}$. In the pessimistic case where $s=0.8$ and peak luminosity $L_{\rm bol}\sim 10^{42}\rm\,erg\,s^{-1}$, the all-sky rate is about $10^{-2}\rm\,yr^{-1}$.

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 $5000\,\rm{km\,s^{-1}}$ (e.g., Campanelli et al., 2007; Lousto et al., 2012), 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 second-generation mergers in clusters (see, e.g., Rodriguez et al., 2018b).

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 $\mathcal{E}\in(0,\mathcal{E}_{\rm{max}})$, where $\mathcal{E}_{\rm{max}}=v_{\rm{max}}^{2}/2$ is given by

[TABLE]

where we have used $r_{\rm TDE}=R_{\star}(M_{\rm{BH}}/M_{\star})^{1/3}$. The escape velocity of the star is given by $v_{\rm{esc}}=\sqrt{2GM_{\star}/R_{\star}}$. The maximum speed of the unbound debris is then $v_{\rm{max}}=(M_{\rm{BH}}/M_{\star})^{1/6}v_{\rm{esc}}$. The total linear momentum, $P$, carried away by the unbound debris depends on the mass distribution over specific energy, $\rm{d}M/\rm{d}\mathcal{E}$, but a rough estimate is that $P\simeq(M_{\star}/3)v_{\rm{max}}$ (this expression is exact for a flat distribution of $\rm{d}M/\rm{d}\mathcal{E}$). Then the BH receives a kick in the $-P$ direction of velocity

[TABLE]

For a typical TDE with a MS star of mass of $0.5\,M_{\odot}$ (with escape velocity of $\sim 600\,$km/s) and a BH of mass of $20\,M_{\odot}$, we obtain $v_{\rm{kick}}\simeq 10\,$km/s. These kicks are low compared to both typical cluster escape velocities ($\sim 50-100\,\rm{km\,s}^{-1}$) as well as typical dynamical recoil kicks attained from small-$N$ 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 spin-up as a result of these TDEs is beyond the present study. Follow-up 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., Hansen & Murali (1998) and Perets et al. (2016). In particular, Hansen & Murali (1998) 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 1-24. Roughly, this translates to a galactic rate of $\sim 10^{-8}\,\rm{yr}^{-1}$ per Milky Way-like 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 (e.g., Ye et al., 2019; Belczynski et al., 2018), 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 core-collapsed cluster with relatively higher stellar densities and $\sim 0$ BHs (see, e.g., Kremer et al., 2019), NS–MS TDEs may become more frequent.

7.5 Link to $r$-process enrichment?

Observational evidence of $r$-process enrichment exists in several Milky Way GCs, including M5, M15, M92, and NGC 3201 (e.g., Roederer, 2011; Bekki, 2018). The LIGO/Virgo detection of the binary NS merger GW170817 (Abbott et al., 2017) and the follow-up electromagnetic observations showed that NS mergers produce large amounts of $r$-process elements (Kasen et al., 2017). However, explaining the observed $r$-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., Gratton et al., 2012). For a particular binary NS merger event (or series of events) to enrich a GC’s stars with $r$-process material, the event(s) must occur while star formation is still taking place (see, e.g., Bekki (2018); Zevin et al. (2019) for further discussion).

Merging binary NSs are expected to form in GCs through two mechanisms: binary evolution of primordial binaries (e.g., Ivanova et al., 2003; Dominik et al., 2012; Tauris et al., 2017) and dynamical formation of NSs at late times (e.g., Ye et al., 2019). In the first scenario, where the NS components typically form through iron core-collapse supernovae (which are expected to lead to large natal kicks; Hobbs et al., 2005), 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 (Zevin et al. (2019); also see Safarzadeh et al. (2018) 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 mass-segregate to the cluster core, which takes $\sim\,$Gyrs, due to the relatively low NS masses (see, e.g., Ye et al., 2019; Zevin et al., 2019). 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 $r$-process material produced in these late-time 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 neutron-rich material, these events may provide a way to enrich GCs with $r$-process elements at early times. As shown in Kremer et al. (2018c), the cluster NGC 3201, one of the MW clusters in which $r$-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 Gratton et al. (2012) 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) $r$-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 $r$-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 (Siegel et al., 2018) may serve as viable alternatives for $r$-process production in GCs. A complete exploration of $r$-process production in GCs is beyond the scope of this paper; see Zevin et al. (2019) 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., Samsing et al., 2017, 2018) has shown that the inclusion of tidal coupling within small-$N$ dynamical encounters may have a significant effect on tidal disruptions. Many of these complexities will require detailed 3-D (magneto)hydrodynamic calculations of these encounters. Our first exploration of some of these complexities will be presented in a forthcoming paper (Fixelle et al., 2019).

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 binary-mediated resonant encounters, is necessary to explore outcomes of these events in greater detail. In a future study, we hope to explore binary-mediated close encounters between BHs and stars using a small $N$ integrator that self-consistently incorporates relevant hydrodynamic effects to examine some of these complexities in more detail (for previous work on the topic, see, for example, Goodman & Hernquist, 1991).

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 1). This is perhaps surprising, given earlier work showing that the number of close encounters can increase substantially for binary-mediated encounters relative to single–single encounters alone (e.g., Bacon et al., 1996; Fregeau et al., 2004; Chatterjee et al., 2013). For example, Chatterjee et al. (2013) 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 Kremer et al. (2018b), 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 Chatterjee et al. (2013), the models with the most marked increase in binary-mediated close encounters in that analysis are those with initial global binary fractions as high as $27\%$ (and as high as roughly $40\%$ in the core). Here we consider models with lower binary fraction ($10\%$). 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.

In this analysis, we have focused primarily on stellar-mass BHs, however a large body of literature has also explored the possible existence of IMBHs in GCs from both observational (e.g., Gebhardt et al., 2005; Lanzoni et al., 2007; Lützgendorf et al., 2011; Feldmeier et al., 2013) and theoretical perspectives (e.g., Portegies Zwart et al., 2004; Freitag et al., 2006; Giersz et al., 2015; Mapelli, 2016; Antonini et al., 2018). TDEs (particularly TDEs of WDs) by IMBHs have been examined in various recent analyses including Rosswog et al. (2008, 2009); MacLeod et al. (2014, 2016) and these events have been proposed as the mechanisms that may have produced a number of observed high-energy events (e.g., Krolik & Piran, 2011; Jonker et al., 2013). The existence of IMBHs in clusters remains contested. Nonetheless, the role that IMBHs, if present, may play in the production of high-energy transients like TDEs and also GC dynamics more broadly is a rich topic that we hope to explore in more detail within the scope of our Monte Carlo code in a later study.

We thank Sasha Tchekhovskoy, Pablo Marchant, Johan Samsing, Michael Zevin, and Enrico Ramirez-Ruiz 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 AST-1716762. KK acknowledges support by the National Science Foundation Graduate Research Fellowship Program under Grant No. DGE-1324585. WL is supported by the David and Ellen Lee Fellowship at Caltech.