Fast Optical Transients from Stellar-mass Black Hole Tidal Disruption Events in Young Star Clusters

We first discuss the case of TDEs occurring during single–single encounters between a BH and an MS star. For a given cluster with N total stars and half-mass radius r_h, the rate of TDEs occurring through single–single dynamical encounters can be estimated as

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{\mathrm{ss}}\approx {n}_{\star }{{\rm{\Sigma }}}_{\mathrm{ss}}{\sigma }_{v}{N}_{\mathrm{BH}},\end{eqnarray} \tag{ 3 }$$

where N_BH is the total number of BHs in the cluster (which, in general, are all found within the half-mass radius due to mass segregation; e.g., Morscher et al. 2015), ${n}_{\star }\approx N/{r}_{{\rm{h}}}^{3}$ is the number density, and σ_v is the cluster's velocity dispersion. Σ_ss is the cross section for a single–single TDE, given by

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{\mathrm{ss}}=\pi {R}_{\mathrm{TD}}^{2}\left[1+\displaystyle \frac{2G({m}_{\mathrm{BH}}+{m}_{\star })}{{R}_{\mathrm{TD}}{\sigma }_{v}^{2}}\right],\end{eqnarray} \tag{ 4 }$$

where m_BH and m_⋆ are the typical BH and stellar masses and R_TD is the tidal disruption radius given by Equation (1). For a cluster with a Kroupa (2001) IMF, we can take m_⋆ ≈ 0.6 M_⊙. For high-metallicity clusters (Z ≈ Z_⊙), we can take m_BH ≈ 10 M_⊙, while for low-metallicity clusters (Z ≲ 0.1 Z_⊙), m_BH ≈ 25 M_⊙ is more appropriate (e.g., Kremer et al. 2020b).

Assuming the cluster is initially in virial equilibrium (${\sigma }_{v}\approx \sqrt{{{GM}}_{\mathrm{cl}}/{r}_{{\rm{h}}}}$, where M_cl ≈ m_⋆ N is the total cluster mass), assuming all encounters occur in the gravitational focusing regime (the second term in the brackets of Equation (4) dominates), and taking m_BH ≫ m_⋆, we can rewrite Equation (3) as

$$\begin{eqnarray}\begin{array}{rcl}{{\rm{\Gamma }}}_{\mathrm{ss}} & \approx & 0.2\,{\mathrm{Gyr}}^{-1}{\left(\displaystyle \frac{{M}_{\mathrm{cl}}}{{10}^{4}{M}_{\odot }}\right)}^{1/2}{\left(\displaystyle \frac{{r}_{{\rm{h}}}}{1\mathrm{pc}}\right)}^{-5/2}\\ & & \times \ {\left(\displaystyle \frac{{m}_{\star }}{0.6{M}_{\odot }}\right)}^{-4/3}{\left(\displaystyle \frac{{m}_{\mathrm{BH}}}{10{M}_{\odot }}\right)}^{4/3}\left(\displaystyle \frac{{R}_{\star }}{0.6\,{R}_{\odot }}\right)\left(\displaystyle \frac{{N}_{\mathrm{BH}}}{10}\right).\end{array}\end{eqnarray} \tag{ 5 }$$

For a Kroupa (2001) IMF, we expect roughly 10⁻³ BHs to form per star, giving us N_BH ≈ 10⁻³ N. In this case, we obtain the analytic scaling

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{\mathrm{ss}}\approx 0.2\,{\mathrm{Gyr}}^{-1}{\left(\displaystyle \frac{{M}_{\mathrm{cl}}}{{10}^{4}{M}_{\odot }}\right)}^{3/2}.\end{eqnarray} \tag{ 6 }$$

In Figure 1, we show as open circles the rates of such encounters as a function of cluster mass as determined from our N-body models. To isolate specifically the effect of cluster mass, we utilize only the first eight sets of simulations listed in Table 1 (models a-h), which have fixed metallicity (Z_⊙) and virial radius (1 pc). To calculate the model rates, we simply count the total number of single–single TDEs in all models of a given cluster mass, then divide by the total number of models and by the total integration time for the given cluster mass (see Table 1). Error bars denote 2σ from the mean, assuming a Poisson distribution.

Zoom In Zoom Out Reset image size

From a least squares fit, we find these data are best fit by the power-law relation

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{\mathrm{ss}}^{\mathrm{model}}=0.07\pm 0.01\,{\mathrm{Gyr}}^{-1}{\left(\displaystyle \frac{{M}_{\mathrm{cl}}}{{10}^{4}{M}_{\odot }}\right)}^{1.6\pm 0.16}.\end{eqnarray} \tag{ 7 }$$

This fit is shown as the blue curve in Figure 1, with the blue bands denoting the 90% confidence interval from the least squares fit. For comparison, we show as the gray dashed line in Figure 1 the ${\rm{\Gamma }}\propto {M}_{\mathrm{cl}}^{3/2}$ scaling relation derived from the simple analytic estimate in Equation (6).

In Equation (6), we have assumed a constant r_h ≈ 1 pc is typical for all YSCs (see, e.g., Portegies Zwart et al. 2010). Alternatively, as discussed in Section 2, r_h < 1 pc may be more appropriate and furthermore, r_h may exhibit a weak dependence on total cluster mass. To explore the possibility, we ran an additional set of models (group l in Table 1) with initial r_v = 0.42 pc, reflective of the r_h–M_cl relation of Marks & Kroupa (2012). Given the phenomenological relation from Marks & Kroupa (2012) predicts cluster radii a factor of roughly 2 lower than the r_h = 1 pc assumption, the following rate, ${{\rm{\Gamma }}}_{\mathrm{ss}}^{\mathrm{phenom}}$, is higher than that estimated from Equation (7).

Combining the phenomenological r_h–M_cl relation of Marks & Kroupa (2012) (Equation (2)) with Equations (3) and (4), we expect ${{\rm{\Gamma }}}_{\mathrm{ss}}^{\mathrm{phenom}}\propto {M}_{\mathrm{cl}}^{1.175}$. Scaling to the rate identified from the models in group l, we can then write

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{\mathrm{ss}}^{\mathrm{phenom}}\approx 0.63\,{\mathrm{Gyr}}^{-1}{\left(\displaystyle \frac{{M}_{\mathrm{cl}}}{{10}^{4}{M}_{\odot }}\right)}^{1.175}.\end{eqnarray} \tag{ 8 }$$

