Searching for Intermediate-mass Black Holes in Globular Clusters through Tidal Disruption Events

Intermediate-mass black holes (IMBHs) may be the link between stellar mass holes and the supermassive variety in the nuclei of galaxies, and globular clusters (GCs) may be one of the most promising environments for their formation. Here, we carry out a pilot study of the observability of tidal disruption events (TDEs) from 103 M ⊙ < M • < 105 M ⊙ IMBHs embedded in stellar cusps at the center of GCs. We model the long super-Eddington accretion phase and ensuing optical flare, and derive the disruption rate of main-sequence stars as a function of black hole mass and GC properties with the help of a 1D Fokker–Planck approach. The photospheric emission of the adiabatically expanding outflow dominates the observable radiation and peaks in the near-ultraviolet/optical bands, outshining the brightness of the (old) stellar population of GCs in Virgo for a period of months to years. A search for TDE events in a sample of nearly 4000 GCs observed at multiple epochs by the Next Generation Virgo Cluster Survey yields null results. Given our model predictions, this sample is too small to set stringent constraints on the present-day occupation fraction of GCs hosting IMBHs. Naturally, better simulations of the properties of the cluster central stellar distribution, TDE light curves, and rates, together with larger surveys of GCs are all needed to gain deeper insights into the presence of IMBHs in GCs.


INTRODUCTION
There is currently no unambiguous confirmation for intermediate mass black holes (IMBHs) with masses 10 2 M ⊙ ≲ M • ≲ 10 5 M ⊙ (for a recent review, see Greene et al. 2020).If one were to extrapolate the known M • −σ ⋆ relation down to low-mass stellar systems, globular clusters (GCs) would be expected to hosts IMBHs with M • ∼ 10 3 − 10 4 M ⊙ (Lützgendorf et al. 2013).Finding objects and characterizing their mass function in this range would provide unique insight into the nature and growth of massive black hole seeds in the early Universe and the dynamical evolution of dense stellar systems, and provide key input into event predictions for gravitational wave facilities.Thus far there is some circumstantial evidence for black holes below 10 5 M ⊙ in galaxy nuclei (NGC 205, Nguyen et al. 2019) and in hyperluminous X-ray sources (HLX-1, Webb et al. 2012), but there are still no compelling candidates with M • ∼ 10 3 M ⊙ .
Different mechanisms have been proposed for the formation of IMBHs in dense stellar clusters, from "slow" scenarios where the black hole remnant of a massive star sinks to the center of the cluster and grows over time through mergers of mass-segregated lighter black holes (Miller & Hamilton 2002), to "fast" collisional runaway of massive stars during an early phase of cluster core collapse, the product of which may eventually lead to the formation of an IMBH (Portegies Zwart & McMillan 2002;Giersz et al. 2015).It has been shown by Miller & Davies (2012) that all clusters above a central velocity dispersion of σ ⋆ ∼ 40 km s −1 will necessarily form an IMBH through some mechanism at any cosmic epoch, since above this dispersion primordial binaries cannot support the system against deep core collapse.The case is more complicated for lower-dispersion systems such as GCs, where -regardless of how IMBHs form -it may be a challenge to retain them in the face of a repeated onslaught of gravitational wave kicks from mergers with other black holes (Holley-Bockelmann et al. 2008;Fragione et al. 2018a).
To make progress, it is key to develop alternative methods for IMBH searches in dense stellar systems.Tidal disruption events (TDEs) may offer an independent, promising probe.A TDE occurs when the tidal force of the black hole exceeds the star self-gravity and rips it apart (Rees 1988).While about half of the stellar debris is ejected at high speed, the remainder gets accreted, producing an optical/UV flare accompanied by thermal X-ray emission.One of the most compelling IMBH TDE candidates to date is 3XMM J215022.4-055108(Lin et al. 2018), an X-ray outburst with a luminosity that peaked at 10 43 erg s −1 and decayed systematically over 10 years, hosted in a star cluster of mass ∼ 10 7 M ⊙ and plausibly powered by an IMBH of mass Wen et al. 2021).This event suggests that an effective means of detecting IMBHs may be through a search of optical flares in a large sample of dense star clusters (e.g., Ramirez-Ruiz & Rosswog 2009;Chen & Shen 2018;Fragione et al. 2018b).As we shall see, this is particularly true for searches in the NUV, where a disruption event is predicted to outshine the old stellar population of the average GC for a period of months to years.
Large surveys with precise multi-wavelength photometry are becoming increasingly available for public use.As a pilot experiment, we take advantage of the Next Generation Virgo Cluster Survey (NGVS, Ferrarese et al. 2012), a comprehensive optical imaging survey of the Virgo galaxy cluster, together with its near-infrared counterpart (NGVS-IR, Muñoz et al. 2014), to set constraints on the fraction of GCs that may be harboring IMBHs.The NGVS-IR covers a total area of 4 deg 2 centered on Virgo's core region in the u * grizK s bandpasses, and we capitalize on the sample of nearly 4,000 Virgo GCs provided by Peng et al. (in prep.).We use theoretical optical light curves to search for potential TDE candidates, and the estimated event rates of TDEs for a given GC population, stellar density profile, and cadence of NGVS observations, to attempt to constrain the fraction of GCs hosting IMBHs with stellar cusps at the present epoch.While there are considerable uncertainties in the modeling of TDE events around IMBHs in GCs, we hope that our pilot investigation will provide a blueprint for future searches of optical flares in dense stellar systems.

