# The growth of intermediate mass black holes through tidal captures and tidal
disruption events

Francesco Paolo Rizzuto1,2, Thorsten Naab2, Antti Rantala2, Peter H.
Johansson1, Jeremiah P. Ostriker3,4, Nicholas C. Stone5, Shihong Liao1,
Dimitrios Irodotou1
1Department of Physics, University of Helsinki, Gustaf Hällströmin katu 2,
FI-00014, Helsinki, Finland
2Max-Planck Institute for Astrophysics, Karl-Schwarzschild-Str. 1,D-85741 ,
Garching, Germany
3Department of Astronomy, Columbia University, New York, NY 10027, USA
4Department of Astrophysical Sciences, Princeton University, Princeton, NJ
08544, USA
5Racah Institute of Physics, The Hebrew University, Jerusalem, 91904, Israel
Contact e-mail<EMAIL_ADDRESS>
###### Abstract

We present $N\mathrm{-body}$ simulations, including post-Newtonian dynamics,
of dense clusters of low-mass stars harbouring central black holes (BHs) with
initial masses of 50, 300, and 2000 $\mathrm{M_{\odot}}$. The models are
evolved with the $N\mathrm{-body}$ code bifrost to investigate the possible
formation and growth of massive BHs by the tidal capture of stars and tidal
disruption events (TDEs). We model star-BH tidal interactions using a
velocity-dependent drag force, which causes orbital energy and angular
momentum loss near the BH. About $\sim 20-30$ per cent of the stars within the
spheres of influence of the black holes form Bahcall-Wolf cusps and prevent
the systems from core collapse. Within the first 40 Myr of evolution, the
systems experience 500 – 1300 TDEs, depending on the initial cluster
structure. Most ($>95$ per cent) of the TDEs originate from stars in the
Bahcall-Wolf cusp. We derive an analytical formula for the TDE rate as a
function of the central BH mass, density and velocity dispersion of the
clusters ($\dot{N}_{\mathrm{TDE}}\propto M\mathrm{{}_{BH}}\rho\sigma^{-3}$).
We find that TDEs can lead a 300 $\mathrm{M_{\odot}}$ BH to reach $\sim
7000\mathrm{M_{\odot}}$ within a Gyr. This indicates that TDEs can drive the
formation and growth of massive BHs in sufficiently dense environments, which
might be present in the central regions of nuclear star clusters.

###### keywords:

galaxies: nuclei – quasars: supermassive black holes – black hole mergers –
galaxies: kinematics and dynamics – methods: numerical

††pubyear: 2020††pagerange: The growth of intermediate mass black holes
through tidal captures and tidal disruption events–References Figure 1: Left
panel: schematic illustration of a TDE. A star approaches a black hole moving
in a wide orbit with a pericentre smaller than the tidal radius
$R_{\mathrm{t}}$ (identified with the blue sphere in the plot). As soon as the
star - BH separation is smaller than $R_{\mathrm{t}}$ the star experiences a
tidal pull strong enough to tear it apart. Right panel: schematic illustration
of a TCE. In this case, the star follows an unbound orbit with pericentre
$\gtrsim R_{\mathrm{t}}$. Therefore the star is not destroyed. However, the BH
is sufficiently close to cause internal oscillations in the star, converting a
fraction of the initial orbital energy of the star into internal energy. The
orbit becomes bound: the star is "captured". If the same process repeats at
every pericentre passage, the stellar orbit will shrink over time, and can
eventually lead to most of the star’s mass accreting onto the BH.

##  1 Introduction

It is well established that most massive galaxies host in their centre black
holes with masses above $10^{6}\mathrm{M_{\odot}}$, also known as supermassive
black holes (SMBHs) (Magorrian et al., 1998; Kormendy & Ho, 2013).
Historically, the first strong indications of their presence were active
galactic nuclei (AGNs), for which the high luminosity and electromagnetic
spectrum were explainable only by the accretion of matter by supermassive
compact objects (Lynden-Bell, 1969). In recent years, observers utilising the
Event Horizon Telescope have even been able to produce the first direct image
of the shadow of the $\approx 6\times 10^{9}\mathrm{M_{\odot}}$ SMBH at the
centre of the Messier 87 galaxy (Event Horizon Telescope Collaboration et al.,
2019) as well as Sagittarius A*, the $\approx 4\times
10^{6}\mathrm{M_{\odot}}$ SMBH in the centre of the Milky Way (Event Horizon
Telescope Collaboration et al., 2022). Similarly, the existence of stellar
mass BHs $M\mathrm{{}_{BH}}<100\mathrm{M_{\odot}}$, the end product of massive
star evolution, has been confirmed by several X-ray and optical observations
(Casares et al., 2014; El-Badry et al., 2022; Remillard & McClintock, 2006)
and more recently even by gravitational wave (GW) detections (Abbott et al.,
2016, 2017; The LIGO Scientific Collaboration et al., 2021).

On the other hand, black holes with masses between
$10^{2}\mathrm{M_{\odot}}\lesssim M\mathrm{{}_{BH}}\lesssim
10^{6}\mathrm{M_{\odot}}$, referred to as intermediate-mass BHs (IMBHs),
remain more elusive. Over the last decades, several observations seem to
indicate the presence of IMBHs in galactic nuclei. More than thirty years ago,
Kunth et al. (1987) pointed out an AGN in the dwarf galaxy POX 52, which might
be powered by a $\sim 10^{5}\mathrm{M_{\odot}}$ BH (Barth et al., 2004).
Similarly, a $10^{4}-10^{5}\mathrm{M_{\odot}}$ BH could be responsible for the
Seyfert activity detected in the nuclear star cluster of the dwarf spiral
galaxy NGC 4395 (Filippenko & Sargent, 1989; Shih et al., 2003; Peterson et
al., 2005; den Brok et al., 2015) ) l.77. Many objects with similar masses
have been found in several subsequent surveys (see Greene et al., 2020;
Reines, 2022, and references therein). Recently, Gültekin et al. (2022) showed
that at least five sources located in low-mass galaxies (with stellar masses
$<3\times 10^{9}\mathrm{M_{\odot}}$), are consistent with active BHs with
masses between
$10^{4.9}\mathrm{M_{\odot}}<M\mathrm{{}_{BH}}<10^{6.1}\mathrm{M_{\odot}}$. The
presence of such BHs might not be limited to the centres of galactic nuclei; a
few of the most massive globular clusters might also be hosting objects of
similar masses (Farrell et al., 2014; Pechetti et al., 2022). Recently, a
tidal disruption event (TDE) was discovered (Lin et al., 2018) in a dense
stellar object of mass $\sim 10^{7}\mathrm{M_{\odot}}$, likely either a large
globular cluster or a tidally stripped satellite nucleus. X-ray continuum
fitting indicates that the TDE was powered by an IMBH of mass $\sim
10^{4}\mathrm{M_{\odot}}$ (Wen et al., 2021).

The physical processes leading to the formation of these BHs are still debated
in the literature (see Volonteri et al., 2021, for a comprehensive review).
One of the most promising scenarios for growing such massive BHs is the
collisional scenario. Low-mass BHs (a few tens of solar masses) located in
dense stellar environments could grow in mass through stellar and compact
object mergers (see, for example Stone et al., 2017). Gravitational wave
observations could have already revealed the tip of the iceberg of this
process. The LIGO-Virgo-KAGRA interferometers have thus far detected about
ninety BH coalescences, a few of which generated BHs with masses above
$\gtrsim 100\mathrm{M_{\odot}}$ (The LIGO Scientific Collaboration et al.,
2021). Being the most massive members of many star clusters, stellar-mass BHs
are expected to sink in the very centre of the host stellar system. There, if
the density is high enough, a chain of hierarchical mergers could produce BHs
with masses around $10^{2}\mathrm{M_{\odot}}-10^{4}\mathrm{M_{\odot}}$
(Atallah et al., 2022; Arca-Sedda et al., 2021; Fragione et al., 2022; Mapelli
et al., 2021; Rizzuto et al., 2021) or even $10^{5}\mathrm{M_{\odot}}$
according to Antonini et al. (2019). The main obstacle to this hierarchical
formation path is recoil caused by the anisotropic emission of gravitational
radiation. The final product of compact object mergers is expected to receive
GW-induced velocity kicks that range from a few tens up to $5000$ km/s
depending on the mass ratio and the spins of the merging bodies (Campanelli et
al., 2007; Lousto et al., 2010; Lousto & Healy, 2019). Systems with low escape
velocities, such as globular clusters, are unlikely to retain the final
product of repeated BH mergers (Arca-Sedda et al., 2021; Gerosa & Berti,
2019). Only dense stellar environments with an escape velocity in excess of
$v_{\mathrm{esc}}\gtrsim 100$ km/s have a non-negligible chance to grow IMBHs
through hierarchical collisions (Gerosa & Berti, 2019; Mapelli et al., 2021).
Nuclear star clusters are the only class of stellar systems that fulfill this
condition (Neumayer et al., 2020).

In star cluster environments, stellar BHs coexist with stars, which provide an
alternative merger channel for BH mass growth, especially in young compact
clusters which are populated by many massive stars. Direct $N\mathrm{-body}$
simulations presented in Portegies Zwart et al. (2004) showed that young dense
low-metallicity star clusters, if compact enough, can trigger repeated stellar
mergers and generate very massive stars (VMSs) with masses up to $\sim
1000\mathrm{M_{\odot}}$, which collapse directly into massive BH without
losing mass in the process. Recent direct $N\mathrm{-body}$ simulations evolve
similar systems using up-to-date stellar evolution recipes and confirm the
collisional formation scenario of VMSs (Di Carlo et al., 2021; Rizzuto et al.,
2022). However, they also indicate that such a process can lead at most to BHs
with masses of a few $100\mathrm{M_{\odot}}$ in agreement with the latest
Monte-Carlo simulations (Kremer et al., 2020; González Prieto et al., 2022).
The evolution of VMSs is highly uncertain. Glebbeek et al. (2009) argues that
the final product of massive stellar mergers is affected by enhanced mass loss
that prevents VMSs from generating BH more massive than
$10^{2}\mathrm{M_{\odot}}$.

Young massive star clusters may have another channel for massive BHs to grow:
mergers between BHs and massive stars. Such events have been reported in both
direct $N\mathrm{-body}$ simulations (Mapelli, 2016; Rizzuto et al., 2021;
Rizzuto et al., 2022) and Monte-Carlo simulations (Giersz et al., 2015). In
Rizzuto et al. (2021) we show that BHs up to $\sim 3\times
10^{2}\mathrm{M_{\odot}}$ can form even when the VMS direct collapse channel
is suppressed111With the stellar evolution recipes adopted in Rizzuto et al.
(2021) VMSs above $>10^{2}\mathrm{M_{\odot}}$ form BHs lighter than
$30\mathrm{M_{\odot}}$. as long as stellar BHs can accrete a significant
fraction of the stellar material when colliding with a massive star.

Young massive star clusters can generate massive BHs only during the early
stage of their evolution when they are still dense, and massive stars are
still present in the system. As soon as stars undergo supernova explosions,
the clusters lose a large fraction of their central mass and subsequently
experience a rapid expansion. The most massive stars terminate their life,
collapsing into stellar black holes. As these BHs are the heaviest objects
left in the system, they segregate in the cluster core and form a compact
subsystem. It is well established in the literature that the BH subsystem
prevents the cluster core-collapse, drives stars away from the centre, and can
even remove them from the cluster(Merritt et al., 2004; Breen & Heggie, 2013).
At the same time, also BHs can be dynamically ejected during repeated few-body
close encounters. Because of these interactions, over time, star clusters can
lose most of their BHs. Which of the two components, BHs or stars, will
dissolve first depends on whether the cluster is tidally filling222A cluster
is tidally filling when it contains stars all the way to its tidal radius,
otherwise it is tidally underfilling. Typically compact and dense systems are
tidally underfilling while more dilute clusters with a moderate initial
central density are tidally filling. at the beginning or during its evolution
(Giersz et al., 2019). In tidally filling clusters, stars escape faster than
BHs. These systems evolve toward a state in which the BHs constitute the
majority of their mass. The low-density globular cluster Palomar 5 has most
likely entered this phase since direct $N\mathrm{-body}$ simulations suggest
that about $20$ per cent of its mass is in the form of BHs (Gieles et al.,
2021).