In addition to TDEs occurring through single–single encounters, TDEs may also take place during binary-mediated resonant encounters involving at least one BH and one star. As discussed in Section 2, we do not include stellar binaries in this study. The only binaries formed are BH binaries assembled through three-body encounters (for simplicity, three-body binary formation is allowed only for BHs in our models; e.g., Morscher et al. 2015). In this case, the binary-mediated TDEs discussed in this subsection are those that occur specifically during binary–single resonant encounters between a BBH and a single MS star.

Using a similar calculation to that performed for the single–single rate estimate, the rate of TDEs during BBH–star binary–single encounters can be written as

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{\mathrm{bs}}\approx {n}_{\star }\,{{\rm{\Sigma }}}_{\mathrm{bs}}\,{\sigma }_{v}\,{N}_{\mathrm{BBH}}\,{P}_{\mathrm{TD}}.\end{eqnarray} \tag{ 9 }$$

Here, Σ_bs is the cross section for binary–single encounters

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{\mathrm{bs}}=\pi {a}_{\mathrm{BBH}}^{2}\left[1+\displaystyle \frac{2G({m}_{\mathrm{BBH}}+{m}_{\star })}{{a}_{\mathrm{BBH}}{\sigma }_{v}^{2}}\right],\end{eqnarray} \tag{ 10 }$$

where a_BBH is the BBH semimajor axis and m_BBH ≈ 2m_BH is the mass of the BBH. N_BBH is the total number of BBHs in the cluster. As shown in a number of recent analyses (e.g., Morscher et al. 2015; Chatterjee et al. 2017; Banerjee 2018), N_BBH is expected to be roughly independent of the total number of BHs in the cluster as well as the cluster's total mass, such that the total number of dynamically formed BBHs present at any given time never exceeds a few. Here, we assume N_BBH ≈ 2.

Finally, P_TD is the probability that a given BBH–star resonant encounter leads to a TDE. We follow the results of Darbha et al. (2018), which showed that for asymmetric mass ratio binary–single encounters, P_TD is roughly proportional to R_TD/a_BBH. Note that this same scaling is found in the equal mass case (e.g., Samsing et al. 2017; Samsing 2018). Here, we take P_TD = 2R_TD/a_BBH as in Samsing et al. (2019). We can then rewrite Equation (9) as

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{\mathrm{bs}}\approx 0.3\,{\mathrm{Gyr}}^{-1}{\left(\displaystyle \frac{{M}_{\mathrm{cl}}}{{10}^{4}{M}_{\odot }}\right)}^{1/2},\end{eqnarray} \tag{ 11 }$$

where, as before, we have assumed r_h ≈ 1 pc, m_BH ≈ 10 M_⊙, m_⋆ ≈ 0.6 M_⊙, and R_⋆ ≈ 0.6 R_⊙, independent of the cluster mass.

Comparing to Equation (6), we find ${{\rm{\Gamma }}}_{\mathrm{bs}}/{{\rm{\Gamma }}}_{\mathrm{ss}}\propto {N}_{\mathrm{BH}}^{-1}\propto {M}_{\mathrm{cl}}^{-1}$. Thus, lower-mass clusters with fewer BHs feature a higher BBH TDE rate per BH (and therefore also per star) compared to higher-mass clusters with more BHs. Thus, given that lower-mass clusters also dominate by number the overall cluster mass function, they are an ideal environment to find BBH TDEs. We return to the question of overall and relative rates in Section 3.3.

In Figure 2, we show the rates of these events as a function of cluster mass as determined from our N-body models. As in Figure 1, the rate is computed as the total number of binary–single TDEs in all models of a given cluster mass, divided by the total number of models and by the time duration, Δt. Unlike in the single–single case, where the necessary dynamical encounters begin immediately, in the binary–single case we must first wait for BHs to mass segregate to the center so that the target BBHs can form through three-body encounters. For each model, we define this timescale, t_3bb, simply as the moment the first BBH forms in that model. Then, ${\rm{\Delta }}t={t}_{\max }-{t}_{3\mathrm{bb}}$. Typically, t_3bb is of order 100 Myr (see also Sigurdsson 1993; Morscher et al. 2015).

Zoom In Zoom Out Reset image size

Again performing a least squares fit, we find these data are best fit by the power-law relation

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{\mathrm{bs}}^{\mathrm{model}}=0.69\pm 0.36\,{\mathrm{Gyr}}^{-1}{\left(\displaystyle \frac{{M}_{\mathrm{cl}}}{{10}^{4}{M}_{\odot }}\right)}^{0.3\pm 0.09},\end{eqnarray} \tag{ 12 }$$

which can be compared with the simple analytic estimate of Equation (11). As with the single–single case, we find the analytic estimate recovers reasonably well the rate inferred from the N-body modeling.

As before, we can also estimate a phenomenological rate, ${{\rm{\Gamma }}}_{\mathrm{bs}}^{\mathrm{phenom}}$, assuming the relation of Marks & Kroupa (2012). Combining Equations (2) and (11), we expect ${{\rm{\Gamma }}}_{\mathrm{bs}}^{\mathrm{phenom}}\propto {M}_{\mathrm{cl}}^{0.175}$. Again normalizing to the binary–single TDE rate estimated in the models from group l in Table 1, we can write

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{\mathrm{bs}}^{\mathrm{phenom}}\approx 4.4\,{\mathrm{Gyr}}^{-1}{\left(\displaystyle \frac{{M}_{\mathrm{cl}}}{{10}^{4}{M}_{\odot }}\right)}^{0.125}.\end{eqnarray} \tag{ 13 }$$

Because we assume that BH binaries form exclusively through three-body encounters, Equations (12) and (13) are relevant only for those clusters with at least three BHs at birth, so that at least one BBH can form. Again, assuming roughly 1 BH forms per 1000 stars (Kroupa 2001), this requires N ≳ 3000 or M_cl ≳ 2000 M_⊙. For clusters with M_cl ≲ 2000 M_⊙, Equations (12) and (13) are no longer applicable and the binary–single TDE rate is zero. We return to this point in Section 3.3 when we compute the total TDE rate by integrating over the full cluster mass function.