LIGHT CURVE MODELING
While the evolution of a TDE light curve has been extensively modeled in the literature (e.g., Strubbe & Quataert 2009;Lodato & Rossi 2011;Strubbe & Quataert 2011;Metzger & Stone 2016;Mockler et al. 2019;Ryu et al. 2020;Lu & Bonnerot 2020, and references therein), the details of the evolution of the debris stream, the efficiency of the process of circularization, and the emission mechanisms responsible for the optical/UV emission in TDE flares all remain an open question (for recent reviews see Roth et al. 2020;Rossi et al. 2021).For illustrative purposes, we shall follow here the simplified analysis of Strubbe & Quataert (2009) and Lodato & Rossi (2011).A star of mass M ⋆ = m ⋆ M ⊙ and radius R ⋆ = r ⋆ R ⊙ on a very eccentric orbit is torn apart when its pericenter r p is within the black hole's tidal sphere of radius For lower main-sequence stars, masses and radii are related by r ⋆ = m 0.8 ⋆ (Kippenhahn & Weigert 1990).The impact parameter of the encounter, β = r t /r p , measures the strength of the tidal interation.Initially, approximately half of the initial stellar debris becomes bounded (Evans & Kochanek 1989), and after a fallback time, most bound material starts coming back to pericenter at the rate While the early behavior of the fallback rate may be influenced by stellar properties (Lodato et al. 2009;Guillochon & Ramirez-Ruiz 2013;Law-Smith et al. 2020), the late-time accretion rate onto the black hole from a TDE always declines as t −5/3 if the star is completely disrupted (Coughlin & Nixon 2019).
For an IMBH, the mass fallback rate predicted by Equations ( 2) and (3) exceeds the Eddington rate for a timescale of years or longer.Here, κ es = 0.34 cm 2 g −1 is the electron scattering opacity, and we assume accretion process at 10% radiative efficiency.In this regime, only a fraction of the fallback material joins the newly-formed Keplerian disk and accretes onto the black hole at the rate while the remaining stellar debris are launched in a radiation pressure-driven outflow during circularization.The effective temperature of the accretion flow is that of a slim advective disk (Strubbe & Quataert 2009), where f ≡ 1− 3R S /R and R S = 2GM • /c 2 is the Schwarzschild radius.The disk extends from 3R S to the circularization radius R c = 2r p = 2r t /β, and the monochromatic disk luminosity as a function of time is that of a multicolour blackbody, When required, our standard choice of parameters for the stellar properties is β = 1, m ⋆ = 1, and r ⋆ = 1. Figure 1 shows the bolometric light curve for the disk emission from the disruption of solar-type star by a 3,000 M ⊙ IMBH.Note how, during the super-Eddington phase, viscosity-generated heat does not have sufficient time to be radiated away, and is instead advected into the hole.The disk effective temperature remains then constant with time even as the fallback rate declines, and the radiated luminosity saturates at a few times L Edd .
A simple model for the radiation-driven wind approximates the geometry as spherical and assumes that stellar debris falls back at close to the escape velocity v esc = 2GM • /R c and shocks at R c , converting kinetic energy into radiation.Photons are trapped in an outflow promptly launched from the circularization radius with Ṁout (t) = f out Ṁfb and terminal velocity v w = f v v esc .
Our standard choice for the wind velocity parameter is The radiation temperature at the base of the wind, T c , can be derived from energy conservation in the wind 4πR The outer radius of the ejecta grows as , and the outflow is optically thick to electron scattering out to a photospheric radius R ph given by R ph ρ ph κ es = 1.The gas density profile for r < R ej follows from matter conservation, ρ(r, t) = Ṁout (t)/(4πr 2 v w ).The temperature of the advected radiation decreases adiabatically as T ∝ ρ 1/3 , and photons escape with a blackbody spectrum at the photospheric temperature of where ρ c = Ṁfb /(4πR 2 c v esc ) is the gas density at the base of the flow.
At the earliest times, the fallback rate can be so large and the density so high that the location of the photosphere is just inside R ej .In this "edge-limited" phase, we set R ph = R ej and ρ ph = (R ej κ es ) −1 .The specific outflow luminosity is given by Figure 1 shows the resulting wind bolometric light curve, which rises as t 11/9 during the edge-limited initial phase, when the photosphere expands and T ph decreases as t −7/36 .As time passes and the density of the outflow subsides, the photosphere sinks inward as t −5/3 , its temperature rises as t 25/36 , and the luminosity declines as t −5/9 .The total radiation luminosity of the wind is of order L Edd in the f out = 0.9 case, decreasing for lower outflow rates as f out during the edge-limited phase, and as f 5/3 out afterwards (cf.Strubbe & Quataert 2009).These relations only apply for R ph ≥ R c , because otherwise the outflow is optically thin.We impose this limit by setting the wind luminosity to zero when the photosphere shrinks below R c .In the f out = 0.1 case the wind becomes optically thin when the fallback rate is still super-Eddington.
While the disk always dominates the total radiated power, the bulk of the disk emission occurs in the soft X-rays.The optical/NUV flash plotted in Figure 2 is instead produced by the adiabatically expanding outflow, which largely outshines the disk during the optically thick phase.For comparison, we also show in the figure the mean NUV luminosity of GCs detected in the NGVS survey, which argues for the detectability of a TDE optical transient against the brightness of a dense star cluster.
Depending on the uncertain geometry of the outflow and the viewing angle of the observer to the source, it is possible that both the accretion disk and the outflow may be visible at early times.Here, we shall neglect obscuration effects and sum up the two contributions to produce a total TDE light curve.