On the other hand, in tidally underfilling clusters, the BH component
evaporates first, leaving behind only a few black holes immersed in a sea of
low-mass stars. An example of such a system is the $\sim
10^{5}\mathrm{M_{\odot}}$ globular cluster NGC 6397, which hosts no more than
a dozen of BHs in its core (Weatherford et al., 2020). The presence of an IMBH
would further catalyse the removal of the BH subcluster, as simulations show
that the massive compact object is likely to rapidly expel most BHs and BH
binaries during few-body interactions (Giersz et al., 2015). In this case the
star cluster would quickly evolve into a system consisting of low-mass stars
enclosing a central massive BH. During this stage the BH can slowly consume
the surrounding stars through tidal disruption events or tidal capture events
(TCEs). The former occurs when the BH is close enough to a star to exert a
tidal pull that rips it apart (left panel of Fig. 1). The latter case occurs
when an unbound star is sufficiently near to a BH to experience a strong
enough tidal perturbation to become bound. The tidal force deforms the nearby
star, so part of its initial orbital energy is transferred to its internal
energy. This mechanism can convert unbound stars onto a bound orbit. When this
happens, the star is said to be "captured" by the BH (see right panel of Fig.
1 as an illustration of a TCE). After being captured, the star can feed the
compact object in various ways, such as, partial disruption events at each
pericentre passage (Kremer et al., 2022), mass transfer events, or even
dynamically induced mergers. TCEs were highlighted for the first time in
Fabian et al. (1975); Press & Teukolsky (1977) provided the first mathematical
description of these events, later on, refined by Lee & Ostriker (1986a).

Since the tidal capture mechanism has a larger cross-section than the tidal
disruption process, it is likely to impact BH mass growth in dense stellar
environments significantly. The analytical work led by Stone et al. (2017)
provided, for the first time, a comprehensive study on the critical role
played by TCEs and TDEs in forming massive BHs in galactic nuclei. Their work
showed that a low mass BH located in the centre of a dense stellar environment
can reach, within a Hubble time, up to $10^{6}\mathrm{M_{\odot}}$ through
tidal disruptions and tidal capture runaway collisions. According to their
calculation, the rapid growth can be triggered only in clusters with a core
density of $n_{\mathrm{c}}>10^{7}/$pc3 and a core velocity dispersion of
$\sigma_{\mathrm{c}}>40$ km/s.

In this work, we explore the tidal collisional runaway scenario using direct
$N\mathrm{-body}$ simulations. Therefore, we evolve star clusters initialised
according to the density and velocity dispersion criteria given in Stone et
al. (2017), and study the dynamic mechanisms responsible for the growth of
massive BH. In Section 2, we describe in detail the numerical integrator and
the tidal interactions prescriptions adopted for our investigation. In Section
3, we discuss the intial conditions of our clusters. In Section 4 we present
the results of our simulations and in Section 5 we offer a summary of our main
results.

##  2 Methods

In order to investigate the growth of massive BHs in compact stellar
environments through TDEs and TCEs, we run five direct $N\mathrm{-body}$
simulations of dense star clusters consisting of a central BH surrounded by
low-mass stars. We evolve the systems for at least $40$ Myr using the direct
$N\mathrm{-body}$ integrator bifrost (Rantala et al., 2022), which we describe
in detail in the following Subsection. We include in bifrost prescriptions for
modelling the effect of tidal interactions during close encounters using a
drag force as described in Subsection 2.4.

###  2.1 The N-body code

In this study we used the novel GPU-accelerated direct-summation
$N\mathrm{-body}$ simulation code bifrost (Rantala et al., 2022) which is
based on the earlier FROST code (Rantala et al., 2021). The code uses a
hierarchical implementation (HHS-FSI, Rantala et al. 2021) of the fourth-order
forward symplectic integrator (FSI, e.g. Chin 1997; Chin & Chen 2005; Dehnen &
Hernandez 2017). In addition to the Newtonian accelerations, FSI uses
additional so-called gradient accelerations to cancel second-order error terms
with strictly positive integrator sub-steps. This is different to the widely
used Yoshida-type symplectic integrators (Yoshida, 1990) which always contain
negative sub-steps in integrator orders higher than two. For the Kepler
problem it has been shown that fourth-order symplectic integrators with
strictly positive sub-steps outperform the integrators that include negative
sub-steps (Chin, 2007). Compared to the common block time-step scheme widely
used in $N\mathrm{-body}$ simulations (Aarseth, 2003), HHS-FSI is manifestly
momentum-conserving due to pair-wise acceleration calculations and
synchronized kick operations. In addition, rapidly evolving parts of the
simulated systems decouple from slowly evolving parts in the hierarchical
integration, making the approach especially efficient for $N\mathrm{-body}$
systems with a large dynamical range (Pelupessy et al., 2012).
Besides the forward integrator, bifrost includes regularised (Rantala et al.,
2020) and secular integration techniques for subsystems, i.e. binaries, triple
systems and small clusters around massive black holes. The equations of motion
of the particles in the subsystems are post-Newtonian (PN) up to order PN3.5
using the formulas of Thorne & Hartle (1985); Blanchet et al. (2006). For
binary systems the post-Newtonian equations of motion enable relativistic
precession effects due to conservative even PN terms, as well as orbit
circularisation and inspiral due to radiation-reaction forces due to
gravitational wave emission. bifrost can evolve massive star clusters with
high numerical accuracy containing up to a few million stars with an arbitrary
fraction of primordial binaries (Rantala et al., 2022), a feature which only
few current codes share (Wang et al., 2020).
The basic bifrost version has four prescriptions for mergers. First, two
particles are merged if their gravitational-wave inspiral timescale becomes
shorter than their current time-steps. We also merge particles if their mutual
separation is shorter than the relativistic innermost stable circular orbit.
bifrost also includes a simple prescription for tidal disruption events which
is significantly extended for this study as described in the next sections.
Finally, two stars are merged if they directly collide, i.e. their radii
overlap. Table 1 reports the main bifrost code parameters used in our runs.

bifrost user-given parameter | symbol | value
---|---|---
forward integrator time-step factor |  $\eta_{\mathrm{ff}}$, $\eta_{\mathrm{fb}}$, $\eta_{\mathrm{\nabla}}$ | $0.2$
subsystem neighbour radius | $r_{\mathrm{rgb}}$ |  $0.008$ pc
regularization GBS tolerance | $\eta_{\mathrm{GBS}}$ | $10^{-10}$
regularization GBS end-time tolerance | $\eta_{\mathrm{endtime}}$ | $10^{-4}$
regularization highest PN order |  | PN3.5
secular integration threshold | $N_{\mathrm{orb,sec}}$ | $10$
secular highest PN order |  | PN2.5
Table 1: The main user-given code bifrost parameters used in the simulation
runs in this study. The parameter definitions correspond to the ones in
Rantala et al. (2022). Figure 2: Left panel: Evolution of an initially unbound
orbit of a tidally affected $1.5\mathrm{M_{\odot}}$ star around a
$50\mathrm{M_{\odot}}$ black hole. Due to tidal energy losses (implemented in
the simulation as a drag force), the star loses orbital energy and its orbit
circularises. Right panel: Time evolution (in log-scale) of the orbital
eccentricity (black) and the pericentre distance (solid red line) of the star.
In the final stages the orbit circularises and the pericentre distance
shrinks. Once the pericentre becomes smaller than the tidal radius (red dot-
dashed line) the star is disrupted in a TDE.

###  2.2 Prescription for tidal disruption and mass accretion

In bifrost, a BH destroys a star when $r$, the distance between the two
interacting objects is smaller than the tidal radius $R_{\mathrm{t}}$. To
compute $R_{\mathrm{t}}$ we adopt the criterion given in Kochanek (1992):

$R_{\mathrm{t}}=1.3R_{*}\left(\frac{M\mathrm{{}_{BH}}+m_{*}}{2m_{*}}\right)^{1/3},$
(1)

where $R_{*}$ and $m_{*}$ are the radius and the mass of the destroyed star,
respectively, while $M\mathrm{{}_{BH}}$ is the BH mass. To ensure the codes
does not miss any TDEs, we decrease the time steps of each interaction with
pericentre $R_{\mathrm{p}}<3R_{\mathrm{t}}$, for the code to check the
condition $r<R_{\mathrm{t}}$ exactly at pericenter passage. With this
restriction, the regularised integrator is less efficient, but fortunately,
these close interactions are infrequent. Therefore they have a negligible
impact on the overall performance of the code.

After a TDE, a fraction of the stellar material is expected to be ejected at
high velocity, while the remaining mass remains bound to the BH and eventually
forms an accretion disk around it. Recent hydrodynamic simulations of tidal
disruption events between stellar BHs and main-sequence stars show that the
fraction of stellar material bound to the BH could be close to unity (Kremer
et al., 2022). However, these simulations do not model the accretion phase in
which a significant fraction of the initially bound gas may get unbound,
especially if the accretion occurs at a super-Eddington rate (Ayal et al.,
2000; Bonnerot & Lu, 2020; Dai et al., 2018; Metzger & Stone, 2016; Steinberg
& Stone, 2022; Toyouchi et al., 2021). In general, many aspects of the
accretion phase are poorly constrained. For example, it is not very well
understood how much stellar mass ends up into the BH and how much is lost, and
the rate at which the accretion should occur is still debated in the
literature. In this work we follow the simple estimate in Rees (1988) and
assume that in a tidal disruption event, $50$ per cent of the stellar mass is
accreted by the BH instantly333During these instantaneous accretion events we
assume that linear momentum is conserved., and the other $50$ per cent is
instantly removed from the cluster. We discuss how this simplification affects
our results in Subsection 4.6.

###  2.3 Tidal interaction energy loss

When a star of mass $m_{*}$ moves in an orbit with a pericentre slightly
larger than the tidal radius $R_{\mathrm{p}}\gtrsim R_{\mathrm{t}}$, tidal
forces are not strong enough to rip the star apart, but they can still deform
the star triggering internal oscillations. Consequently, a fraction of the
orbital energy of the star is lost in the internal energy of the star. In
other words, tidal interactions force stars to deviate from their original
orbits: bound eccentric orbits become more bound and circular, and unbound
parabolic or hyperbolic orbits become bound.

To estimate the energy lost during a parabolic encounter, Press & Teukolsky
(1977) approximated the oscillations amplitude of the perturbed star through
an expansion in spherical harmonics. With their procedure, the fraction of
orbital energy deposited into the star is:

$\Delta
E=\int_{-\infty}^{+\infty}\frac{\mathrm{d}E}{\mathrm{d}t}\mathrm{d}t\approx\frac{Gm_{*}^{2}}{R_{*}}\left(\frac{M\mathrm{{}_{BH}}}{m_{*}}\right)^{2}\sum_{l=2}^{\infty}\left(\frac{R_{*}}{R_{\mathrm{p}}}\right)^{2l+2}T_{l}(\eta).$
(2)

Here $T_{l}$ are dimensionless coefficients that depend on the internal
structure of the deformed object. They are functions of the parameter $\eta$
defined as:

$\eta\coloneq\left(\frac{m_{*}}{m_{*}+M\mathrm{{}_{BH}}}\right)^{1/2}\left(\frac{R_{\mathrm{p}}}{R_{*}}\right)^{3/2}.$
(3)

As follows from its definition, $\eta$ indicates the duration of the
pericentre passage with respect to the hydro-dynamical timescale of the star
$m_{*}$. In practical applications, only $l=2$ and $l=3$ are taken into
account. $T_{0}=0$ and higher terms give a negligible contribution to the
final value of $\Delta E$. The values of $T_{2}$ and $T_{3}$ have been
calculated explicitly for objects with polytropic indices $n=1.5,2$ and $n=3$
(see Lee & Ostriker, 1986b; Ray et al., 1987). Based on these values,
Portegies Zwart & Meinen (1993) provide fitting functions to rapidly estimate
$T_{2}$ and $T_{3}$ during close encounters in $N\mathrm{-body}$ simulations.