3.2.1. BBH Orbital Separations

Given that the binary–single TDE rate is expected to be roughly independent of the binary orbital separation, the distribution of a_BBH for BBHs that undergo TDEs is expected to follow the semimajor axis for all BBHs in a cluster of a given mass (e.g., Kremer et al. 2019b; Samsing et al. 2019). In Figure 3, we show the distribution of a_BBH for all BBHs found in models a-h (various initial masses, but fixed r_v and Z). These distributions are determined by two primary physical processes:

Zoom In Zoom Out Reset image size

Three-body formation: The maximum semimajor axis of a binary formed through a three-body encounter is determined by the hard-soft boundary ${a}_{\mathrm{HS}}\propto {m}_{\mathrm{BH}}/{\sigma }_{v}^{2}$ (e.g., Morscher et al. 2015). Assuming ${\sigma }_{v}^{2}\propto {M}_{\mathrm{cl}}$, we expect, in general, more massive clusters will produce more compact BBHs compared to lower-mass clusters (assuming a fixed cluster virial radius).

Ejection from dynamical recoil: Once a BBH is formed, it will (on average) harden through subsequent binary-mediated encounters with other BHs and stars in the cluster core (e.g., Sigurdsson 1993; Morscher et al. 2015; Rodriguez et al. 2016). Following a dynamical encounter, the BBH will receive a dynamical recoil kick with a magnitude comparable to the BBH orbital velocity, ${v}_{\mathrm{recoil}}^{2}\propto {a}_{\mathrm{BBH}}^{-1}$. Thus, as a BBH hardens, it attains increasingly large dynamical recoil kicks. Eventually, v_recoil is sufficiently large for the binary to be ejected from the cluster. This is set by the cluster's escape velocity, ${v}_{\mathrm{esc}}^{2}\propto {M}_{\mathrm{cl}}/{r}_{{\rm{h}}}$. As a result, in lower-mass clusters, BBHs will be dynamically ejected before they can harden as far as is possible in higher-mass clusters.

As a consequence of these two processes, we expect the BBH semimajor axis distribution to shift toward lower values in increasingly massive clusters. Indeed, this is shown in Figure 3.

One exciting possibility proposed in Lopez et al. (2020), Samsing et al. (2019), and Kremer et al. (2019b) is that these binary-mediated TDEs may be used to indirectly probe properties of the underlying BBH population, if the corresponding EM signal can be detected. The basic idea is the second BH produces breaks in the light curve on timescale comparable to the BBH orbital period. This idea has been illustrated in the supermassive BH (SMBH) regime using numerical techniques (e.g., Liu et al. 2009; Coughlin et al. 2017), and one SMBH candidate has proposed (TDE J1201+30), from which the authors were able to constrain the SMBH binary orbital period (Liu et al. 2014). This process likely requires a_BBH to be comparable to the disk radius, which in turn will be comparable to the tidal disruption radius. For reference, we show r_TD ≈ 10 R_⊙ as a solid black line in Figure 3 (see Equation (1)). As shown in the figure, this process is likely only possible in the most massive clusters explored here (M_cl ≳ 10⁵ M_⊙). Even for our massive cluster simulations, we find that only ≈0.1% of all BBHs meet this criterion. Thus, we conclude that the presence of a second BH is unlikely to significantly affect the TDE dynamics and subsequent light-curve evolution.

Although it appears this possibility is not relevant in typical YSCs, we can speculate that more massive clusters such as nuclear star clusters with masses of 10⁷ M_⊙ or larger may be ideal environments for this processes, given that more massive clusters should be able to host even more compact BBHs. We reserve a more detailed study of this possibility for future study, and direct the reader to Fragione et al. (2020) for a discussion of TDEs in the nuclear star cluster regime.

The functional form for the initial mass function of YSCs is expected to be well represented by a power-law distribution (e.g., Lada & Lada 2003) with a possible exponential truncation above cluster masses of roughly M_cut ≈ 10⁶ M_⊙ (e.g., Portegies Zwart et al. 2010):

$$\begin{eqnarray}&&\displaystyle \frac{{dN}}{{{dM}}_{\mathrm{cl}}}\propto {M}_{\mathrm{cl}}^{-2}\exp (-{M}_{\mathrm{cl}}/{M}_{\mathrm{cut}}).\end{eqnarray} \tag{ 14 }$$

As we are interested here primarily in the low-mass tail of the mass function (≲10⁵ M_⊙), the specific value of M_cut is not relevant to this study.

We can then compute the total TDE rate from a realistic population of YSCs in Milky Way–like galaxies by integrating the rate per cluster (Equation (5)) over the cluster mass function:

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{\mathrm{tot}}={\int }_{{M}_{\mathrm{low}}}^{{M}_{\mathrm{high}}}\displaystyle \frac{{{\rm{\Gamma }}}_{\mathrm{cl}}}{{M}_{\mathrm{cl}}}\,\displaystyle \frac{{dN}}{{{dM}}_{\mathrm{cl}}}\,{\rm{\Delta }}t\,{\rho }_{\mathrm{SF}}\,{f}_{\mathrm{SF}}\,{{dM}}_{\mathrm{cl}}.\end{eqnarray} \tag{ 15 }$$

The integration limits represent the assumed range in cluster masses; we assume M_low = 100 M_⊙ and M_high = 10⁵ M_⊙ (Lada & Lada 2003). By definition, in order for TDEs to occur through the binary–single channel discussed in Section 3.2, BBHs must have formed in the cluster. As discussed in Section 3.2, because we assume that BBHs form exclusively through three-body encounters, this requires at least three BHs be present in the cluster at birth, which in turn requires M_cl ≳ 2000 M_⊙. Thus, to compute the rate of TDEs occurring through binary–single encounters, we use M_low = 2000 M_⊙ (keeping the overall normalization of Equation (14) the same as before). Note that this mass requirement automatically takes care of the additional requirement that the cluster not dissolve before the first BBHs begin to form.