TDE RATES: ESTIMATES
To provide estimates and trends of the number of TDE flares expected in a local GC survey like the NGVS, we present in this section a simple and intuitive configuration-space approach following Syer & Ulmer (1999).Our treatment is informed by the more rigorous numerical integration of the one-dimensional, orbit-averaged Fokker-Planck equation for the evolution of the stellar distribution function, which we perform and discuss in the next section.
Stars within a cluster of current mass M GC are assumed to follow a Plummer density profile with 1D velocity dispersion Here, ρ c = 3M GC /(4πa 3 ) and σ 2 c = GM GC /(6a) are the density and velocity dispersion in the cluster core, and a is the Plummer scale parameter that sets the size of the cluster core.We have fit the masses and structural parameter data of Milky Way's GCs of Baumgardt & Hilker (2018) with the core density-cluster mass relation: For a GC with mass M GC = 5×10 5 M ⊙ , the expressions above yield ρ c = 4.3 × 10 4 M ⊙ pc −3 , a = 1.4 pc, and σ c = 16 km s −1 .At the center of this Plummer cluster we place an IMBH with M • ≪ M GC .Following mass segregation (Bahcall & Wolf 1977), low mass main-sequence stars and their remnants relax to a Bahcall-Wolf r −γ density cusp near the black hole, with 3/2 < γ < 7/4 (Alexander & Hopman 2009).We thus write the stellar density near the center of the cluster as where R i ≪ a is the radius of influence of the black hole.This is defined as the location where the interior stellar mass equals M • (e.g., Binney & Tremaine 1987;Merritt 2004;Wang & Merritt 2004), i.e.
where g(γ) = (3 − γ)/(6 − γ).In our set-up, R i is typically a few times larger than the radius r • ≡ GM • /σ 2 c defined kinematically in terms of the intrinsic stellar velocity dispersion.For IMBHs weighing a few percent of the cluster core mass, collisional N -body simulations show that the ρ ⋆ ∝ r −γ cusp extends all the way up to R i > r • (Baumgardt et al. 2004).This is also confirmed in our direct numerical integration of the Fokker-Plank equation for the time-dependent stellar distribution in the vicinity of an IMBH (see next Section).The Jeans equations associate an r −γ density cusp with a velocity dispersion that near the black hole approaches ] (e.g., Zhao 1996).We therefore approximate the cluster central velocity dispersion as Hereafter, we shall assume a slope of γ = 7/4 for the central stellar cusp.
Stars on nearly parabolic orbits are tidally disrupted when their specific angular momentum J is smaller than J lc = (2GM • r t ) 1/2 .The ensemble of nearly-radial orbits with J < J lc forms the so-called "loss cone", the set of velocity vectors at some distance r from the black hole that lie within a cone of half-angle θ lc (Frank & Rees 1976;Lightman & Shapiro 1977) Stars in the loss cone are disrupted at the first periapsis passage, and the continued supply of stars to the IMBH is driven by gravitational deflections that repopulate the loss cone on a timescale that is related to the two-body relaxation time 1987).Here, n ⋆ = ρ ⋆ /⟨M ⋆ ⟩ is the number density of main sequence stars and their remnants, and we take ln Λ = ln(0.4M • /⟨M ⋆ ⟩) for the Coulomb logarithm.In a simplified treatment, two dynamical regimes can be defined depending on the ratio of the dynamical timescale, t d = r/σ ⋆ , to the diffusion timescale of angular momentum across the loss cone, t lc = θ 2 lc t R .Close to the IMBH, orbital periods are short and stars diffusing into the loss cone are immediately disrupted.This is the empty loss-cone regime, and the rate of TDEs per star can be estimated as where N MS (< r) is the number of main-sequence stars contained within a radius r (Frank & Rees 1976;Lightman & Shapiro 1977;Syer & Ulmer 1999).This flux is proportional to the relaxation rate but only weakly dependent on the size of the loss cone.In the other regime, t d > t lc , scattering in and out of the loss cone is faster than the orbital time, and the loss cone will always be full, uniformly populated in orbital phase.The fraction of stars in the loss cone at any time is then just θ 2 lc , and the TDE rate per star is given by (Syer & Ulmer 1999) independent of the stellar encounter timescale.After solving for the transition radius r crit where the two perstar disruption rates in Equations ( 19) and ( 20) are equal, we write the total TDE rate of main-sequence stars as In the above equations, the term ⟨M ⋆ ⟩ represents the change in the event rate with the number of stars at fixed total stellar density ρ ⋆ , while the second moment ⟨M 2 ⋆ ⟩ is responsible for the decrease of the diffusion timescale as the gravitational potential becomes more "granular" and stellar-mass black holes dominate the relaxation rates (Kochanek 2016;Bortolas 2022).In general, to compute ⟨M ⋆ ⟩ and ⟨M 2 ⋆ ⟩ one needs to adopt a presentday mass function (PDMF) with an upper truncation at M max ⋆ ≈ 1 M ⊙ to approximate the old stellar population of a GC (Magorrian & Tremaine 1999).For simplicity, and in analogy with the time-dependent Fokker-Planck numerical approach discussed in the next section, we will assume instead an idealized, single-mass stellar system composed by M ⋆ = 1 M ⊙ stars (n MS = n ⋆ ).Compared to a single-mass stellar population, a realistic PDMF is known to enhance TDE rates by several-fold (Stone & Metzger 2016) even after accounting for the smaller tidal radii, , of sub-solar stars. 1 Figure 3 shows the predicted TDE rates for 10 3 < M • /M ⊙ < 104 IMBHs in GCs with a single-mass stellar 1 In our model, the empty loss-cone disruption rate can be shown to scale as Ṅ< ∝ r 4/9 t up to a slowly varying logarithmic factor.population and increasing core densities in the range ρ c = 10 4.5 − 10 6 M ⊙ pc −3 .Most events are sourced at r crit < R i , and both the empty and full loss-cone regimes contribute an O(1) fraction of the total disruption rate.The rates of TDEs increase with black hole mass and cluster core density approximately as ṄTDE ∼ (ρ c M • ) 0.5 , and typically exceed 10 −5 yr −1 for ∼ > 10 3 M ⊙ IMBHs harbored in dense GCs.Obviously, these high rates cannot be sustained indefinitely because of stellar depletion.In the next section we will compute TDE rates in a time-evolving stellar density profileincluding depletion -using a Fokker-Plank approach.
We have tested the simplified model of the previous section against the results of the PhaseFlow Fokker-Plank integrator of Vasiliev (2017).The publicly available PhaseFlow code solves the coupled Poisson and orbit-averaged 1D Fokker-Planck equations, evolving a spherically symmetric distribution of stars under twobody relaxation and loss-cone effects in the neighbourhood of a central black hole, the mass of which is allowed to grow with time following stellar captures (Bortolas 2022).We explore a one-parameter family of models, varying the initial core density of the initial Plummer profile between 10 4.5 − 10 6 M ⊙ pc −3 and adjusting the total mass of the cluster according to Equation ( 13) (M GC = 4.3 × 10 5 − 2.4 × 10 6 M ⊙ ).In all configurations, a central IMBH of initial mass M • (t = 0) = 600 M ⊙ is embedded in a r −7/4 initial stellar cups that extends all the way up to the influence radius R i .Following each capture event, 30% of the mass of the disrupted star is accreted by the IMBH.For simplicity, we assume here an idealized, single-mass stellar system composed by M ⋆ = 1 M ⊙ stars, postponing a numerical treatment that includes a complete, time-dependent stellar mass function and the effect of mass segregation on TDE rates to future work.
Figure 4 shows the evolving stellar density profiles of four IMBH+GC systems having different initial core densities in the range ρ c = 10 4.5 − 10 6 M ⊙ pc −3 (as in Figure 3).As shown by Vasiliev (2017), the cusp is not a steady-state structure: the gray line shows the initial Plummer+Bahcall-Wolf cusp profile (t = 0), the red line corresponds to the time when the cusp amplitude attains its maximum value, and the blue line depicts the profile at the end of the integration, t = 14 Gyr, when the cusp density normalization has decreased in response to the heat source at the center.In the figure, the vertical lines mark the influence radius R i containing a mass of stars equal to M • (t), while the starred points denote the locus where most TDEs are sourced.
The evolution of all clusters follows a similar route.At first the system expands, powered by an outward heat conductive flux driven by two-body relaxation.The density in the cusp decreases, roughly maintaining an r −7/4 profile inside the influence radius.The initial expansion phase lasts ∼ 0.1 Gyr, which is of order the two-body relaxation timescale at the influence radius R i (t = 0), and is followed by a phase of secular gravo-thermal contraction leading to a maximum cusp density at time t peak ≃ 8 Gyr; as expected, the latter timescale is of order the relaxation timescale measured in the outer regions of the cluster.IMBH-star interactions eventually generate enough heat to prevent further core collapse and cause a late re-expansion of the cluster.During the 14 Gyr evolution, the IMBH acquires a mass of 2, 919 − 48, 720 M ⊙ depending on the chosen initial stellar density, with a more efficient growth in the case of denser and more massive systems (Figure 5) and a faster rise at time t = 7 − 9 Gyr during the gravo-thermal cluster contraction.The three main phases of the GC density profile progression -fast expansion, contraction of the core, and slow re-expansion -are reflected in the evolution of the TDE rate (Figure 6), which is simultaneously modulated by the growth of the central black hole.Disruption rates reach a minimum after 0.1 − 0.2 Gyr, then climb dramatically by more than 2 orders of magnitude during core contraction, peak at t peak ≃ 8 Gyr, and slowly decay at late times.The peak in the rates is clearly sourced by the collapse of the outer regions of the cluster delivering a substantial amount of new stars near the IMBH and promoting its growth.Peak rates are more pronounced in initially denser clusters where stellar captures lead to more massive IMBHs.The orbital radius where most disruption events originate (starred point in Figure 4) moves outward during the core expansion phases, and inward when the core contracts.
Because of the intrinsic time-dependence of the problem, a detailed comparison between our steady-state analytical estimates in Section 3 and the numerical results obtained with the PhaseFlow code is far from trivial.The Fokker-Planck integration validates qualitatively: • the assumed steady-state inner density profile, which actually evolves in amplitude in Figure 4 but approximately retains the same r −7/4 functional form, with little dependence on initial conditions; • the IMBH mass range used to estimate light curves and TDE rates in GCs, which is consistent with the mass accreted in ∼ 7 − 8 Gyr by sufficiently massive black hole seeds (with seed masses above a few hundred solar masses) via stellar captures in dense stellar systems (see Figure5); and • the trends of increasing TDE rates with M • and ρ c (Figure 3), which are also seen in our numerical integration model (Figure 6).
When compared with our steady-state analytical estimates for a single-mass stellar population of 1 M ⊙ stars (Figure 3), the peak rates observed in the Fokker-Planck integration appear to be a few times lower than the corresponding rates of the simplified model.
In our analysis below of the NGVS GC sample, we shall adopt for simplicity the analytical rates of Figure 3, as they are expected to be intermediate between the numerical results and steady-state estimates that include a complete, evolved stellar mass function.A useful powerlaw fit to these TDE rates as a function of M • and ρ c is given by with parameters a = (2.2 ± 0.3) × 10 −6 yr −1 , b = 0.44 ± 0.002, and c = 0.71 ± 0.005.2