The formula shown in Eq. 2, has been derived exclusively for parabolic
encounters ($e=1$). To extend this formulation for eccentricities larger or
smaller than $1$, we utilise the prescription given in appendix A of Mardling
& Aarseth (2001), which provides a generalisation of Eq. 2 by introducing a
new expression for $\eta$ (indicated with $\zeta$) that depends explicitly on
the eccentricity of the orbit:

$\zeta\coloneq\eta\left(\frac{2}{1+e}\right)^{\alpha(\eta)/2}$ (4)

where $\alpha(\eta)=1+\frac{1}{2}\left|\frac{\eta-2}{2}\right|$, therefore for
$e=1$ we restore the original formulation. Numerical tests show that using
$T_{l}(\zeta)$ instead of $T_{l}(\eta)$, the tidal evolution of hyperbolic
orbits ($e\gtrsim 1$) is consistent with a more accurate model for tidal
interactions presented in Mardling (1995). In addition, $T_{l}(\zeta)$ should
give more reliable results for bound orbits with $e\gtrsim 0.5$ (see Mardling
& Aarseth, 2001). To estimate the tidal energy loss during a BH - star close
encounter in our simulations we first compute $\zeta$ using Eq. 4 and then
evaluate $T_{2}(\zeta)$ and $T_{3}(\zeta)$ using the interpolation formulas
given in Portegies Zwart & Meinen (1993). Finally, we use Eq. 2 to estimate
the energy loss in tidal interactions. For orbits with $e<0.5$, we assume the
tides switch from the dynamical to the equilibrium regime. The amount of tidal
energy dissipation, therefore, decreases with decreasing eccentricity. To take
into account this effect, we multiply the energy dissipation $\Delta E$ by the
eccentricity $e$ ($\Delta E=e\times\Delta E$) similarly to what was done by
Baumgardt et al. (2006).

###  2.4 Modelling tidal interactions with a drag force

In the previous subsection, we described the procedures adopted to estimate
$\Delta E$, the amount of energy dissipated during tidal interactions. There
remains the problem of how to use $\Delta E$ to change the orbit of the two
closely interacting objects. The standard approach assumes that all the energy
$\Delta E$ is instantly emitted at each pericentre passage. Then under the
assumption that the angular momentum is conserved, one can derive the new
eccentricity and semi-major axes from the present one (Mardling & Aarseth,
2001) and updates the orbit accordingly. The entire procedure must be
implemented outside the regularised integrator. In practice, to resolve the
trajectory, the integrator must evolve the orbit to the pericentre, apply the
orbital change due to tidal energy loss, then continue the evolution to the
next pericentre passage. This procedure is relatively straightforward to
implement in isolated two-body encounters, but it becomes very cumbersome when
three or more objects are involved in the interaction. Integrating the effects
of tidal interactions directly in the regularisation would be significantly
more convenient. For this reason, we decided to model tidal interactions using
a drag force following Samsing et al. (2018):

$\vec{F}_{\text{t}}(r,v)=-C\frac{\vec{v}}{r^{4}}$ (5)

where $r$ and $v$ are the relative separation and relative velocity of the two
tidally interacting bodies, respectively, while $C$ is a normalisation factor:

$\int_{\mathrm{orb}}F_{t}(r,v)dr=\Delta E.$ (6)

Here $\Delta E$ is the same quantity computed in the instant emission
treatment (Eq. 2). With this definition, $C$ ensures the two method dissipate
the same amount of energy per orbit.

The force expression of Eq. 5 closely resembles the formula of the dissipative
PN 2.5 radiation-reaction force. For this reason, Eq. 5 can be included with
little effort in the bifrost regularised integrator alongside the PN terms.
The strong dependence of $F_{\text{t}}$ on $r$, the distance of the two
interacting bodies, ensures that tidal interactions are effectively activated
during very close encounters and are negligible at large distances.

Figure 3: Left panel: Evolution of an initially unbound orbit of a tidally
affected $0.5\mathrm{M_{\odot}}$ star around a $50\mathrm{M_{\odot}}$ black
hole. The evolution is similar to the $1.5\mathrm{M_{\odot}}$ star but on
shorter timescales. Right panel: Time evolution (in linear-scale) of the
orbital eccentricity (black) and the pericentre distance (solid red line) of
the star. In the final stages the orbit circularises and the pericentre
distance shrinks. Once the pericentre becomes smaller than the tidal radius
(red dot-dashed line) the star is disrupted in a TDE.

###  2.5 Tidal interactions drag force and orbital evolution

The instant emission prescription for tidal interactions assumes conservation
of angular momentum; therefore the pericentre increases until the orbit
circularizes (see Mardling & Aarseth, 2001, for more details). On the
contrary, the tidal interaction drag force does not conserve angular momentum.
Therefore, with this method the pericentre of tidally captured stars decreases
over time until $R_{\mathrm{p}}<R_{\mathrm{t}}$. In other words, in the
absence of external perturbers, the tidal interaction drag force drives every
TCE into a TDE . The time required to bring a captured star to tidal
disruption depends on the star’s initial pericentre, initial eccentricity and
internal structure. In this work, we model low-mass stars
($m_{*}<0.7\mathrm{M_{\odot}}$) with the polytropic index $n=1.5$ and massive
stars ($m_{*}>0.7\mathrm{M_{\odot}}$) with polytropic index $n=3.0$.
Consequently, massive stars are more compact and less affected by tidal
perturbations than low-mass stars. Fig. 2 shows the orbital evolution (left
panel) of a $1.5\mathrm{M_{\odot}}$ star captured by a $50\mathrm{M_{\odot}}$
BH. The star initially moves in a quasi-parabolic orbit with $e=1.0001$ and
$R_{\mathrm{p}}=1.5R_{\mathrm{t}}$. In the beginning, the eccentricity
decreases while the pericentre stays constant. In the last stage of the
evolution, $R_{\mathrm{p}}$ drops rapidly and reaches the tidal radius. It
takes about $100$ yr for the star to be destroyed (see right panel of Fig. 2).
The $0.5\mathrm{M_{\odot}}$ star shown in Fig. 3 experiences a similar
evolution. However, the effect of the drag force is much stronger when acting
on the $0.5\mathrm{M_{\odot}}$ star. Consequently, the low-mass star spiral
towards the BH in just 0.6 years.

As shown in Figs. 2 and 3, the stars, after being captured, undergo a intial
phase of circularisation (the eccentricity drops) while the pericentre remains
unchanged. Only in the last phase of evolution, also the pericentre begins to
decrease. We will use this feature of the drag force to identify all the TCEs
that occurred in our simulations (see Subsection 4.4).

To summarise, in our simulations, we model the effect of tidal interactions
using the drag force described by Eq. 5. With such a prescription, every TCE
will produce a TDE in less than $\lesssim 10^{3}$ yr. It must be said that the
outcome following a TC event is highly debatable. In this paper, we assume
that every tidal capture leads rapidly to tidal disruption. However, this may
not be true. A star, after being captured, if not affected by external
perturbations, follows a bound but very eccentric orbit. If internal
mechanisms (such as viscosity) of the star are adequately efficient to rapidly
dissipate the oscillatory motion of the star, at the second pericentre
passage, the star will experience again a strong tidal interaction that leads
to orbital energy loss and consequently to the shrinking of the semi-major
axis. The process repeats until the star gets destroyed and swallowed by the
BH, or the two bodies circularise. If the spin-orbit coupling is inefficient,
orbital angular momentum is conserved, and the star circularises around the
BH. If the spin-orbit interaction is efficient, a fraction of orbital angular
momentum can be transferred efficiently into the stellar spin.

Another thing to consider is the possible expansion of the stars after the
pericentre passage due to the gain in internal energy. In this scenario, the
radius of the star $R_{*}$ will increase, resulting in an increased tidal
radius: the star will be tidally disrupted after a few pericentre passages.
Another possible scenario is described in Mardling & Aarseth (2001), where
they argue that if the oscillations induced by tidal interaction are not
damped efficiently, TCEs might lead to chaotic evolution. In fact, the
pericentre passages that follow the capture can further excite but also damp
oscillation modes in the deformed star. In other words, the orbit and the star
can randomly exchange energy packages in both directions. The resulting
trajectories of the object are unpredictable and resemble a random walk.
Recent hydrodynamical simulations indicate that captured stars that experience
partial tidal disruptions to be completely destroyed after a few pericentre
passages (Kremer et al., 2022). Unfortunately, the simulations sample is too
small to constrain the fate of tidally captured stars extensively.

##  3 Initial Conditions

Name | $M\mathrm{{}_{BH}}$ | $R_{\mathrm{h}}$ | $r_{\mathrm{c}}$ | $n_{\mathrm{c}}$ | $\sigma_{\mathrm{c}}$
---|---|---|---|---|---
- | $\mathrm{M_{\odot}}$ | pc | pc |  $\\#$ stars / pc 3 | km/s
R04M300 | 300 | 0.4 | 5.4e-03 | 1.4e+09 | 42.3
R06M300 | 300 | 0.6 | 7.5e-03 | 1.7e+08 | 32.9
R08M300 | 300 | 0.8 | 1.1e-02 | 6.8e+07 | 28.5
R06M50 | 50 | 0.6 | 1.7e-02 | 5.5e+07 | 26.3
R06M2000 | 2000 | 0.6 | 3.3e-03 | 1.7e+09 | 90.4
Table 2: Initial clusters and BH properties: Name: name of the model; $M\mathrm{{}_{BH}}$: initial BH mass; $R_{\mathrm{h}}$: initial half mass radius; $r\mathrm{{}_{c}}\,$: initial core radius; $n\mathrm{{}_{c}}\,$: initial central particle density; $\sigma_{\mathrm{c}}$: initial velocity dispersion.  Name | t | $M\mathrm{{}_{BH}}$ | N${}_{\text{TDE}}$ | $R_{\mathrm{h}}$ | $r_{\mathrm{c}}$ | $\rho_{\mathrm{c}}$ | $\sigma_{\mathrm{c}}$
---|---|---|---|---|---|---|---
- | Myr | $\mathrm{M_{\odot}}$ | - | pc | pc |  $\mathrm{M_{\odot}}$ / pc 3 | km/s
R04M300 | 41.2 | 1152 | 1239 | 0.5 | 2.0e-02 | 8.8e+06 | 33.3
R06M300 | 148.0 | 1329 | 1487 | 0.9 | 6.6e-02 | 8.5e+05 | 23.6
R08M300 | 41.1 | 686 | 532 | 0.9 | 5.3e-02 | 1.6e+06 | 21.4
R06M50 | 41.9 | 556 | 665 | 0.7 | 2.9e-02 | 5.0e+06 | 25.1
R06M2000 | 41.7 | 2844 | 1342 | 0.8 | 1.6e-02 | 1.3e+07 | 51.3
Table 3: Cluster and BHs properties at the end of each simulation: Name: name
of the model; t: simulation run time. $M\mathrm{{}_{BH}}$: final BH mass;
N${}_{\text{TDE}}$: total number of TDE; $R_{\mathrm{h}}$: final half mass
radius; $r\mathrm{{}_{c}}\,$: final core radius; $\rho\mathrm{{}_{c}}\,$:
final central density; $\sigma$: final velocity dispersion.

We generated the initial conditions for five very dense star clusters with
$256000$ single stars and no primordial binaries using the code mcluster
(Küpper et al., 2011). The masses of the stars are sampled from Kroupa (2001)
initial mass function with a range from $0.08\mathrm{M_{\odot}}$ up to
$2.00\mathrm{M_{\odot}}$444Having chosen the most massive star to be
$2\mathrm{M_{\odot}}$, allows us to neglect the effect of stellar evolution,
and focus only on the dynamical evolution of the systems.. We located at the
centre of the cluster single BHs with initial masses equal to
$50\mathrm{M_{\odot}}$, $300\mathrm{M_{\odot}}$ and $2000\mathrm{M_{\odot}}$.
The five models are initialized with primordial mass segregation and at
starting time they follow a King (1966) density profile with W${}_{0}=9$;
three of these have an initial half mass radius of $R_{\mathrm{h}}=0.6$ pc,
the other two have a half mass radius respectively equal to $0.4$ pc and $0.8$
pc. From now on, we will refer to the five simulations using the labels
reported in Table 2555The core radius reported in the table is computed using
$r\mathrm{{}_{c}}=\sqrt{(\sum_{i}\rho_{i}^{2}r_{i}^{2})/(\sum_{i}\rho_{i}^{2})}$,
where $\rho_{i}$ is the local density around the $i$-particle and $r_{i}$ is
its distance from the centre..