Δt is the cluster disruption timescale, t_dis, which depends upon the location of the cluster in its host galaxy's tidal field, as well as on more complex phenomena such as, e.g., tidal shocks and interactions with giant molecular clouds as mentioned in Section 2. Here, we adopt the following relation from Lamers et al. (2005):

$$\begin{eqnarray}&&{t}_{\mathrm{dis}}\simeq 810{\left(\displaystyle \frac{{M}_{\mathrm{cl}}}{{10}^{4}{M}_{\odot }}\right)}^{0.62}{\left(\displaystyle \frac{{\rho }_{\mathrm{amb}}}{{M}_{\odot }{\mathrm{pc}}^{-3}}\right)}^{-0.5}\,\mathrm{Myr},\end{eqnarray} \tag{ 16 }$$

where ρ_amb ≈ M_⊙ pc⁻³ is the assumed local ambient density.

Again, as we are specifically interested in young clusters with ages less than roughly 500 Myr, we take

$$\begin{eqnarray}&&{\rm{\Delta }}t=\min \left({t}_{\mathrm{dis}},500\,\mathrm{Myr}\right).\end{eqnarray} \tag{ 17 }$$

For TDEs occurring during binary–single encounters, we must also incorporate the timescale for BBH formation to begin. In Section 3.2, we computed this timescale directly from the models. Here, for simplicity (and as motivated by the results from the models), we assume BBH formation occurs after roughly 100 Myr for all cluster masses (see also, e.g., Sigurdsson 1993). In this case, for binary–single TDEs, ${\rm{\Delta }}t=\min ({t}_{\mathrm{dis}},500\,\mathrm{Myr})-100\,\mathrm{Myr}$.

ρ_SF is the assumed cosmological density of the star formation rate (SFR). We adopt the SFR of Hopkins & Beacom (2006). Specifically, this study finds ρ_SF/[M_⊙ yr⁻¹ Mpc⁻³] = 1.5 × 10⁻², 0.1, and 0.2 at redshift z = 0, 1, and 2.5 (peak star formation), respectively. f_SF is the fraction of the SFR assumed to occur in star clusters. For low-mass clusters (M_cl < 10⁵ ⊙ ) we assume f_SF = 0.8 (Lada & Lada 2003).

Finally, Γ is the TDE rate per cluster of a given mass, M_cl. For this we adopt the scaling relations derived in Section 3 (Equations (7) and (12)), for the single–single and binary–single cases, respectively, assuming constant r_h = 1 pc). Additionally, we use Equations (8) and (13) to compute the rates assuming higher-density YSCs as in the phenomenological fits of Marks & Kroupa (2012).

We present in Table 2 the rate estimates obtained by integrating Equation (15). We find that the binary–single channel dominates over the single–single channel by a factor of roughly a few to 10, depending upon the assumptions made regarding r_h. If we adopt the more conservative choice of constant r_h = 1 pc, we estimate a combined TDE rate of roughly 20 Gpc⁻³ yr⁻¹ in the local universe. Adopting the phenomenological assumption from Marks & Kroupa (2012), we find a combined TDE rate of roughly 200 Gpc⁻³ yr⁻¹.

For reference, Kremer et al. (2019c) predicted a TDE rate of roughly 10 Gpc⁻³ yr⁻¹ for old globular clusters. In the more massive nuclear star cluster regime, Fragione et al. (2020) predicted a stellar-mass BH TDE rate of roughly 10⁻⁷–10⁻⁶ yr⁻¹ per galaxy. Assuming a galactic density of roughly 10⁻² Mpc⁻³, this corresponds to a rate of roughly 1–10 Gpc⁻³ yr⁻¹ in the local universe. Thus, we conclude YSCs may dominate the overall stellar-mass BH TDE rate in the local universe by a factor of a few to more than an order of magnitude compared to more massive clusters.

In the previous subsections, we have explored specifically models a-h of Table 1, which assume solar metallicity, reflective of YSCs born recently in the local universe. However, for YSCs found at higher redshifts, assuming a lower metallicity is more appropriate. Given the overall TDE rate to be substantially higher at high redshift (see Section 3.3), a careful investigation of metallicity is warranted.

To explore the effect of cluster metallicity on TDE properties, we have run several additional sets of models with Z = 0.1 Z_⊙. In Figure 4, we show the distribution of BH masses that undergo TDEs. The black and hatched gray histograms show BH masses found in models assuming r_h = 1 pc and the Marks & Kroupa (2012) r_h–M_cl relation, respectively. In the top (bottom) panel we show the results for models assuming solar (10% solar) metallicity.

Zoom In Zoom Out Reset image size

The first general feature we see is that higher-metallicity clusters yield lower-mass BH TDEs. This is anticipated: at higher metallicity, stellar winds are expected to lead to increased mass loss prior to stellar core collapse, which in turn is expected to reduce the mass of the BH ultimately formed (e.g., Vink et al. 2001; Fryer et al. 2012; Belczynski et al. 2016). For the r_v = 1 pc models, we find median BH TDE masses of 13 M_⊙ and 27 M_⊙ for Z_⊙ and 0.1 Z_⊙, respectively.

We also see that, for a given metallicity, increasing the initial cluster density yields an extended tail in the upper part of the BH mass spectrum. This high-mass tail is populated by BHs formed through stellar collisions, which occur at an increased rate in higher-density clusters. A number of recent analyses (Spera & Mapelli 2017; Di Carlo et al. 2019a, 2019b; Kremer et al. 2020a) have demonstrated that dynamically mediated stellar collisions occurring within the first roughly 5 Myr of cluster evolution (before formation of BHs) may lead to formation of massive stars that may ultimately collapse to form high-mass BHs. In particular, this process may permit formation of BHs with masses occupying the pair-instability mass gap from roughly 40–120 M_⊙ expected to arise through (pulsational) pair-instability supernovae (e.g., Belczynski et al. 2016; Spera & Mapelli 2017; Woosley 2017). Additionally, this process may be closely related to collisional runways, which have been touted as a potential formation mechanism for intermediate-mass BHs (IMBHs) with masses in excess of roughly 100 M_⊙ (e.g., Portegies Zwart et al. 2004; Gürkan et al. 2006; Giersz et al. 2015; Mapelli 2016). Here, we adopt the same prescriptions implemented in Kremer et al. (2020a) to treat stellar collisions and the subsequent evolution of the collision products, as described in Section 2. We also assume the maximum BH mass formed through single-star evolution (i.e., unaffected by any dynamical processes) is 40.5 M_⊙, as determined by our (pulsational) pair-instability supernova treatment (see Belczynski et al. 2016, for details).