NGVS GC SAMPLE
To compare our theoretical predictions with observations, we search for optical flares from TDEs in a robust sample of GCs in u-and g-band.These clusters are in the core of the nearby Virgo galaxy cluster within a 2 deg 2 region centered on M87.The photometric data are part of the NGVS (Ferrarese et al. 2012) and its near-infrared counterpart, NGVS-IR (Muñoz et al. 2014).GC candidates were selected using extreme deconvolution (Bovy et al. 2011) to model the distribution of foreground stars, GCs, and background galaxies in a multidimensional parameter space of color and morphology (concentration parameter), which determined the probability of a given source to be a GC (Peng et al., in prep.).Our sample of NGVS globulars consists of CFHT MegaCam time-series photometry in the ugiz bandpasses spanning 5 years with a cadence timescale ranging from hours to months and sometimes several years.The cadence was determined by the need to make large dithers for studying the low surface brightness outskirts of massive galaxies (Ferrarese et al. 2012).Most GCs have about 5-20 measurements in each of the four filters with 0.05−0.2mag photometric precision depending on the apparent magnitude of the system.After removing those with only single-epoch observations, our final NGVS sample of GCs with u-band photometry includes a total of 3,849 sources.

Photometric Variability
A temporary brightening of a NGVS GC could indicate a potential TDE, a false positive due to cosmic ray hits, or other image artifacts.We have checked for large photometric variations, applying a sigma-clipping technique to find outliers in subsets of data.Four outliers  (Jordán et al. 2007).The vertical red dashed line marks the mean value of the distribution, ⟨MGC⟩ = 1.1 × 10 6 M⊙.In the u-band, most GCs are fainter than 10 39 erg s −1 .at 3.5σ were found in the u−band, and image inspection showed that cosmic rays were causing the increase in brightness.Our search for detectable TDE events in the NGVS GC sample yields null results.