Thus initialised, the systems have very high core densities (from $\sim
6\times 10^{7}$ particles per pc3 for R06M50 up to $\sim 2\times 10^{9}$
particles per pc3 for R06M2000) and high central velocities dispersion
(ranging from $30$ km/s up to $90$ km/s) and they might resemble the
conditions inside the cores of nuclear star clusters. The initial conditions
we chose for our cluster match the criteria indicated by Stone et al. (2017)
to trigger a tidal capture runaway collision. As discussed in the introduction
the work by Stone indicate that stellar environments with
$n_{\mathrm{c}}>10^{7}$ stars per pc3 and a central velocity dispersion
$\sigma_{\mathrm{c}}>40$km/s are expected to trigger tidal capture runaway
collisions.

##  4 Results

All five models maintain a significant tidal disruption rate
$\dot{N}_{\mathrm{TDE}}$ throughout the simulation. The final number of TDEs
depends strongly on the initial size ($R_{\mathrm{h}}$) of the cluster. For
instance, R04M300, the most compact model, registers a total of 1239 TDEs in
about 40 Myr. On the other hand, in the same evolution time, the least compact
model, R08M300 records less than half the TDEs (see Table 3).

The primary mechanism that leads stars being close enough to the BH to be
destroyed or captured can be broken down into two parts. First, dynamical
interactions in the vicinity of the BH build up a cloud of bound orbits, which
accompany the BH until the end of the simulation. Once these bound orbits are
formed, since they are located in the inner part of the cluster, they are
constantly perturbed by many gravitational scattering events. As a consequence
of these interactions the bound orbit become over time more bound and very
eccentric ($e\sim 1$), resulting in their tidal capture or disruption. In
simple words, the BH first forms rapidly a bound cloud in its vicinity and
later feeds on it slowly over time.

We provide a comprehensive description of the main properties of the bound
cloud in the different simulations in Subsections 4.1 and 4.2. The bound cloud
does not only influence the BH mass growth, but it also plays an important
role in the evolution of the cluster core as we discuss in Subsection 4.3. In
Subsection 4.4, we show the BHs grow primarily through direct TDEs666Stars are
destroyed at first pericentre passage.. However, also TCEs contribute to the
BHs mass growth, as they trigger about $10$ per cent of the total TDEs. We
study how the TDE rate changes over time in Subsection 4.5. There we provide
fitting formulas and an analytical approximation to estimate
$\dot{N}_{\mathrm{TDE}}$, which also helps to understand quantitatively the
dynamical processes that drive stars to disruption. Finally, in Subsection
4.6, we briefly discuss the phase of BH gas accretion after a star is tidally
destroyed. Moreover, we investigate how the BH mass growth changes varying the
fraction of stellar mass that ends up into the BH after a TDE.

Figure 4: Time evolution of the fraction of stars $(f_{\mathrm{b}})$ within
the BH influence radius ($R_{\mathrm{inf}}$) which are bound to the central
BH. Models R04M300, R06M300, and R08M300 start with about $\sim 25$ per cent
of bound stars, and they maintain this fraction until the end of the
simulations. R06M2000 starts with a fraction of about $\sim 50$ per cent that
drops to $\sim 30$ per cent. A bound cloud of stars in R06M50 forms within the
first 5 Myr reaching $\sim 20$ per cent after $\sim 10$ Myr.

###  4.1 Formation and evolution of the bound cloud

Our simulations reveal that, in the vicinity of the BHs, inside a region
enclosed within the influence radius $R_{\mathrm{inf}}$777The influence radius
$R_{\mathrm{inf}}$ defines a sphere centred on the BH that encompasses a
number of stars with total mass equal to $2M\mathrm{{}_{BH}}$., a significant
fraction of the stars is bound to the BH i.e. they have negative orbital
energy $E=\frac{m_{*}v^{2}}{2}-\frac{Gm_{*}M\mathrm{{}_{BH}}}{r}<0$, where
$m_{*}$ is the mass of the star. In Fig. 4, we plot, as a function of time,
$f_{\mathrm{b}}$ the fraction of stars within the influence radius bound to
the BH888By definition
$f_{\mathrm{b}}=\frac{N_{\mathrm{b}}}{N_{\mathrm{inf}}}$, where
$N_{\mathrm{b}}$ and $N_{\mathrm{inf}}$ are respectively the number of bound
stars and the total number of stars within the influence radius.. All the
models that started with a $M\mathrm{{}_{BH}}=300\mathrm{M_{\odot}}$ BH
display a stable bound fraction approximately equal to $\approx 0.2$. Fig. 4
shows that this value stays constant for at least $40$ Myr. Since we evolved
R06M300 for $150\,\mathrm{Myr}$, we know that the bound fraction is very
likely to remain constant for a significantly longer time (see Fig. 16 in
Appendix A). On the other hand, model R06M2000 initially has a fraction of
bound stars close to $50$ per cent, which drops rapidly in the first $\sim 10$
Myr and stabilises then at $30$ per cent, as indicated by the dark blue line
in Fig. 4.

In contrast to the other models, R06M50 starts its evolution without bound
orbits around the BH. However, the bound cloud forms rapidly in the first
$\sim 4$ Myr (orange line in Fig. 4). At the end of this build-up phase, the
fraction of bound stars levels out at about $20$ per cent, a similar value to
the one observed in R04M2000, R06M2000 and R08M2000.

Figure 5: Total stellar number density profiles (thin green line) for the
model R06M300 with a 300 $\mathrm{M_{\odot}}$ central BH at 5 Myr (left), 45
Myr (middle) and 100 Myr (right). The power-law fits indicate shallow profiles
with slopes of $\beta\sim 1.08-1.2$. The number density profiles of stars
bound to the central black hole (bold green lines) are significantly steeper
with power-law exponents of $\beta\sim 1.71-1.84$ (orange dot-dashed lines).
These slopes are consistent with the Bahcall & Wolf (1976) $\beta=7/4$ density
profile (see Fig. 6 for the time evolution of the power-law indices).  Figure
6: Time evolution of the power-law index $\beta$ for the stellar number
density (see Fig. 5) of all stars (orange) and stars bound to the central BH
(blue). The central BH mass is increasing from $50\mathrm{M_{\odot}}$ (R06M50,
left) to 300 $\mathrm{M_{\odot}}$ (R06M300, middle) and 2000
$\mathrm{M_{\odot}}$ (R06M2000, right). The total density profile slope
increases with increasing BH initial mass. It varies from $\beta\sim 1.0$ for
R06M50 (left panel) to $\beta\sim 1.4$ for R06M2000 (right panel). The slope
of the bound stars is independent of the black hole mass and mostly consistent
with a Bahcall-Wolf cusp ($\beta\sim 1.75$). This result also holds for all
other models and also on longer timescales (see Appendix A). The uncertainty
on the fitted slope, indicated in the plots with shadow regions, is estimated
using one hundred consecutive simulation snapshots.

###  4.2 Stellar distribution around the central black hole

When a single-mass stellar system with a massive central BH evolves for more
than a relaxation time, gravitational interactions drive the stellar
distribution in the vicinity of the BH to form a cusp with a particle density
profile: $n\propto r^{-1.75}$, where $r$ is the distance from the BH and $n$
is the stellar particle density. This cusp, known as a Bahcall-Wolf (BW) cusp,
was predicted for a stellar system with equal mass stars by Bahcall & Wolf
(1976) and thereafter it was generalised for systems with unequal mass stars
(Bahcall & Wolf, 1977; Alexander & Hopman, 2009). In particular, Alexander &
Hopman (2009) show that heavy stars, when they are relatively common, tend to
follow the original BW profile with slope $\beta=1.75$, while lighter stars
form a shallower cusp with $1.5\leq\beta\leq 1.75$. Direct $N\mathrm{-body}$
simulations of star clusters with a central SMBH confirm that stars very close
to the massive object follow a BW cusp and indicate that the cusp extends for
a tenth of the influence radius $\sim 0.1R_{\mathrm{inf}}$ (Preto et al.,
2004). Our runs show somewhat different results. The entire stellar population
neighbouring the BH follows a significantly shallower distribution than the BW
cusp. However, the distribution of the bound stars surrounding the BH is very
much consistent with a BW profile that extend for about an influence radius
$r<R_{\mathrm{inf}}$.

To visualise the radial disposition of the stars around the BH in our
simulations, we plot the radial stellar density for R06M300 at three different
times, displayed with the green thin line in Fig. 5. Within the influence
radius, this line is well fitted by a power-law $n\propto r^{-\beta}$ with
slope that oscillates around $\beta\sim 1.1$ (see middle panel of Fig. 6). On
the other hand, the distribution of the bound component (green thick line in
Fig. 5) is remarkably similar to the Bahcall-Wolf density profile. The
relaxation timescale at $R_{\mathrm{inf}}$ for R06M300 is very short:
gravitational scattering processes are expected to form a Bahcall-Wolf in less
than $\lesssim 5$ Myr999Here the relaxation timescale is estimated accounting
for both bound and unbound stars.. However, only the bound orbits are affected
by relaxation processes as they stay in the vicinity of the BH for a prolonged
time.
This result holds for R06M50 and R06M300, as shown in the left and the central
panels of Fig. 6, which illustrate that the power-law index $\beta$ for the
bound cloud (in blue) oscillates around the Bahcall-Wolf predicted value
throughout the simulation. Although the bound clouds in R06M50 and R06M300
consist of unequal mass stars, they are mainly composed of the heaviest stars
in the cluster (due to primordial mass segregation), which, as predicted by
Alexander & Hopman (2009), follow a BW density profile. This is a robust
result that does not depend on the initial half-mass radius of the cluster
(see Fig. 17).
On the other hand, the power-law index $\beta$ of R06M2000 bound cloud is
systematically lower than $1.75$ throughout the simulation, as shown in the
right panel of Fig. 6. R06M2000, having a larger influence radius than R06M50
and R06M300, also has a more extended bound cloud containing many more low-
mass stars. The latter tend to form a density profile shallower than the BW
cusp as indicated by Alexander & Hopman (2009).

###  4.3 Evolution of the cluster core

Figure 7: Time evolution of the central density (red) for the star cluster
model R04 without a central black hole. The system rapidly undergoes core
collapse and develops gravothermal oscillations driven by dynamically forming
binary stars. The same cluster model with a central black hole of 300
$\mathrm{M_{\odot}}$ (R04M300, blue) does not undergo core collapse but
expands gradually. The prevention of core collapse and the steady expansion is
driven by the already existing population of "binaries" which the central
black hole forms with many stars.

The BHs in our runs constantly occupy the very centre of the cluster.
Therefore, their mass growth is heavily influenced by the cluster core
properties (such as the density $\rho_{\mathrm{c}}$ and the velocity
dispersion $\sigma_{\mathrm{c}}$) and their evolution over time. At the same
time, the BHs, exerting violent dynamical interactions in their surroundings,
are expected to affect substantially the evolution of the cluster inner
region.