In our Z_⊙ and 0.1 Z_⊙ models, which adopt the Marks & Kroupa (2012) r_h relation, we find roughly 10% of all TDEs have BH masses within the assumed pair-instability gap. In the 0.1 Z_⊙ models specifically, we find an additional 5% of TDEs occur with a BH mass in excess of 120 M_⊙ (the assumed upper limit to the pair-instability gap). These fractions are consistent with the rates of formation of massive BHs shown in previous studies of YSCs (Di Carlo et al. 2019a, 2019b). Given the cosmological rates predicted in Section 3.3, these mass-gap and IMBH TDEs may constitute a non-negligible fraction of all TDEs occurring in YSCs. Furthermore, these TDEs may provide a potentially novel way to probe the formation of pair-instability gap BHs, similar to the recent LIGO/Virgo detection GW190521 (Abbott et al. 2020a, 2020b).

Finally, comparing models i,j,k with models e,f,g, we see that lower-metallicity clusters exhibit a moderate increase (a factor of roughly 1.4) in the total number of TDEs occurring through both single–single and binary–single encounters. This slight increase in the rate is anticipated: as shown in Equation (5), the TDE rate scales with the BH mass, m_BH, through the influence of gravitational focusing and through tidal disruption radius calculation (Equation (1)). Thus, we expect that metallicity leads to a moderate difference in the TDE rate for a given cluster mass.

Following the disruption of a star of mass m_⋆ and radius R_⋆, the timescale, t_fb for (bound) orbiting material to fallback to the disruption point, R_TD is given by

$$\begin{eqnarray}\begin{array}{rcl}{t}_{\mathrm{fb}} & = & \displaystyle \frac{\pi {R}_{\mathrm{TD}}^{3}}{\sqrt{2{{Gm}}_{\mathrm{BH}}{R}_{\star }^{3}}}\\ & & \approx \ 1.1\times {10}^{4}\,{\rm{s}}{\left(\displaystyle \frac{{m}_{\star }}{{M}_{\odot }}\right)}^{-1}{\left(\displaystyle \frac{{R}_{\star }}{{R}_{\odot }}\right)}^{3/2}{\left(\displaystyle \frac{{m}_{\mathrm{BH}}}{10{M}_{\odot }}\right)}^{1/2}\end{array}\end{eqnarray} \tag{ 18 }$$

(e.g., Perets et al. 2016).

For simplicity, we assume a disk is formed promptly at radius r_d ≃ R_TD and take t_fb as the characteristic timescale for disk formation. This is motivated by the fact that the orbits of the bound debris are only weakly eccentric. Once a disk is formed, the timescale for the debris to accrete is set by the viscous timescale (Shakura & Sunyaev 1973). For a thick disk, with disk height ratio h = H/r_d (where H is the disk scale height), the viscous accretion timescale is

$$\begin{eqnarray}\begin{array}{rcl}{t}_{\mathrm{acc}} & \approx & {\left[{h}^{2}\alpha {{\rm{\Omega }}}_{{\rm{K}}}({r}_{d})\right]}^{-1}\\ & & \approx \ 6\times {10}^{4}\,{\rm{s}}{\left(\displaystyle \frac{h}{0.5}\right)}^{-2}{\left(\displaystyle \frac{\alpha }{0.1}\right)}^{-1}{\left(\displaystyle \frac{{m}_{\star }}{{M}_{\odot }}\right)}^{-1/2}{\left(\displaystyle \frac{{R}_{\star }}{{R}_{\odot }}\right)}^{3/2}.\end{array}\end{eqnarray} \tag{ 19 }$$

Because t_acc > t_fb, the subsequent light-curve evolution is viscosity driven (i.e., determined by accretion timescale). This marks one key difference from SMBH TDEs, where t_acc ≪ t_fb and thus, the disk evolution is likely dominated by the fallback and the accretion rate is believed to follow the standard t^−5/3 power law (e.g., Rees 1988; Evans & Kochanek 1989; Phinney 1989).

For viscosity-driven accretion, assuming roughly half of the stellar material is bound to the BH following the TDE, the peak accretion rate can be approximated as

$$\begin{eqnarray}&&\begin{array}{l}{\dot{M}}_{p}\approx \displaystyle \frac{{m}_{\star }/2}{{t}_{\mathrm{acc}}}\\ \quad \approx \ 2.6\times {10}^{2}\,{M}_{\odot }\,{\mathrm{yr}}^{-1}{\left(\displaystyle \frac{h}{0.5}\right)}^{2}\left(\displaystyle \frac{\alpha }{0.1}\right){\left(\displaystyle \frac{{m}_{\star }}{{M}_{\odot }}\right)}^{3/2}{\left(\displaystyle \frac{{R}_{\star }}{{R}_{\odot }}\right)}^{-3/2}.\end{array}\end{eqnarray} \tag{ 20 }$$

The maximum possible luminosity is ${L}_{\max }\sim \epsilon {\dot{M}}_{p}{c}^{2}$, where ∼ 0.1 is the accretion efficiency near the innermost stable circular orbit. This case corresponds to the most efficient energy release (possibly when a jet is formed), and we estimate ${L}_{\max }\sim {10}^{48}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$ for typical TDE parameters. As pointed out in Perets et al. (2016), in this extreme case the TDE may power an ultra-long gamma-ray burst (e.g., Gendre et al. 2013; Levan et al. 2014).