Mass and Luminosities
Stellar masses of NGVS GCs were computed from their AB magnitudes m AB as The constant C includes the logarithm of a mass-to-light ratio Υ in any arbitrary band-pass and the distance d to the Virgo cluster [16.5 ±0.1 (random) ±1.1 (systematic) Mpc, Mei et al. 2007], and is given by From the reported MegaCam photometry, we take the average g-and z-band magnitudes and compute for each GC the color (g − z).The colors are then used to obtain z-band mass-to-light ratios by adopting the corresponding values predicted by the metallicity-dependent population-synthesis model PEGASE (Jordán et al. 2007).The resulting NGVS GC mass distribution is shown in Figure 7.The GC sample has mean stellar mass ⟨M GC ⟩ = 1.1 × 10 6 M ⊙ , and 31% of all globulars are more massive than 10 6 M ⊙ .We used z-band rather than g−band for mass determinations because mass-tolight ratios at redder wavelengths are less sensitive to the exact age or star formation history of a stellar population, minimizing systematic uncertainties.Figure 7 also shows the u-and g-band monochromatic luminosities of our GC sample.Globulars are typically fainter in u, with NUV luminosities rarely exceeding 10 39 erg s −1 .

Detectability of TDEs in NGVS GCs
To assess quantitatively the observability of TDEs in our GC sample, and in preparation for future multiepoch surveys of extra-galactic star clusters, we assume here that each NGVS GC hosts an IMBH of mass 0.01 M GC .The input parameters for each cluster are its stellar luminosity, the corresponding TDE rate (Eq.22), and the theoretical light curves in the two cases, f out = 0.1 and f out = 0.9, chosen to bracket the strength of the radiation-driven wind component.For each GC we then compute the time interval ∆t TDE during which the TDE outshines the steady stellar GC component in the u-band, as well as the total time span between the first and last epoch of cluster observation, ∆t S .The number of observable TDEs in the sample is then Here, the sum is over all GCs, f occ < 1 is the occupation fraction, and the first term on the RHS corrects for the events that went off before the start of the survey but whose light-curve tail would still be visible at survey times.Note that typically ∆t TDE > ∆t S , and that ∆t S can be as long as 5.1 yr.Our model predicts N TDE /f occ = (10, 4) events for f out = (0.1, 0.9), respectively.Monte Carlo simulations of TDEs occurring randomly at the rate ṄTDE,i and observed for a timescale ∆t S at the cadence of the NGVS confirm these basic estimates.Figure 8 shows the distribution of ∆t TDE for our GC sample at varying f out , with the weaker winds+disk systems outshining the steady stellar component for longer timescales because of late time accretion.Figure 9 shows three examples of TDE mock detections in the Monte Carlo simulations, with the TDE light curve (the sum of the disk and wind components) compared to the steady stellar luminosity in the u-band.In the figure, the total time span between the first and last epoch of cluster observation is shown as the gray region.