To understand how the central compact object impacts the evolution of the
cluster core, we rerun model R04M300 without a BH for 15 Myr. Fig. 7 compares
the evolution of the R04M300 core density with and without a BH. R04M300
without a BH (red line) undergoes gravothermal oscillations101010 Gravothermal
oscillations are well understood and have been studied extensively
theoretically. The first work that revealed such oscillations goes back to
Sugimoto & Bettwieser (1983). For a comprehensive overview of this topic see
Heggie (1993).: its core collapses and expands several times. Consequently,
the central density oscillates from $10^{8}$ up to $10^{11}$ particles per
pc3. The collapse phase is driven by two-body gravitational interactions,
which randomly scatter particles into the centre, causing the core to contract
(Spitzer, 1987). When the core is dense enough three-body, close interactions
can efficiently form hard binaries, which release energy into the core,
forcing it to expand. The expansion phase ends when the binaries escape the
centre due to violent dynamical interactions. Two-body interactions force the
core to re-collapse as soon as all the binaries are ejected away from the
centre, and the cycle restarts. Conversely, the core in R04M300 with a BH
(blue line in Fig. 7) experiences a monotonic expansion. In this simulation,
stars near the BH quickly bind to the compact object BH quickly bind to the
compact object and form BH-star binaries (many stars form a binary with the
central BH), which cause the core to expand. Due to the high BH mass, BH-stars
binaries are unlikely to be ejected from the core. Instead, they can provide a
steady energy flow to the centre until the end of the cluster evolution. In
other words, the central BH and the bound close stars act as a dynamical heat
source and force the core to a steadily expand (as already shown in previous
works Shapiro, 1977; Heggie et al., 2007).

Figure 8: Time evolution of the central number density (top panel) and core
radius (bottom panel) for the models as indicated in the legend. The initially
very dense systems expand rapidly and the central density drops by almost two
orders of magnitude. Stars bound to the central black holes drive the
expansion. This effect prevents the expected core collapse in the absence of a
central black hole (see Fig. 7). For the low mass (50 $\mathrm{M_{\odot}}$)
black hole simulation (orange) the onset of core collapse is seen but is
terminated once a bound nuclear subsystem has formed after a few million years
(see Fig. 4).

Fig. 8 illustrates that R06M2000, R06M300 and R08M300 register early core
expansion similar to R04M300. The expansion in these models starts right at
the beginning of the simulation. On the other hand, in the first few million
years, the core in R06M50 contracts (see orange line in Fig. 8). The
contraction halts after $\sim 4$ Myr, and the core expands until the
simulation ends. The $50\mathrm{M_{\odot}}$ BH in R06M50 cannot immediately
reverse the core contraction because, at $t=0$ Myr, it does not have a bound
cloud. The bound subsystem builds up in about $4$ Myr (orange line in Fig. 4)
right before the core starts to expand (orange line in Fig. 8). This fact
further indicates that the bound stars around the BH are responsible for the
dynamical heating of the cluster centre.

Figure 9: The number of tidal disruption events for the first 40 Myr of the
simulations. The direct disruption events are dominating in all simulations
(shaded regions) and the cicularised TDEs are sub-dominant. As explained in
the text circularised TDEs correspond to TCEs. They contribute between $9$ per
cent and up to $15$ per cent depending on the model. For comparison we show
the R06M50 simulation without a tidal drag force (blue) for which the number
of TDEs is reduced by about $10$ per cent. More massive central black holes
(R06M2000) as well as more compact clusters (R04M300) increase the number of
disruption events. Figure 10: The eccentricity distribution of the stellar
orbits right before tidal disruption for the M300 models with varying initial
half mass radii (left panel) and the R06 models with varying central black
hole masses (right panel). For all models most disruption orbits are very
radial with high eccentricities $e\gtrsim 0.99$ (these events are indicated in
the plot as direct TDEs). In a few cases, the distruption orbit is even
hyperbolic (unbound) with $e\gtrsim 1.0$ (not shown here). About $\sim 10$ per
cent of the stars have been tidally captured and have undergone partial or
total circularisation by drag forces before disruption (these events are
indicated as circularized TDEs and correspond, as explained in the text, to
TCEs). The size of the star cluster and the mass of the central black holes
only has a weak impact on the eccentricity distribution (see Section 2.4.)

###  4.4 Direct tidal disruption events and tidal capture events

Overall, the early expansion of the core registered in all our models does not
prevent the BH from experiencing a large number of TDEs. Fig. 9 shows that
only after 40 Myr two of our models (R04M300 and R06M2000) register more than
$\gtrsim 1200$ TDEs, and even our least dense system, R08M300, recorded more
than $\gtrsim 500$ TDEs in the same period of time.

We are interested in understanding how many of these TDEs were induced by
TCEs. Unfortunately, as we have modeled tidal interactions with a drag force,
directly integrated into the equations of motion, it is not possible to record
TCEs directly. Nevertheless, we can deduce indirectly how many of the TDEs
were triggered by TCEs by looking at the eccentricity of a star right before
its disruption. With the prescription adopted in this study, every captured
orbit undergoes a phase of circularisation followed by a phase of pericentre
shrinking (see Subsection 2.5). Thus, a tidally captured star always
experiences partial or total circularisation before being tidally disrupted.
On the contrary, direct TDEs tend to have eccentricities very close to unity.
The two panels in Fig. 10 display the eccentricity distribution of all TDEs.
The distributions are bimodal with a broad peak at $\log(1-e)\lesssim-2.0$ and
a narrow peak that extends between $\log(1-e)\sim-1.0$ and $\log(1-e)\sim
0.0$. Such plots provide a clear cut separation between direct TDEs ($e\sim
1.0$) and circularised TDEs ($e\lesssim 0.9$).

In principle, dynamical interactions could decrease the eccentricity of some
direct TDEs. However, in practice, this never occurs. Only the tidal
interaction drag force can significantly impact the eccentricity evolution of
these events. For instance, when we deactivate the tidal interaction drag
force, the model R06M50 registered only direct TDEs (no events with $e\lesssim
0.9$ see blue line in Fig. 10). With the drag force, R06M50 records 87
circularised TDEs. In summary, we can estimate the number of TCEs by counting
the number of circularised TDEs ($e\lesssim 0.9$).

Fig. 9 indicates that the BHs feed mainly on direct TDEs. However, TCEs have
also a relevant impact on the total number of TDEs: they trigger between $10$
per cent and $15$ per cent of the TDEs. Their contribution to the BH mass
growth is even more substantial. Contrary to a direct TDE, in a TCE, the
stellar material is tightly bound to the BH. Hence, the BH could absorb a
larger fraction of the stellar material. We discuss this in more detail in
Section 4.6.

Figure 11: Left panel: distribution of stellar masses $m_{*}$ tidally
disrupted (TDE) by the central BH for models R04M300, R06M300 and R08M300
(same initial BH mass and increasing half mass radii). Due to mass segregation
the most massive stars in the simulations are more likely to be disrupted.
Stars above $\sim 2$ $\mathrm{M_{\odot}}$ have formed in stellar mergers.
Right panel: distribution of stellar masses $m_{*}$ tidally disrupted by the
central BH for models R06M50, R06M300 and R06M2000 (same initial half mass
radius and increasing black hole masses). The distribution of R06M2000 looks
markedly flatter and does not peak at $\sim 1.7$ $\mathrm{M_{\odot}}$. This
model has about three times more TDEs than the other models and the most
massive stars are consumed in the first few million years. Thereafter
increasingly lower mass stars are disrupted by the central black hole.

###  4.5 Tidal disruption rate as a function of time

Figure 12: Semi-major axis vs. eccentricity about $10$ kyr before a star is
disrupted by the central black hole for models R06M50, R06M300, and R06M2000
(from top left clockwise). The R06M50 model without a drag force is also shown
for comparison at top right. The colored circles reported in each panel
represent all the stars that merged with the BH. The orbital elements are
computed before they could have been affected by tidal energy loss (here
indicated with the blue area). More than $90$ per cent of the colliding stars
were moving on a bound orbit.

Unsurprisingly, in our simulations, the average mass of the tidally disrupted
stars is significantly higher than the mean stellar mass; it varies from
$m_{*}\sim 1.3$ $\mathrm{M_{\odot}}$ up to $m_{*}\sim 1.6$
$\mathrm{M_{\odot}}$ depending on the initial conditions. In fact, due to mass
segregation, the BH preferentially devours the most massive stars in the
system. The mass distribution of these stars is shown in Fig. 11. For four of
our models, the distribution peaks at $2\mathrm{M_{\odot}}$, which are the
most massive stars of the clusters. R06M2000 is the only exception; it has a
rather flat mass distribution. The $2000\mathrm{M_{\odot}}$ BH destroys with
the same probability stars between $\sim 0.5\mathrm{M_{\odot}}$ and
$2.0\mathrm{M_{\odot}}$. This exception can be explained by looking at Fig.
12, which shows that most TDEs originate from bound orbits. Most of the stars
that contribute to the growth of the BH come from the bound cloud. In other
words, the BH feeds mainly on the bound subsystem. Due to its large BH mass,
R06M2000 is the model with the most extended bound cloud since bound clouds
extend for about an influence radius ($r<R_{\mathrm{inf}}$) (see Subsection
4.1). Therefore, compared to the other models, its bound subsystems encompass
many more low-mass stars. To recap, the reservoir of stars that sustain the
growth of the BH coincides mostly with the bound subsystem. Since R06M300 has
the most extended bound cloud, it has a higher probability of consuming
lighter stars.

Figure 13: Cumulative number of TDEs fitted using equation 10 for all the five
models. After calibrating the constant factor $F$ for each model, Eq. 10
predicts very accurately the number of TDEs as a function of time. The
histogram in the right bottom panel reports the distribution of the
$F$-values, as well as $F$ mean value $\pm$ the standard deviation. The
evolution of the cumulative number of TDEs is displayed up to 40 Myr for all
models, except for the R06M300 model which we evolved for $150\,\mathrm{Myr}$.

As we have already mentioned, the bound cloud plays a key role in driving most
of the TDEs. Equally important is the role played by the unbound component of
the cluster, which, first of all, acts as a reservoir for the bound subsystem,
ensuring that the fraction of bound stars is constant over time (see Fig. 4).
In addition, unbound stars constantly scatter with the bound orbits, and
through complex dynamical interactions, they also trigger TDEs and TCEs. Note
that orbits in the bound cloud are typically too wide to be affected by the
tidal drag force: most of the orbits are beyond the drag force region of
influence (blue area in Fig. 12 ). Gravitational scattering events are
necessary to lead stars close enough to be tidally captured or tidally
destroyed. In other words, the two-body relaxation process is indispensable
for driving the bound cloud stars to tidal disruption. The time scale at which
orbits near the BH change their energy and angular momentum is consistent with
the non-resonant relaxation time scale $t\mathrm{{}_{rel}}$. Since we include
post-Newtonian corrections in the equation of motion of the simulation
particles, resonant effects are suppressed in agreement with what previous
studies have found (Merritt et al., 2011; Brem et al., 2013). From these
considerations it follows that we can estimate the TDE rate using:

$\dot{N}_{\mathrm{TDE}}=F\frac{N_{b}}{t\mathrm{{}_{rel}}}$ (7)

where $N_{b}$ is the number of bound stars within the BH influence radius, $F$
is a numerical pre-factor that is calibrated for the different models and
$t\mathrm{{}_{rel}}$ is the relaxation time scale of the cluster calculated at
$R_{\mathrm{inf}}$ which we estimate following Spitzer (1987):

$t\mathrm{{}_{rel}}=\frac{1.8\times
10^{3}}{\ln(\Lambda)}\left(\frac{10^{7}\mathrm{M_{\odot}}/\mathrm{pc^{3}}}{\rho}\right)\left(\frac{\sigma}{100\mathrm{km/s}}\right)^{3}\left(\frac{1\mathrm{M_{\odot}}}{m_{*}}\right)\mathrm{Myr}$
(8)

where $\sigma$, $\rho$ and $m_{*}$ are respectively the velocity dispersion,
stellar density and the average stellar mass computed within
$r<R_{\mathrm{inf}}$. The Coulomb logarithm, $\Lambda$, can be approximated
with $\Lambda=0.11N$ (Giersz & Heggie, 1994). In our case $N$ is the number of
particles within the influence radius, therefore $N=2M\mathrm{{}_{BH}}/m_{*}$.
Similarly we can now rewrite $N_{b}=2f_{\mathrm{b}}M\mathrm{{}_{BH}}/m_{*}$,
where $f_{\mathrm{b}}$ is the fraction of bound stars at $r<R_{\mathrm{inf}}$.
It then follows that we can rewrite Eq. 7 as:

$\displaystyle\dot{N}_{\mathrm{TDE}}=$ $\displaystyle
1.1Ff_{\mathrm{b}}\ln\left(0.22\frac{M\mathrm{{}_{BH}}}{m_{*}}\right)\left(\frac{M\mathrm{{}_{BH}}}{10^{3}\mathrm{M_{\odot}}}\right)\left(\frac{\rho}{10^{7}\mathrm{M_{\odot}}/\mathrm{pc^{3}}}\right)$
(9)
$\displaystyle\times\left(\frac{100\mathrm{km/s}}{\sigma}\right)^{3}\mathrm{Myr^{-1}}\
$

where the stellar velocity dispersion $\sigma$, the stellar density $\rho$ and
the average stellar mass $m_{*}$ are computed at the influence radius
$(r=R_{\mathrm{inf}})$. To estimate the cumulative number of TDEs as a
function of time we integrate Eq. 9 obtaining:

$\displaystyle N_{\mathrm{TDE}}(t)=1.1F\int_{0}^{t}$ $\displaystyle
f_{\mathrm{b}}(t^{\prime})\left(\frac{M\mathrm{{}_{BH}}(t^{\prime})}{10^{3}\mathrm{M_{\odot}}}\right)\ln\left(0.22\frac{M\mathrm{{}_{BH}}(t^{\prime})}{m_{*}(t^{\prime})}\right)$
(10)
$\displaystyle\times\left(\frac{\rho(t^{\prime})}{10^{7}\mathrm{M_{\odot}}/\mathrm{pc^{3}}}\right)\left(\frac{100\mathrm{km/s}}{\sigma(t^{\prime})}\right)^{3}dt^{\prime}.\
$

Here, to stay as general as possible, we assumed that the variables
$f_{\mathrm{b}}$, $M\mathrm{{}_{BH}}$, $\rho$, $m_{*}$ and $\sigma$ functions
of time. Calibrating the variable $F$ on each model Eq. 10 predicts fairly
accurately the number of TDEs as a function of time as displayed in Fig. 13.
The bottom right panel of the same figure shows the best fitting value for $F$
in the different simulations and it indicates that $0.6\lesssim F\lesssim
1.0$.

Figure 14: The mass of the BH as a function of time for the simulation R06M300
(green line) and using equation 12 (black region). The black area perfectly
reproduces the simulation result until $150\,\mathrm{Myr}$ (the end of the
simulation). It then predicts $M\mathrm{{}_{BH}}$ for later times. The
uncertainty in the prediction is calculated including the fitting error for
$\rho$ and $m_{*}$ (see Subsection 4.6 for more details). Figure 15: Left
panel:Mass of the BH as a function of time computed by integrating Eq. 12.
When solving the equation we assume $m_{*}$ to be constant and equal to
$m_{*}=1.5\mathrm{M_{\odot}}$; We vary $f_{\mathrm{c}}$ from $1.0$ (full
accretion scenario, see dark blue region) to $0.2$ ($80$ per cent stellar
material is lost, see magenta region). Right panel: BH mass accretion as a
function of time for the five scenarios presented in the left panel. In one
scenario (red line) we limit the accretion to the Eddington rate and we assume
all the stellar materia in TDEs accumulate around the BH without mass loss.

###  4.6 Extrapolating the mass growth of the central black hole

This subsection describes a method to predict the BH mass growth beyond the
simulation time. We cannot evolve our systems for more than
$150\,\mathrm{Myr}$ as our simulations are computationally costly, but we can
use the results of our runs to extrapolate the BH mass for a longer time. We
start by connecting the tidal disruption rate $\dot{N}_{\mathrm{TDE}}$
estimated in equation 9 with $\dot{M}_{\mathrm{BH}}$ the BH accretion rate.
For this purpose, we need to estimate how much stellar mass remains bound to
the BH during a TDE and how much of the bound material ends up in the compact
object.

The simple argument presented in Rees (1988) indicates that if a BH destroys a
star in a parabolic orbit, half of the stellar mass is lost during the
disruption process. The remaining portion is expected to form an accretion
disk around the BH, eventually falling on the compact object on a time scale
that depends on the accretion process (Shiokawa et al., 2015; Bonnerot et al.,
2016; Hayasaki et al., 2016). If the accretion disk injects gas into the BH at
a rate larger than the Eddington rate, the disc likely becomes weakly bound
and a significant fraction of the gas is ejected through winds (Strubbe &
Quataert, 2011).

In our simulations, we model neither the accretion disc formation phase nor
the actual accretion phase. Instead we average out all these processes by
introducing the parameter $f_{\mathrm{c}}$, which represent the fraction of
stellar mass accreted by the BH during a TDE. It follows from the definitions
of $f_{\mathrm{c}}$ that:

$\displaystyle\dot{M}_{\mathrm{BH}}=f_{\mathrm{c}}m_{*}\dot{N}_{\mathrm{TDE}}.$
(11)

We can therefore estimate the BH mass growth rate $\dot{M}_{\mathrm{BH}}$ by
computing $\dot{N}_{\mathrm{TDE}}$ using Eq. 9. Observing that the velocity
dispersion to estimate $\dot{N}_{\mathrm{TDE}}$ is calculated at
$r=R_{\mathrm{inf}}$ we can further simplify the previous equations by
rewriting $\sigma$ as a function of the BH mass $M\mathrm{{}_{BH}}$ and the
density $\rho$ using the virial theorem (see Eq. 15). It thus follows:

$\dot{M}_{\mathrm{BH}}=Cf_{\mathrm{c}}f_{\mathrm{b}}\left(\frac{m_{*}}{\mathrm{M_{\odot}}}\right)\sqrt{\frac{\rho}{10^{7}\mathrm{M_{\odot}}/\mathrm{pc}^{3}}}\ln\left(0.22\frac{M\mathrm{{}_{BH}}}{m_{*}}\right)\frac{\mathrm{M_{\odot}}}{\mathrm{Myr}}$
(12)

where we have incorporated all the constants in $C=1.1\times 10^{6}F/C_{0}$
(see Appendix B for the definition of $C_{0}$). Since our simulations are
computationally very expensive, we evolved them for no more than
$150\,\mathrm{Myr}$. However, as long as we find a reliable extrapolation of
$f_{\mathrm{b}}$, $\rho$ and $m_{*}$, we can use Eq. 12 to predict the mass of
the black hole beyond our simulation time. It turns out that all three of
these quantities, the density, the average stellar mass and the fraction of
bound stars, manifest regular and predictable behaviour. As transpires from
the left panel of Fig. 16 it is reasonable to assume $f_{\mathrm{b}}$ does not
change over time. Hence we consider $f_{\mathrm{b}}\approx 0.2$ to be
constant. The density $\rho$ in the vicinity of the BH (upper panel in Fig. 8)
displays a consistent trend: after an initial phase of core expansion, all
models, independently of their initial BH mass and initial concentration,
converge to an almost identical evolution. Consequently, we can use a fitting
function to robustly extrapolate the later evolution of $\rho$, as we
illustrate comprehensively in Appendix B. In the same appendix, we also show
that the change in time of $m_{*}$ can be modelled using a linear fit (see the
left panel in Fig. 18).

To summarise, we provided fitting formulas to predict the evolution of $\rho$
and $m_{*}$, which, together with Eq. 12 allow us to extrapolate
$M\mathrm{{}_{BH}}$ as a function of time. Substituting in Eq. 12
$f_{\mathrm{c}}=0.5$, which is the fiducial value we assumed in our runs, we
obtain the results displayed in Fig. 14. The time evolution of the R06M300 BH
mass is perfectly consistent with the simulation results (green line in Fig.
14) reaching $2500\mathrm{M_{\odot}}$ in about $\sim 1$ Gyr.

Thanks to its simplicity, we can use Eq. 12 to explore systems that are
slightly different from the ones we simulated. For instance, a star cluster
with a larger number of particles ($N>10^{6}$) has a larger number of massive
stars at its disposal that the BH can consume. This implies that $m_{*}$ stays
constant for much longer. Simply substituting the star’s mass to a constant
value of $m_{*}=1.5\mathrm{M_{\odot}}$, leads to a more significant mass
growth: the BH in this scenario is expected to reach $3500\mathrm{M_{\odot}}$
in 1 Gyr (see the purple line in Fig. 15). By means of equation 12 we can also
explore how modifying $f_{\mathrm{c}}$ would affect our findings. We are not
only limited to changing the value of $f_{\mathrm{c}}$ but we can also make a
distinction between the accretion fraction for direct TDEs
($f_{\mathrm{c}}^{\mathrm{TDE}}$) and for TCEs
($f_{\mathrm{c}}^{\mathrm{TCE}}$). Since our models registered $\sim 90$ per
cent of direct TDEs and about $\sim 10$ per cent of TDEs triggered by tidal
capture. We can thus link $f_{\mathrm{c}}$ with
$f_{\mathrm{c}}^{\mathrm{TDE}}$ and $f_{\mathrm{c}}^{\mathrm{TCE}}$ as follow:

$f_{\mathrm{c}}=(0.9f_{\mathrm{c}}^{\mathrm{TDE}}+0.1f_{\mathrm{c}}^{\mathrm{TCE}})/2.$
(13)

We use the equality above to explore the scenario
$f_{\mathrm{c}}^{\mathrm{TCE}}=1.0,f_{\mathrm{c}}^{\mathrm{TDE}}=0.5$, which
is motivated by the fact that, in TDEs induced by tidal captures, stellar gas
is more tightly bound to the BH. The left panel of Fig. 15 outlines the main
results. Not surprisingly
$f_{\mathrm{c}}^{\mathrm{TDE}}=f_{\mathrm{c}}^{\mathrm{TCE}}=0.2$ fades the
growth of the BH, which reaches only $\sim 1500\mathrm{M_{\odot}}$ in 1 Gyr.
On the other hand, the scenario with
$f_{\mathrm{c}}^{\mathrm{TDE}}=f_{\mathrm{c}}^{\mathrm{TCE}}=1.0$ brings the
BH to $\sim 7600\mathrm{M_{\odot}}$. This scenario might require an initial
accretion rate of $\sim 15\dot{M}_{\mathrm{EDD}}$ as the dark blue line in the
right panel of Fig. 15 indicates. Observations indicate that super-Eddington
accretion can be possible during TDEs (see Komossa, 2015, and references
therein). However, theoretical studies show that jets and violent outflows
will probably unbind between $30$ per cent and up to $75$ per cent of the gas
around the BH (Ayal et al., 2000; Metzger & Stone, 2016; Toyouchi et al.,
2021). For this reason, the BH in our model might not be able to sustain an
accretion rate with $f_{\mathrm{c}}=1$, at least during the first $20-30$ Myr.
Nevertheless, the BH can still reach $\sim 7000\mathrm{M_{\odot}}$ in 1 Gyr,
even limiting the accretion rate to the Eddington accretion as shown (in red)
in the left panel of Fig. 15. If all the stellar gas generated in TDEs
accumulates around the BH and stays bound to it, the BH will accrete it within
1 Gyr, even restricting its growth rate to
$\dot{M}_{\mathrm{BH}}\lesssim\dot{M}_{\mathrm{EDD}}$. The accretion rate for
this scenario, indicated in the right panel of Fig. 15 in red, is constant
$\dot{M}_{\mathrm{BH}}=\dot{M}_{\mathrm{EDD}}$ for the first $\sim 100$ Myr;
as soon as the BH consumes all the gas accumulated during the past TDEs, the
accretion drops rapidly to significantly lower values. In conclusion, our
analysis indicates that $300\mathrm{M_{\odot}}$ BH can grow up to $\sim
10^{4}\mathrm{M_{\odot}}$ as long as, after TDEs, all stellar material remains
bound. This scenario might be unfeasible if the tidally destroyed stars move
in open orbits, but as we showed in the previous subsection that the BH feeds
primarily on the bound cloud. To put it in another way, the great majority of
the destroyed stars come from bound orbits. Therefore, most stellar material
formed during TDEs might have a high chance of remaining bound to the BH.

##  5 Discussion and conclusions