In the more widely accepted adiabatic inflow–outflow solution model (Blandford & Begelman 1999), the mass inflow of a super-Eddington disk onto a BH is nonconservative and only a small fraction of the mass supplied at large radii is actually accreted. The accretion rate is expected to be reduced by factor ${(10{r}_{g}/{r}_{d})}^{s}$, where r_g = Gm_BH/c² is the BH gravitational radius, r_d is the disk radius, and the power-law index s ∈ (0, 1). For the most pessimistic case of s = 1, we can estimate the accretion luminosity of the inner disk, a significant fraction of which should be observable in the X-ray band for favorable viewing angles close to face-on,

$$\begin{eqnarray}&&\begin{array}{l}{L}_{\min }\approx \epsilon \left(\displaystyle \frac{10{r}_{g}}{{r}_{d}}\right){\dot{M}}_{p}{c}^{2}\\ \quad \approx \ 1.5\times {10}^{44}\,\mathrm{erg}\,{{\rm{s}}}^{-1}{\left(\displaystyle \frac{{m}_{\star }}{{M}_{\odot }}\right)}^{7/6}{\left(\displaystyle \frac{{m}_{\mathrm{BH}}}{10{M}_{\odot }}\right)}^{4/3}{\left(\displaystyle \frac{{R}_{\star }}{{R}_{\odot }}\right)}^{-5/2},\end{array}\end{eqnarray} \tag{ 21 }$$

where we have once again assumed h = 0.5 and α = 0.1. For 0 < s < 1, we expect the accretion power to be somewhere in between 10⁴⁴ and 10⁴⁸ erg s⁻¹. For instance, based on numerical simulations of adiabatic accretion flows, Yuan et al. (2012) argued for s ≈ 0.5, which corresponds to L ∼ 10⁴⁶ erg s⁻¹.

At later time t ≫ t_acc, the disk radius increases as r_d ∝ t^2/3 due to viscous spreading, the mass inflow rate drops as $\dot{M}\propto {t}^{-2(s+2)/3}$, and the BH accretion rate drops as ${\dot{M}}_{\mathrm{BH}}\propto {t}^{-4(s+1)/3}$ (e.g., Kumar et al. 2008). Thus, the late-time X-ray light curve is a power law between L_X ∝ t^−4/3 and t^−8/3.

The majority of the mass inflow is lost from the disk in the form of a radiatively driven wind. The energy generated by accretion in the inner disk is expected to be reprocessed by this wind and released at the photon trapping radius, ${r}_{\mathrm{tr}}$, where the radiative diffusion time equals the expansion time (typically occurring a few to 10 days after the TDE; Kremer et al. 2019c). A rough estimate for the trapping radius is

$$\begin{eqnarray}\begin{array}{rcl}{r}_{\mathrm{tr}} & \approx & \min \left({v}_{{\rm{w}}}t,\displaystyle \frac{\dot{M}(t)\kappa }{4\pi c}\right)\\ & & \approx \ {10}^{15}\,\mathrm{cm}\,\min \left(\displaystyle \frac{{v}_{{\rm{w}}}}{{10}^{9}\mathrm{cm}\,{{\rm{s}}}^{-1}}\displaystyle \frac{t}{{10}^{6}\,{\rm{s}}},\displaystyle \frac{\dot{M}(t)}{17{M}_{\odot }\,{\mathrm{yr}}^{-1}}\right),\end{array}\end{eqnarray} \tag{ 22 }$$

where v_w ≈ 10⁹ cm s⁻¹ is the typical wind speed (Kremer et al. 2019c) and κ = 0.34 cm² g⁻¹ is the electron scattering opacity for solar composition material.

For a detailed discussion of the radiation hydrodynamics of the accretion disk and wind, we direct the reader to Kremer et al. (2019c) and Piro & Lu (2020) and references therein. Here, we summarize the key points. As a result of adiabatic loss, the emerging luminosity is smaller than the accretion luminosity by a factor of ${({r}_{d}/{r}_{\mathrm{tr}})}^{2/3}\sim {10}^{-2}$ for ${r}_{\mathrm{tr}}\sim {10}^{15}\,\mathrm{cm}$ and r_d ∼ 10¹² cm roughly at t ∼ 10 days. Kremer et al. (2019c) considered s ∈ (0.2, 0.8) and found the peak bolometric luminosity to be in the range 10⁴¹–10⁴⁴ erg s⁻¹ (depending also upon the assumed stellar parameters, such as masses) and the spectrum to be in the optical/UV, for typical stellar-mass BH TDEs. This optical emission is expected to last roughly 10–100 days until the mass inflow rate drops below roughly (r_d /r_g )L_Edd/c², at which point the disk is expected to transition to a geometrically thin state and the accretion rate may drop by many orders of magnitude (Shen & Matzner 2014).

In Figure 5, we summarize the key evolutionary features of the first roughly 100 days following the tidal disruption.

Zoom In Zoom Out Reset image size

Recent high-cadence surveys have uncovered a growing number of fast-evolving transients with a wide range of observed properties. One class of particular interest is the FBOTs (Drout et al. 2014), also known as fast-evolving luminous transients (Rest et al. 2018). Although a clear understanding of FBOTs remains elusive, this class of transients is generally defined by rise times (of order one to a few days) and peak luminosities (roughly 10⁴¹–10⁴⁴ erg s⁻¹) that are too fast and too luminous to be explained by the radioactive decay of ⁵⁶Ni. The majority of FBOTs are found in star-forming galaxies of roughly solar metallicity. Furthermore, the explosion sites span a range of offsets from the galaxy centers, most closely resembling the offset distribution of core-collapse supernovae (e.g., Drout et al. 2014).

The majority of FBOTs have been identified via archival searches of various optical surveys, including the Pan-STARRS1 Medium Deep Survey (Drout et al. 2014), the Dark Energy Survey (Pursiainen et al. 2018), Kepler (Rest et al. 2018), the Supernova Legacy Survey (Arcavi et al. 2016), and the Zwicky Transient Facility (ZTF; Ho et al. 2020). In addition to the archival searches, a handful of FBOTs have been discovered while still active, notably AT2018cow (also ATLAS18qqn; Prentice et al. 2018; Rivera Sandoval et al. 2018; Smartt et al. 2018; Margutti et al. 2019), ZTF18abvkwla (Ho et al. 2020), and CSS161010 (Coppejans et al. 2020), enabling X-ray and/or radio follow-up.