SUMMARY AND DISCUSSION
We have conducted a pilot study to search for IMBHs in dense stellar systems, one that uses TDEs as a probe of the presence of 10 3 M ⊙ ≲ M • ≲ 10 5 M ⊙ black holes embedded in stellar cusps at the center of massive GCs.Following previous work, we have modelled the long super-Eddington accretion phase in the slim advective disk regime together with the accompanying adiabatically expanding radiation-driven outflow.The ensuing optical/UV flare easily outshines the brightness of the (old) stellar population of globulars for a period of months to years, making TDEs triggered by IMBHs in principle detectable in a large sample of GCs.The disruption rate of main-sequence stars as a function of black hole mass and GC properties was estimated with a simple model of loss-cone dynamics and the help of a numerical 1D Fokker-Planck approach.
Large surveys with precise multi-wavelength photometry are becoming increasingly available for public use.As an illustrative example, we have taken advantage of the NGVS, an optical-near IR imaging survey of the Virgo galaxy cluster, and of its robust sample of GCs observed in the u-and g-band.We have checked for the presence of large photometric variations in the u-band induced by potential TDEs flares and found no obvious candidates.Since our model predicts as many as 10×f occ detectable events in the NGVS sample, the lack of recognizable candidates in the data implies that the fraction of GCs harboring IMBHs must be f occ ∼ < 10%.This is not very constraining, as post-merger recoil kicks originated by anysotropic GW emission may make it hard for IMBHs to be retained in lower-dispersion parent clusters (e.g.Arca Sedda et al. 2023).
Naturally, better modeling of the properties of the cluster central stellar distribution and of TDE light curves and rates are all needed to gain deeper insights into the presence of IMBHs in GCs.Large samples of extra-galactic star clusters, like those that will be assembled by the Vera C. Rubin Observatory (Usher et al. 2023), should enable significant progress in the search for TDEs in dense stellar systems.While there are considerable uncertainties in the modeling of TDE events in GCs, we hope that our pilot investigation will provide a blueprint for future searches of optical flares triggered by IMBHs.