In this work, we evolved and analyzed five direct $N\mathrm{-body}$
simulations of compact ($R_{\mathrm{h}}\leq 0.8$ pc) star clusters
($N=256000$) for at most $150\,\mathrm{Myr}$ to investigate the growth of
massive BHs through repeated TCEs and TDEs. All our systems start with a
single central BHs ($50\mathrm{M_{\odot}}\leq M\mathrm{{}_{BH}}\leq
2000\mathrm{M_{\odot}}$) surrounded by mass-segregated low-mass stars in the
mass range $0.08\mathrm{M_{\odot}}\lesssim m_{*}\lesssim
2.0\mathrm{M_{\odot}}$. The clusters initial central stellar density
$n_{\mathrm{c}}>10^{7}$pc-3 and initial central velocity dispersion
($30\mathrm{km/s}\lesssim\sigma_{\mathrm{c}}\lesssim 90\mathrm{km/s}$) fulfill
the criteria provided by Stone et al. (2017) to trigger tidal capture runaway
collisions. The five models show several similarities in both the evolution of
the cluster core and the growth of the central black hole. The BHs in all
realizations, at the very beginning of the runs, give rise to a subsystem of
bound orbits in its vicinity. Around $20$ per cent to $30$ per cent of the
stars enclosed within the influence radius ($r<R_{\mathrm{inf}}$) stays bound
to the BH throughout the simulation. This cloud of bound stars forms a Bahcall
& Wolf (1976) cusp embedded in a reservoir of unbound particles with a
shallower density distribution. Our analysis indicate that the bound cloud
plays a key role in the evolution of the cluster core, as it constantly
releases energy to the core through dynamical interactions resulting in a
monotonic core expansion. Consequently, the clusters cannot mantain their
original central density which drop by about an order of magnitude in the
first few Myr. Despite this initial expansion, the BHs experienced a sustained
tidal disruption rate until the simulation ended. About $10-15$ per cent of
the TDEs are caused by tidal captures; in all these cases, the orbit of the
destroyed stars underwent partial or total circularisation.
We derived simple equations (see Eq. 9 and 10) that predict accurately the
number of TDEs experienced by the central BHs as shown in Fig. 13. We use
these equations to derive an expression that model the BH mass growth rate as
a function of time. With this expression we explore scenarios that slightly
deviate from our initial conditions and predict the time evolution of BH
masses up to 1 Gyr. We predict that a few hundred solar masses BH can easily
reach $\sim 2000-3000\mathrm{M_{\odot}}$ through repeated TDEs and TCEs within
1 Gyr. The BH can even grow up to $\sim 10^{4}\mathrm{M_{\odot}}$ if all the
stellar material remains bound to the BH and no mass loss occur during the
destruction and accretion phase (see Fig. 15). We argue that this scenario is
plausible because, as revealed by our simulations (see Fig. 12), the BHs feed
primarily on bound orbits. As long as the accretion proceeds at sub-Eddington
rate, all the stellar gas should remain bound to the BH. Numerous observations
indicate the presence of $10^{4}\mathrm{M_{\odot}}-10^{5}\mathrm{M_{\odot}}$
IMBHs in galactic nuclei. One of the most convincing observations of this type
reveals that the nuclear star cluster of the galaxy NGC 4395 hosts a BH with a
mass between $9\times 10^{3}\mathrm{M_{\odot}}$ (Woo et al., 2019) up to $\sim
4\times 10^{5}\mathrm{M_{\odot}}$ (Peterson et al., 2005; den Brok et al.,
2015). Our results show that tidal captures and tidal disruption events might
contribute significantly to the growth of the many observed IMBHs inside
nuclear star clusters. TCEs and TDEs might even be the most dominant IMBH
growth channel in clusters with moderate escape velocities ($<100$ km/s),
where the hierarchical formation scenario is suppressed by relativistic
recoils (Mapelli et al., 2021).
In this work, we did not intend to conduct realistic simulations of NSCs.
Instead, we developed idealised systems that contain all the key ingredients
to address the questions we aim to answer while being an approximate version
of real systems. There are two main simplifications adopted in this study.
First of all, we initialised our systems with an old population of stars. In
this way, we were able to neglect effects related to stellar evolution and
focus exclusively on dynamic effects. However, NSCs experience multiple
episodes of star formation that replenish the system of massive stars (Carson
et al., 2015; Kacharov et al., 2018). The latter might contribute
significantly to the BH mass growth and could dramatically alter the scenario
described in this work. In addition, our simulations contain only a single
central BH. In other words, we assume that after the BH subsystem evaporates,
only one BH remain in the cluster. However, the stellar system might keep a BH
binary instead of a single BH which might significantly change the tidal
disruption rate experienced by the compact objects.
In future studies, we plan to include in our initial conditions the central
BHs subcluster, since the latter could significantly influence the IMBH tidal
disruption rate (Teboul et al., 2022). We also plan to locate our systems at
the centre of an external potential to study the effect of the galaxy bulge
around a NSC. In addition, in subsequent work, we intend to make a
comprehensive comparison between our findings with analytical models for loss
cone dynamics (Stone & Metzger, 2016) and tidally-driven runaway growth Stone
et al. (2017).

## Acknowledgements

FPR, PHJ, SL and DI acknowledge the support by the European Research Council
via ERC Consolidator Grant KETJU (no. 818930). PHJ also acknowledges the
support of the Academy of Finland grant 339127. TN acknowledges support from
the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under
Germany’s Excellence Strategy – EXC-2094 – 390783311 from the DFG Cluster of
Excellence ’ORIGINS’. The simulations have been carried out with the Ampere
GPU system of the MPI for Astrophysics cluster "Freya" hosted by the Max
Planck Computing and Data facility. This study used also facilities hosted by
the CSC - IT Centre for Science, Finland.

##  Appendix A The properties of the bound cloud

Figure 16: Left panel: The same description as in Fig. 4. The fraction of
bound stars stays constant for the entire evolution of the model R06M300.
Right panel: The same description as in Fig. 6. The bound cloud density
profile in R06M300 is consistent with a BW density profile throughout the
simulation. Figure 17: The same description as in Fig. 6. Here we compare the
value of the slope parameter $\beta$ as a function of time for two models with
different initial half mass radii: R04M300 (left panel) and R08M300 (right
panel). The two models appear to have almost identical behaviour, and in both
cases, the bound subsystem density profile is consistent with the BW
distribution.

In this paper, we show that the central BHs form a cloud of bound stars within
their influence radius ($r<R_{\mathrm{inf}}$). Apart from an initial phase,
the fraction of bound stars tends to a constant value in all simulations, as
we see Fig. 4. For simplicity, in this figure, we show the evolution of the
fraction of the bound stars $f_{\mathrm{b}}$ only until $40$ Myr. Since we
evolve the model R06M300 for $150\,\mathrm{Myr}$ here (see left panel of Fig.
16), we report that $f_{\mathrm{b}}$ stay roughly constant also for longer
time periods. Moreover, in the main text, we indicate that stars in the bound
cloud assume a spatial distribution compatible with a Bahcall-Wolf density
profile. Also, in this case, for the sake of simplicity and clarity, we show
how the Bahcall-Wolf develops only for the models R06M50, R06M300 and R06M2000
for $t\leq 40$ Myr. In this appendix, we reveal that the cusp persists for a
time scale beyond $40$ Myr as displayed in the right panel of Fig. 16.
Moreover, we show in Fig. 17 that also models with different initial
concentrations develop a Bahcall-Wolf cusp.

##  Appendix B Fitting formulas for the cluster core evolution

Figure 18: Left panel: Average mass of bound particles enclosed within the
influence radius as a function of time (green line). Best linear fit with
parameters $a,b$ calibrated using the first 20 Myr of evolution (continuous
black line). Extrapolation of the linear fit until the end of the simulation
(black dashed line). Central panel: R06M300 density at $r=R_{\mathrm{inf}}$ as
a function of time (green line). Best fit using Eq. 14 with parameters
calibrated using the first 20 Myr of the evolution (continuous black line).
Extrapolation of the fit (black dashed line). Right panel: R06M300 velocity
dispersion within the influence radius as a function of time (green line).
Predicted velocity dispersion using equation 15.

In subsection 4.6, we describe how two-body relaxation processes drive bound
orbits to tidal disruption, and we estimate the rate of TDEs with
$\dot{N}_{\mathrm{TDE}}\propto\frac{N_{\mathrm{b}}}{t\mathrm{{}_{rel}}}$.
Using this relation, we derived the TDE rate as a function of $m_{*}$, $\rho$
and $\sigma$ calculated at the influence radius (see Eq. 9). In this appendix,
we predict the time evolution of these quantities for the model R06M300
employing fitting formulas, which we then use to extrapolate the mass growth
of the BH for about 1 Gyr.

The green line in the right panel of Fig. 18 indicates that $m_{*}$ decrease
linearly over time. We therefore use the fitting formula $m_{*}(t)=a+bt$. We
calibrate the parameters $a=1.38\mathrm{M_{\odot}}$ and $b=-8.4\times
10^{-4}\frac{\mathrm{M_{\odot}}}{\mathrm{Myr}}$ using the first 20 Myr of the
simulation R06M300. The black line in Fig. 18 (right panel) shows the fit,
while the dashed black line shows the extrapolation of the fit, which is also
in good agreement with the simulation. From empirical considerations we
realised that the evolution of $\rho$ as a function of time, can be modelled
using the following fitting function:

$\rho(t)=\frac{\rho_{0}}{[(t_{0}+t)/1\,\mathrm{Myr}]^{\frac{6}{5}}}$ (14)

The best parameters to fit the central density evolution of R06M300 are
$\rho_{0}=1.94\times 10^{8}\frac{\mathrm{M_{\odot}}}{\mathrm{pc^{3}}}$ and
$t_{0}=1.3$ Myr. Also, in this case, the parameters are found using the first
20 Myr of the simulation. The central panel in Fig. 18 indicates that equation
14 does an excellent job of extrapolating the simulation density for a longer
time. Regarding $\sigma$, we derive using the virial theorem that:

$\sigma=\left(C_{0}\left(\frac{M\mathrm{{}_{BH}}}{10^{3}\mathrm{M_{\odot}}}\right)\sqrt{\frac{\rho}{10^{7}\mathrm{M_{\odot}}/\mathrm{pc}^{3}}}\right)^{\frac{1}{3}}\mathrm{km/s}$
(15)

where $C_{0}$ is a dimensionless constant that depends on the distribution of
the stars within the influence radius. For the model R06M300 we find
$C_{0}=3.33\times 10^{4}$.