Various analyses have computed volumetric rates of transients similar to FBOTs, with estimates ranging from roughly 100 Gpc⁻³ yr⁻¹ to more than 1000 Gpc⁻³ yr⁻¹ in the local universe (e.g., Drout et al. 2014; Pursiainen et al. 2018; Coppejans et al. 2020; Ho et al. 2020). Although the precise rate remains uncertain, the general consensus appears to be that FBOTs are roughly two to three orders of magnitude rarer than standard core-collapse supernovae (e.g., Botticella et al. 2008).

The specific origin of FBOTs remains unknown, with a number of channels having been proposed, including TDEs by IMBHs (Perley et al. 2019), massive star collapse and BH formation (Quataert et al. 2019), electron capture collapse following a white dwarf merger (Lyutikov & Toonen 2019), and magnetar formation (Margutti et al. 2019). Additionally, Piro & Lu (2020) pointed out that many FBOT features (specifically in the case of AT2018cow) show similarities to what would be expected for a wind-reprocessed transient, regardless of the specific central engine.

Here, we propose stellar-mass BH TDEs as another possible FBOT progenitor. The rise times and peak optical luminosities predicted for these TDEs (Section 4.1 and Kremer et al. (2019c)) occupy the same region of parameter space expected for FBOTs, as does the estimated X-ray luminosity. ⁸ If these TDEs occur in high-metallicity young stellar clusters (expected to be a dominant site for star formation; Lada & Lada 2003), the host-galaxy type and offset distribution for observed FBOTs may also be recovered. Furthermore, the rate we predict for stellar-mass BH TDEs in YSCs (up to roughly 200 Gpc⁻¹ yr⁻¹ in the local universe) is comparable to the observationally inferred FBOT rates. In light of these similarities, stellar-mass BH TDEs that occur in YSCs may in principle be a viable progenitor for FBOT-like transients. In Table 3, we summarize for comparison the key features of both FBOTs and stellar-mass BH TDEs.

Although it remains to be determined whether bright radio emission is a defining feature of the FBOT class broadly, radio emission (consistent with self-absorbed synchrotron radiation) is observed for the AT2018cow, ZTF18abvkwla, and CSS161010 events. This radio emission suggests the presence of circumstellar medium with densities ranging from roughly 10–10⁶ cm⁻³, for various observation epochs and for a range of microphysics assumptions (Margutti et al. 2019; Coppejans et al. 2020; Ho et al. 2020). If prior to its tidal disruption, a main-sequence star is tidally captured by a BH, such that a small amount of orbital energy is deposited into the star (e.g., Fabian et al. 1975, see left-most panel of Figure 5 for illustration), a small amount of debris may be unbound from the star. This may occur either through partial stripping by the BH at the first pericenter passage or during subsequent pericenter passages if the star's envelope expands due to energy injected through the tidal encounter. The circumstellar medium density inferred from radio emission may be produced by a small amount of ejecta (≲10⁻² M_⊙) homologously expanding for a few years at a speed of 10³ km s⁻¹. More detailed hydrodynamic models of the disruption process are required to calculate the radial density profile and predict the radio emission.

We summarize here the main findings of this work:

• 1.  
  Using a suite of roughly 3000 N-body simulations, we have explored the rates and properties of stellar-mass BH TDEs in YSCs. We derived TDE rates as a function of cluster mass and showed that the rates derived from the models agree closely with rates derived from simple analytic estimates.
• 2.  
  Using the rate scalings derived from our models and integrating over the full cluster mass function, we predict these TDEs occur at a rate of roughly 20–200 Gpc⁻³ yr⁻¹ in the local universe. The range in this estimate is determined primarily by the assumptions made concerning initial cluster densities (e.g., Portegies Zwart et al. 2010; Marks & Kroupa 2012). Overall, we find TDEs occur a factor of a few to 10 times more frequently in binary-mediated dynamical encounters compared to single–single encounters.
• 3.  
  In YSCs of roughly solar metallicity, we predict a median BH mass of roughly 10 M_⊙. In lower-metallicity (Z ≲ 0.1 Z_⊙) clusters that may have formed at higher redshift, we predict a median mass of roughly 30 M_⊙. This difference is primarily a result of the metallicity-dependent stellar wind mass loss (Vink et al. 2001).
• 4.  
  We showed that stellar-mass BH TDEs exhibit the following key features: (i) rise times of roughly a day driven by the viscous accretion timescale, (ii) peak X-ray luminosities of roughly 10⁴⁴–10⁴⁸ erg s⁻¹, and (iii) optical luminosities of up to roughly 10⁴⁴ erg s⁻¹ produced by reprocessing of X-rays by a disk wind on a timescale of a few to 10 days after the TDE. On the basis of all of these features, combined with the estimated TDE rates and host-galaxy properties, we propose stellar-mass BH TDEs as a viable progenitor for the FBOT class of transients (e.g., Drout et al. 2014; Margutti et al. 2019; Coppejans et al. 2020; Ho et al. 2020).

Stellar-mass BH TDEs distinguish themselves from other proposed channels for FBOTs (see Section 4.2) in the combination of X-ray and optical signatures. For instance, in the case of stellar collapse, long-lasting power from a central engine may be supplied by fallback accretion onto a BH or magnetar spindown, where the late-time energy injection rates are expected to be L ∝ t^−5/3 and L ∝ t⁻², respectively. In our model for stellar-mass BH TDEs, the energy injection rate may have a broad range of energy injection rates from L ∝ t^−4/3–L ∝ t^−8/3. This leads to different late-time (t ≳ 10 days) behaviors in the light curve that can in principle to used to distinguish our model from others. Of course, the detection of a host YSC may also hint at the proposed TDE scenario (for example, see Lyman et al. 2020, for a discussion of host features of AT2018cow). However, core-collapse supernovae will also occur in sufficiently YSCs, thus, the presence of a host cluster may not provide indisputable evidence of a TDE origin.

Although the event rates inferred from FBOT observations are uncertain and vary between different analyses, in general, the observed FBOT rate appears higher than the stellar-mass BH TDE rate estimated from our models by a small factor (see Table 3). However, there are several reasons why our estimated rate may underestimate the true TDE rate, possibly explaining this discrepancy:

First, in this study, we have assumed zero primordial binaries. BH binaries were allowed to form exclusively through three-body formation (see Section 2). Stellar binaries are known to increase the rates of dynamical collisions and TDEs (e.g., Fregeau & Rasio 2007). In this case, our predicted TDE rates may be a lower limit. To test this, we ran 50 additional simulations (o in Table 1) that adopt a primordial binary fraction of 100%. Consistent with previous CMC studies (e.g., Kremer et al. 2020b), secondary masses are drawn from a uniform distribution in mass ratio in the range of [0.1, 1] and initial orbital periods are drawn from a log-uniform distribution (e.g., Sana et al. 2012). In total, these models yield a factor of roughly 5 more TDEs per cluster (all of which occur through binary-mediated encounters) compared to the models with comparable mass and zero primordial binaries (f in Table 1). Thus, primordial binaries may increase the TDE rates quoted in Table 2 by a small factor. We reserve for follow-up studies a more expansive investigation of primordial binaries and their implications for binary-mediated TDEs.

Second, as discussed briefly in Section 4, a fraction of TDEs may occur through tidal capture where the main-sequence star may undergo multiple passages before ultimately being disrupted. As discussed in, e.g., Fabian et al. (1975), the cross section for tidal capture may be a factor of roughly a few times larger than that for tidal disruption. If the fate of the majority of tidal captures is a TDE, the total TDE rate may be higher than that estimated here by a factor of a few.

Of course, realistic YSCs exhibit a much larger range in properties than those considered here. For instance, in our N-body models, we have adopted a relatively narrow range in initial cluster sizes, which may in reality extend from roughly 0.1–10 pc or even more, depending on the cluster mass (e.g., Krumholz et al. 2019). Although our models were intended to capture the "mean" of this distribution, this value is uncertain, and of course, lower-density clusters would exhibit a lower TDE rate. For example, from the supplemental models n in Table 1, we find that the TDE rate for clusters with initial r_v = 2 pc is a factor of roughly 6 times lower than than the rates computed in the r_v = 1 pc models. If indeed r_v = 2 pc is more representative of the full distribution of YSCs, the volumetric rates shown in Table 2 may overestimate the true value.

Along similar lines, we have adopted a relatively narrow range in initial cluster masses that does not span the full distribution of observed YSCs. For example, YSCs in the antennae are observed with masses of 10⁶ M_⊙ (e.g., Portegies Zwart et al. 2010) or more, while observed clusters in massive elliptical galaxies can reach 10⁷ M_⊙ (e.g., Bastian et al. 2006). Although of high potential interest (see, e.g., Rodriguez et al. 2020, for a recent study of super-star clusters with CMC), we reserve consideration of these high-mass systems for later study.

In Section 3.4, we showed that, depending on the assumed initial virial radius and cluster metallicity, IMBHs with masses in excess of 100 M_⊙ may form in YSCs through stellar collisions (see also Di Carlo et al. 2019a; Kremer et al. 2020a) and ultimately undergo TDEs. The topic of TDEs by IMBHs in stellar clusters has been examined at length (e.g., Rosswog et al. 2009; MacLeod et al. 2014, 2016; Fragione et al. 2020). We have focused in this study on the more common (by a factor of roughly 10:1 or higher in our models) stellar-mass BH TDE. However, IMBH TDEs likely lead to several notable differences. For example, for sufficiently massive IMBHs, the TDEs will transition from viscosity driven to fallback driven, thus becoming qualitatively more similar to an SMBH TDE (e.g., Rees 1988). We leave for future work a more careful examination of the properties and observational signatures of TDEs by IMBHs.

We have used here the CMC code to model the evolution of YSCs. The Monte Carlo-based approach used in CMC allows us to model large populations of clusters at low computational expense compared to direct N-body models (e.g., Pattabiraman et al. 2013; Rodriguez et al. 2016). However, unlike massive globular clusters, the relatively low-mass YSCs studied here can also be studied efficiently with direct N-body models (e.g., Banerjee 2017; Di Carlo et al. 2019a). Future work may examine the topic of TDEs using direct N-body models, which would allow detailed examination of several regimes that Monte Carlo simulations like CMC are ill-equipped to study. For example, examination of the final stage of cluster's life as it dissolves on a dynamical timescale through its tidal boundary and examination of clusters that contain massive IMBHs that may affect the overall cluster and TDE dynamics. Additionally, direct N-body models are more suited to study the lowest-mass stellar associations (N ≲ 100; Lada & Lada 2003), for which the spherical symmetry assumptions at the heart of the Monte Carlo-based approach break down.

Finally, unlike SMBH TDEs, which have been studied extensively with hydrodynamic simulations, the stellar-mass BH regime has been little explored, with a few exceptions (Perets et al. 2016; Lopez et al. 2020). Ultimately, hydrodynamic models are necessary to understand the detailed features of these events. For example, in the analytic estimates in Section 4, we considered only those interactions near the tidal disruption boundary, r_TD, of the star. However, given that the interaction cross section scales linearly with stellar radius (Equation (4)), a fraction of such BH–star encounters may actually occur in the direct collision regime, where r_p < R_⋆ (e.g., Fryer & Woosley 1998; Hansen & Murali 1998). Although this subclass of "head-on" encounters is unlikely to affect our rate predictions by more than a small factor, these physical collisions may produce a subclass of unique transients. For instance, a direct collision may lead to prompt accretion (i.e., t_fb ≈ 0) onto the BH, which, among other possible consequences, may make the effects of feedback critical on the subsequent evolution, if indeed a disk forms at all. On the other hand, as mentioned briefly in Section 4, even more distant encounters (r_p > r_TD) may lead to tidal capture, possibly resulting in multiple passages that each partially strip the star (e.g., Fabian et al. 1975; Ivanova et al. 2017), again potentially producing EM signatures unique from those presented in the classic TDE regime discussed in Section 4. More careful treatment of the various regimes of BH–star interactions with hydrodynamic models is necessary to explore the potentially broad range of outcomes of these events in greater detail.