Figure 1 .
Figure1.Bolometric light curves for the disk (purple lines) and wind (orange lines) emission resulting from the disruption of a solar-type star by a 3,000 M⊙ IMBH.The time t = 0 corresponds to the pericenter passage of the disrupted star, and the fallback time is t fb = 2.24 days.The solid and dashed lines show, respectively, the fout = 0.1 and fout = 0.9 cases for the fraction of infalling gas ejected in the wind (see text for details).The fout = 0.1 outflow becomes optically thin at t = 469 days.

Figure 4 .
Figure 4. of the stellar density profile of four IMBH+GC systems with initial Plummer core densities ρc = 10 4.5 M⊙ pc −3 (top-left panel), 10 5.0 M⊙ pc −3 (top-right panel), 10 5.5 M⊙ pc −3 (bottom-left panel), and 10 6.0 M⊙ pc −3 (bottomright panel).In all cases, a central IMBH of initial mass M•(t = 0) = 600 M⊙ is embedded in a pre-existing r −7/4 stellar cusp that extends all the way up to the influence radius (vertical lines in the panels).The starred points denote the locus where most TDEs are sourced.The collisional evolution of these spherical isotropic stellar systems under two-body relaxation and loss-cone effects was simulated by integrating the coupled Poisson and orbit-averaged 1D Fokker-Planck equations with the code PhaseFlow (Vasiliev 2017) for a single-mass stellar system composed by M⋆ = 1 M⊙ stars.