## References

  * Aarseth (2003) Aarseth S. J., 2003, Gravitational N-Body Simulations. Cambridge University Press
  * Abbott et al. (2016) Abbott B. P., et al., 2016, Phys. Rev. Lett., 116, 241103
  * Abbott et al. (2017) Abbott B. P., et al., 2017, Phys. Rev. Lett., 118, 221101
  * Alexander & Hopman (2009) Alexander T., Hopman C., 2009, ApJ, 697, 1861
  * Antonini et al. (2019) Antonini F., Gieles M., Gualandris A., 2019, MNRAS, 486, 5008
  * Arca-Sedda et al. (2021) Arca-Sedda M., Rizzuto F. P., Naab T., Ostriker J., Giersz M., Spurzem R., 2021, ApJ, 920, 128
  * Atallah et al. (2022) Atallah D., Trani A. A., Kremer K., Weatherford N. C., Fragione G., Spera M., Rasio F. A., 2022, arXiv e-prints, p. arXiv:2211.09670
  * Ayal et al. (2000) Ayal S., Livio M., Piran T., 2000, ApJ, 545, 772
  * Bahcall & Wolf (1976) Bahcall J. N., Wolf R. A., 1976, ApJ, 209, 214
  * Bahcall & Wolf (1977) Bahcall J. N., Wolf R. A., 1977, ApJ, 216, 883
  * Barth et al. (2004) Barth A. J., Ho L. C., Rutledge R. E., Sargent W. L. W., 2004, ApJ, 607, 90
  * Baumgardt et al. (2006) Baumgardt H., Hopman C., Portegies Zwart S., Makino J., 2006, MNRAS, 372, 467
  * Blanchet et al. (2006) Blanchet L., Buonanno A., Faye G., 2006, Phys. Rev. D, 74, 104034
  * Bonnerot & Lu (2020) Bonnerot C., Lu W., 2020, MNRAS, 495, 1374
  * Bonnerot et al. (2016) Bonnerot C., Rossi E. M., Lodato G., Price D. J., 2016, MNRAS, 455, 2253
  * Breen & Heggie (2013) Breen P. G., Heggie D. C., 2013, MNRAS, 432, 2779
  * Brem et al. (2013) Brem P., Amaro-Seoane P., Spurzem R., 2013, MNRAS, 434, 2999
  * Campanelli et al. (2007) Campanelli M., Lousto C. O., Zlochower Y., Merritt D., 2007, Physical Review Letters, 98, 231102
  * Carson et al. (2015) Carson D. J., Barth A. J., Seth A. C., den Brok M., Cappellari M., Greene J. E., Ho L. C., Neumayer N., 2015, AJ, 149, 170
  * Casares et al. (2014) Casares J., Negueruela I., Ribó M., Ribas I., Paredes J. M., Herrero A., Simón-Díaz S., 2014, Nature, 505, 378
  * Chin (1997) Chin S. A., 1997, Physics Letters A, 226, 344
  * Chin (2007) Chin S. A., 2007, arXiv e-prints, p. arXiv:0704.3273
  * Chin & Chen (2005) Chin S. A., Chen C. R., 2005, Celestial Mechanics and Dynamical Astronomy, 91, 301
  * Dai et al. (2018) Dai L., McKinney J. C., Roth N., Ramirez-Ruiz E., Miller M. C., 2018, ApJ, 859, L20
  * Dehnen & Hernandez (2017) Dehnen W., Hernandez D. M., 2017, MNRAS, 465, 1201
  * Di Carlo et al. (2021) Di Carlo U. N., et al., 2021, MNRAS, 507, 5132
  * El-Badry et al. (2022) El-Badry K., et al., 2022, MNRAS,
  * Event Horizon Telescope Collaboration et al. (2019) Event Horizon Telescope Collaboration et al., 2019, ApJ, 875, L1
  * Event Horizon Telescope Collaboration et al. (2022) Event Horizon Telescope Collaboration et al., 2022, ApJ, 930, L12
  * Fabian et al. (1975) Fabian A. C., Pringle J. E., Rees M. J., 1975, MNRAS, 172, 15
  * Farrell et al. (2014) Farrell S. A., et al., 2014, MNRAS, 437, 1208
  * Filippenko & Sargent (1989) Filippenko A. V., Sargent W. L. W., 1989, ApJ, 342, L11
  * Fragione et al. (2022) Fragione G., Kocsis B., Rasio F. A., Silk J., 2022, ApJ, 927, 231
  * Gerosa & Berti (2019) Gerosa D., Berti E., 2019, Phys. Rev. D, 100, 041301
  * Gieles et al. (2021) Gieles M., Erkal D., Antonini F., Balbinot E., Peñarrubia J., 2021, Nature Astronomy, 5, 957
  * Giersz & Heggie (1994) Giersz M., Heggie D. C., 1994, MNRAS, 268, 257
  * Giersz et al. (2015) Giersz M., Leigh N., Hypki A., Lützgendorf N., Askar A., 2015, MNRAS, 454, 3150
  * Giersz et al. (2019) Giersz M., Askar A., Wang L., Hypki A., Leveque A., Spurzem R., 2019, MNRAS, 487, 2412
  * Glebbeek et al. (2009) Glebbeek E., Gaburov E., de Mink S. E., Pols O. R., Portegies Zwart S. F., 2009, A&A, 497, 255
  * González Prieto et al. (2022) González Prieto E., Kremer K., Fragione G., Martinez M. A. S., Weatherford N. C., Zevin M., Rasio F. A., 2022, arXiv e-prints, p. arXiv:2208.07881
  * Greene et al. (2020) Greene J. E., Strader J., Ho L. C., 2020, ARA&A, 58, 257
  * Gültekin et al. (2022) Gültekin K., et al., 2022, MNRAS, 516, 6123
  * Hayasaki et al. (2016) Hayasaki K., Stone N., Loeb A., 2016, MNRAS, 461, 3760
  * Heggie (1993) Heggie D. C., 1993, arXiv e-prints, pp astro–ph/9312018
  * Heggie et al. (2007) Heggie D. C., Hut P., Mineshige S., Makino J., Baumgardt H., 2007, PASJ, 59, L11
  * Kacharov et al. (2018) Kacharov N., Neumayer N., Seth A. C., Cappellari M., McDermid R., Walcher C. J., Böker T., 2018, MNRAS, 480, 1973
  * King (1966) King I. R., 1966, AJ, 71, 64
  * Kochanek (1992) Kochanek C. S., 1992, ApJ, 385, 604
  * Komossa (2015) Komossa S., 2015, Journal of High Energy Astrophysics, 7, 148
  * Kormendy & Ho (2013) Kormendy J., Ho L. C., 2013, ARA&A, 51, 511
  * Kremer et al. (2020) Kremer K., et al., 2020, ApJ, 903, 45
  * Kremer et al. (2022) Kremer K., Lombardi J. C., Lu W., Piro A. L., Rasio F. A., 2022, ApJ, 933, 203
  * Kroupa (2001) Kroupa P., 2001, MNRAS, 322, 231
  * Kunth et al. (1987) Kunth D., Sargent W. L. W., Bothun G. D., 1987, AJ, 93, 29
  * Küpper et al. (2011) Küpper A. H. W., Maschberger T., Kroupa P., Baumgardt H., 2011, MNRAS, 417, 2300
  * Lee & Ostriker (1986a) Lee H. M., Ostriker J. P., 1986a, ApJ, 310, 176
  * Lee & Ostriker (1986b) Lee H. M., Ostriker J. P., 1986b, ApJ, 310, 176
  * Lin et al. (2018) Lin D., et al., 2018, Nature Astronomy, 2, 656
  * Lousto & Healy (2019) Lousto C. O., Healy J., 2019, Phys. Rev. D, 100, 104039
  * Lousto et al. (2010) Lousto C. O., Campanelli M., Zlochower Y., Nakano H., 2010, Classical and Quantum Gravity, 27, 114006
  * Lynden-Bell (1969) Lynden-Bell D., 1969, Nature, 223, 690
  * Magorrian et al. (1998) Magorrian J., et al., 1998, AJ, 115, 2285
  * Mapelli (2016) Mapelli M., 2016, MNRAS, 459, 3432
  * Mapelli et al. (2021) Mapelli M., Santoliquido F., Bouffanais Y., Arca Sedda M. A., Artale M. C., Ballone A., 2021, Symmetry, 13, 1678
  * Mardling (1995) Mardling R. A., 1995, ApJ, 450, 722
  * Mardling & Aarseth (2001) Mardling R. A., Aarseth S. J., 2001, MNRAS, 321, 398
  * Merritt et al. (2004) Merritt D., Piatek S., Portegies Zwart S., Hemsendorf M., 2004, ApJ, 608, L25
  * Merritt et al. (2011) Merritt D., Alexander T., Mikkola S., Will C. M., 2011, Phys. Rev. D, 84, 044024
  * Metzger & Stone (2016) Metzger B. D., Stone N. C., 2016, MNRAS, 461, 948
  * Neumayer et al. (2020) Neumayer N., Seth A., Böker T., 2020, A&ARv, 28, 4
  * Pechetti et al. (2022) Pechetti R., et al., 2022, ApJ, 924, 48
  * Pelupessy et al. (2012) Pelupessy F. I., Jänes J., Portegies Zwart S., 2012, New Astron., 17, 711
  * Peterson et al. (2005) Peterson B. M., et al., 2005, ApJ, 632, 799
  * Portegies Zwart & Meinen (1993) Portegies Zwart S. F., Meinen A. T., 1993, A&A, 280, 174
  * Portegies Zwart et al. (2004) Portegies Zwart S. F., Baumgardt H., Hut P., Makino J., McMillan S. L. W., 2004, Nature, 428, 724
  * Press & Teukolsky (1977) Press W. H., Teukolsky S. A., 1977, ApJ, 213, 183
  * Preto et al. (2004) Preto M., Merritt D., Spurzem R., 2004, ApJ, 613, L109
  * Rantala et al. (2020) Rantala A., Pihajoki P., Mannerkoski M., Johansson P. H., Naab T., 2020, MNRAS, 492, 4131
  * Rantala et al. (2021) Rantala A., Naab T., Springel V., 2021, MNRAS, 502, 5546
  * Rantala et al. (2022) Rantala A., Naab T., Rizzuto F. P., Mannerkoski M., Partmann C., Lautenschütz K., 2022, arXiv e-prints, p. arXiv:2210.02472
  * Ray et al. (1987) Ray A., Kembhavi A. K., Antia H. M., 1987, A&A, 184, 164
  * Rees (1988) Rees M. J., 1988, Nature, 333, 523
  * Reines (2022) Reines A. E., 2022, Nature Astronomy, 6, 26
  * Remillard & McClintock (2006) Remillard R. A., McClintock J. E., 2006, ARA&A, 44, 49
  * Rizzuto et al. (2021) Rizzuto F. P., et al., 2021, MNRAS, 501, 5257
  * Rizzuto et al. (2022) Rizzuto F. P., Naab T., Spurzem R., Arca-Sedda M., Giersz M., Ostriker J. P., Banerjee S., 2022, MNRAS, 512, 884
  * Samsing et al. (2018) Samsing J., Leigh N. W. C., Trani A. A., 2018, MNRAS, 481, 5436
  * Shapiro (1977) Shapiro S. L., 1977, ApJ, 217, 281
  * Shih et al. (2003) Shih D. C., Iwasawa K., Fabian A. C., 2003, MNRAS, 341, 973
  * Shiokawa et al. (2015) Shiokawa H., Krolik J. H., Cheng R. M., Piran T., Noble S. C., 2015, ApJ, 804, 85
  * Spitzer (1987) Spitzer L., 1987, Princeton University Press. IET
  * Steinberg & Stone (2022) Steinberg E., Stone N. C., 2022, arXiv e-prints, p. arXiv:2206.10641
  * Stone & Metzger (2016) Stone N. C., Metzger B. D., 2016, MNRAS, 455, 859
  * Stone et al. (2017) Stone N. C., Küpper A. H. W., Ostriker J. P., 2017, MNRAS, 467, 4180
  * Strubbe & Quataert (2011) Strubbe L. E., Quataert E., 2011, MNRAS, 415, 168
  * Sugimoto & Bettwieser (1983) Sugimoto D., Bettwieser E., 1983, MNRAS, 204, 19P
  * Teboul et al. (2022) Teboul O., Stone N. C., Ostriker J. P., 2022, arXiv e-prints, p. arXiv:2211.05858
  * The LIGO Scientific Collaboration et al. (2021) The LIGO Scientific Collaboration et al., 2021, arXiv e-prints, p. arXiv:2111.03606
  * Thorne & Hartle (1985) Thorne K. S., Hartle J. B., 1985, Phys. Rev. D, 31, 1815
  * Toyouchi et al. (2021) Toyouchi D., Inayoshi K., Hosokawa T., Kuiper R., 2021, ApJ, 907, 74
  * Volonteri et al. (2021) Volonteri M., Habouzit M., Colpi M., 2021, Nature Reviews Physics, 3, 732
  * Wang et al. (2020) Wang L., Iwasawa M., Nitadori K., Makino J., 2020, MNRAS, 497, 536
  * Weatherford et al. (2020) Weatherford N. C., Chatterjee S., Kremer K., Rasio F. A., 2020, ApJ, 898, 162
  * Wen et al. (2021) Wen S., Jonker P. G., Stone N. C., Zabludoff A. I., 2021, ApJ, 918, 46
  * Woo et al. (2019) Woo J.-H., Cho H., Gallo E., Hodges-Kluck E., Le H. A. N., Shin J., Son D., Horst J. C., 2019, Nature Astronomy, 3, 755
  * Yoshida (1990) Yoshida H., 1990, Physics Letters A, 150, 262
  * den Brok et al. (2015) den Brok M., et al., 2015, ApJ, 809, 101