Figure 5 .
Figure 5.The growth of central IMBHs by stellar captures in the time-evolving stellar density profiles of Figure 4.In all configurations, the initial black hole mass was fixed to M•(t = 0) = 600 M⊙.Plummer core initial densities range from ρc(t = 0) = 10 4.5 M⊙ pc −3 (solid line) at the bottom to 10 6.0 M⊙ pc −3 (dotted line) at the top.

Figure 6 .
Figure 6.TDE rates in the time-evolving stellar density profiles of Figure 4.A single-mass stellar system composed by M⋆ = 1 M⊙ stars was assumed for simplicity in the numerical integration of the coupled Poisson and orbitaveraged 1D Fokker-Planck equations with PhaseFlow.In all cases, a central IMBH of initial mass M•(t = 0) = 600 M⊙ was embedded in a pre-existing r −7/4 stellar cusp (see text for details).

Figure 7 .
Figure 7. Stellar mass (left panel) and mean monochromatic luminosity (right panel) distributions of NGVS GCs.Stellar masses are computed using the z-band mass-to-light ratios predicted by the metallicity-dependent population-synthesis model PEGASE(Jordán et al. 2007).The vertical red dashed line marks the mean value of the distribution, ⟨MGC⟩ = 1.1 × 10 6 M⊙.In the u-band, most GCs are fainter than 10 39 erg s −1 .

Figure 8 .
Figure 8. of ∆tTDE timescales (see text for details) for our sample of GCs at varying fout, with the weaker winds+disk systems outshinining the steady stellar component for longer timescales because of late time accretion.

Figure 9 .
Figure 9. Three examples of TDE detections from our Monte Carlo simulations.Each light curve is compute in the u-band (solid orange line) for fout = 0.1 (top panel) and fout = 0.9 (bottom panel) with IMBH mass M• = 0.01 MGC.The solid blue line shows the GC mean monochromatic luminosity while the gray shaded region marks the total time span between the first and last epoch of cluster observation, ∆tS (see text for details).From top left to bottom right, the masses of the IMBHs responsible for the TDE are 1.39 × 10 4 , 1.26 × 10 5 , 9.25 × 10 4 , 3.02 × 10 4 , 1.26 × 10 5 , and 3.03 × 10 4 M⊙.