# Growth of A Massive Black Hole In A Dense Star Cluster Via Tidal Disruption
Accretion

Seungjae Lee  Center for Theoretical Physics, Department of Physics and
Astronomy, Seoul National University, Seoul 08826, Korea  Ji-hoon Kim Ji-hoon
Kim<EMAIL_ADDRESS>Center for Theoretical Physics, Department of Physics
and Astronomy, Seoul National University, Seoul 08826, Korea  Seoul National
University Astronomy Research Center, Seoul 08826, Korea  Boon Kiat Oh  Center
for Theoretical Physics, Department of Physics and Astronomy, Seoul National
University, Seoul 08826, Korea

(Received November XX, 2022; Revised November XX, 2022; Accepted November XX,
2022)

###### Abstract

Stars that are tidally disrupted by the massive black hole (MBH) may
contribute significantly to the growth of the MBH, especially in dense nuclear
star clusters (NSCs). Yet, this tidal disruption accretion (TDA) of stars onto
the MBH has largely been overlooked compared to the gas accretion (GA) channel
in most numerical experiments until now. In this work, we implement a black
hole growth channel via TDA in the high-resolution adaptive mesh refinement
code Enzo to investigate its influence on a MBH seed’s early evolution. We
find that a MBH seed grows rapidly from $10^{3}\,{\rm M}_{\odot}$ to $\gtrsim
10^{6}\,\mathrm{M}_{\odot}$ in 200 Myrs in some of the tested simulations.
Compared to a MBH seed that grows only via GA, TDA can enhance the MBH’s
growth rate by up to more than an order of magnitude. However, as predicted,
TDA mainly helps the early growth of the MBH (from $10^{3-4}\,{\rm M}_{\odot}$
to $\lesssim 10^{5}\,\mathrm{M}_{\odot}$) while the later evolution is
generally dominated by GA. We also observe that the star formation near the
MBH is suppressed when TDA is most active, sometimes with a visible cavity in
gas (of size $\sim$ a few pc) created in the vicinity of the MBH. It is
because the MBH may grow expeditiously with both GA and TDA, and the massive
MBH could consume its neighboring gas faster than being replenished by gas
inflows. Our study demonstrates the need to consider different channels of
black hole accretion that may provide clues for the existence of supermassive
black holes at high redshifts.

galaxies: supermassive black holes – galaxies: kinematics and dynamics –
galaxies: formation – galaxies: star clusters – star clusters: general –
cosmology: theory – methods: numerical

##  1 Introduction

Discovered at the centers of most massive galaxies are the massive black holes
(MBHs) with masses $\gtrsim 10^{6}\,{\rm M}_{\odot}$. While the physical size
of a MBH is insignificant compared to that of its host galaxy, the dynamical
influence of the MBH extends throughout the entire host galaxy. In our
contemporary understanding, most of the MBH mass is believed to have been
acquired from the accreting gas (Bondi, 1952; Salpeter, 1955, 1964;
Zel’dovich, 1964; Lynden-Bell, 1969; Lynden-Bell & Rees, 1971; Lynden-Bell,
1978). However, the radiation from the accretion disk limits the gas accretion
rate (GAR) to the so-called Eddington rate, making it challenging to explain
the observed massive quasars at high redshifts (for reviews, see Inayoshi et
al., 2020). Thus, many studies have been conducted to understand the formation
and growth history of MBHs (for reviews, see Volonteri, 2010; Sesana, 2012;
Kormendy & Ho, 2013).

In hydrodynamics simulations that probe the formation and evolution of
galaxies, MBH physics has been a crucial component. Ever since the first
implementation of MBHs in a galaxy-scale simulation by Springel et al.
(2005a), simulators have shown that the energy from an accreting MBH is
essential in preventing the overcooling of gas (e.g., Sijacki et al., 2009;
Booth & Schaye, 2009; Dubois et al., 2010). They also have demonstrated that
the MBHs’ evolution is closely linked to the star formation of their host
galaxies (e.g., Springel et al., 2005a, b; Di Matteo et al., 2005; Hopkins et
al., 2006; Johansson et al., 2009; Kim et al., 2011, 2019; Choi et al., 2014).
The growth mechanism of the MBH itself has also been a topic of great
interest. For example, how the infalling gas overcomes the angular momentum
barrier has been discussed by many authors (e.g. Hopkins & Quataert, 2010,
2011; Emsellem et al., 2015). However, MBH physics in a galaxy-scale
hydrodynamics simulation is still far from complete. First and foremost, in
most previous numerical studies, MBHs grow only via the gas accretion (GA)
channel. And these GA models are often resolution-dependent (Booth & Schaye,
2009) and also rely on the assumption of spherically symmetric gas inflows in
many implementations.

Meanwhile, in the vicinity of a MBH seed at the center of a galaxy, a
significant fraction of mass exists in the form of stars, not just in gas.
When a star approaches a compact object such as a BH closer than its tidal
radius, the star is disrupted by the compact object’s tidal force, and a large
fraction of the resulting stellar debris eventually accretes to the compact
object (Rees, 1988; Strubbe & Quataert, 2009). A sufficient number of tidal
disruption events (TDEs) could help the BH to grow more rapidly than was
previously thought.

This may be especially true for the MBHs residing in the very dense nuclear
star clusters (NSCs) often found in the centers of galaxies. The typical
masses of NSCs are in the range $10^{4}\,{\rm
M}_{\odot}<M_{\mathrm{NSC}}<10^{9}\,{\rm M}_{\odot}$ while their central
stellar densities may exceed $10^{6}\,\mathrm{M}_{\odot}\mathrm{pc}^{-3}$ (for
reviews, see e.g., Neumayer et al., 2020). In addition, the NSCs are found to
coexist with MBHs in many galaxies (e.g., Lauer et al., 1998; Schödel et al.,
2018), and the correlations between NSCs and MBHs have been widely discussed
(e.g., Seth et al., 2008; Graham & Spitler, 2009; Nguyen et al., 2018). Some
authors have suggested that dense star clusters such as NSCs could be the
birthplaces of intermediate-mass black holes (IMBHs; $10^{2}\,{\rm
M}_{\odot}<M_{\mathrm{BH}}<10^{5}\,{\rm M}_{\odot}$). The gravothermal
contraction and the ensuing core collapse significantly increase the central
density of a star cluster (Aarseth et al., 1974; Giersz & Spurzem, 1994;
Takahashi, 1995), and then, runaway collisions between stars inside this dense
core could lead to the formation of an MBH seed of mass $\gtrsim 10^{2}\,{\rm
M}_{\odot}$ (Begelman & Rees, 1978; Quinlan & Shapiro, 1987; Ebisuzaki et al.,
2001; Portegies Zwart & McMillan, 2002; Gürkan et al., 2004; Portegies Zwart
et al., 2004; Freitag et al., 2006a, b; Miller & Davies, 2012). Even after the
MBH seed has formed, there are still numerous stars in its vicinity inside the
NSC. While the star-star collision must be nearly head-on for them to merge,
the now massive MBH seed may tidally disrupt stars even when they are not in a
head-on collision course. For example, combining observations and theories,
Wang & Merritt (2004) predicted that the TDEs are prevalent near a MBH, and
the tidal disruption rate (TDR) varies inversely with the MBH mass. Although
it is rather uncertain whether we can generalize their result to the MBHs of
masses $\lesssim 10^{5}\,{\rm M}_{\odot}$, their work showed that a young MBH
may have grown via successive TDEs. Using an analytical approach, Stone et al.
(2017) also argued that MBHs can form and grow inside NSCs.

Despite its importance, no numerical study has considered the tidal disruption
accretion (TDA) channel of MBH growth in a galaxy-scale (or star cluster-
scale) hydrodynamics simulation. One of the reasons is that the scattering
between a star and a MBH is not properly resolved in a typical galaxy-scale
hydrodynamics simulation — unless the spatial resolution is vastly increased
near the MBH. Indeed, given the current computational constraints, stellar
scatterings or TDEs can be explicitly resolved only with a direct summation
code such as $N$-body6++gpu (Wang et al., 2015). Yet, these direct $N$-body
codes tend not to include the hydrodynamics solver required to describe GA.
Therefore, a sub-resolution prescription for TDA in a high-resolution
hydrodynamic simulation can be the first step towards bridging the two
different numerical approaches.

Recently Pfister et al. (2021) devised a TDA model in the Ramses code and
studied the rate of TDEs in galaxy mergers. However, their work did not fully
investigate the MBH evolution itself. In addition, their TDA model computes
the TDR based on the stellar profile found in the simulation. Because a
simulation with insufficient resolution may not follow the scattering process
between stars and MBHs precisely, the stellar profile in the simulation may
give inaccurate TDR estimates (for more discussion, see Section 5.2).
Therefore, in the present study, we implement a TDA model based on a different
approach to mitigate the issue, while focusing on the growth of a MBH seed
after its formation in a NSC.

The paper is organized as follows. In Section 2, we present our TDA model:
theoretical backgrounds, and how we implement the TDA model in the adaptive
mesh refinement (AMR) hydrodynamics code Enzo. In Section 3, we introduce the
initial condition and the parameters of the simulations. The results of our
galaxy–NSC–MBH co-evolution simulations are described in Section 4, and we
discuss the limitations of our work and future studies in Section 5. In
Section 6, we conclude our paper with the key findings.

##  2 Tidal Disruption Accretion Model

###  2.1 Methodology Overview

If a star approaches a compact object such as a BH, a MBH, or a neutron star
closer than the tidal radius, the star is disrupted due to the tidal force of
the compact object. For a star with mass $m_{\star}$ and radius $r_{\star}$
approaching a MBH of mass $M_{\mathrm{BH}}$, we can estimate the tidal radius
$R_{\mathrm{T}}$ of the MBH as

$\displaystyle R_{\mathrm{T}}$
$\displaystyle=r_{\star}\left(\frac{\eta^{2}M_{\mathrm{BH}}}{m_{\star}}\right)^{\frac{1}{3}}$
$\displaystyle\simeq\,2.25\times
10^{-8}\eta^{\frac{2}{3}}\left(\frac{r_{\star}}{{\rm
R}_{\odot}}\right)\left(\frac{M_{\mathrm{BH}}}{m_{\star}}\right)^{\frac{1}{3}}\mathrm{pc},$
(1)

where $\eta$ is the Safronov number that depends on the inner structure of the
star, and ${\rm R}_{\odot}$ is the solar radius. After the TDE, a fraction $f$
of the resulting stellar debris accretes to the MBH ($f\lesssim 1.0$), while
the rest escapes the MBH’s gravity (Frank & Rees, 1976; Strubbe & Quataert,
2009).

In an ideal numerical simulation, one might imagine simply removing a star
particle that encounters the MBH particle within $R_{\mathrm{T}}$ and add the
star’s mass to the MBH’s. However, the length scale of $R_{\mathrm{T}}$ — and
the scattering process within — is too small to be resolved in a typical
galaxy-scale simulation that depicts the co-evolution of a MBH and its host
galaxy. Hence, we have adopted a method that computes the TDR based on a
statistical approach, also known as the loss cone theory (Frank & Rees, 1976;
Syer & Ulmer, 1999; Magorrian & Tremaine, 1999; Wang & Merritt, 2004). The
loss cone is the area in the phase space that is within or adjacent to the
tidal radius of the MBH. If a stellar orbit lies in the loss cone, the star is
thought to be a candidate for TD. Once we estimate the TDR,
$\dot{N}_{\mathrm{TD}}$, we obtain the tidal disruption accretion rate (TDAR)
of the MBH by

$\dot{M}_{\mathrm{BH,\,TDA}}=f\times m_{\star}\times\dot{N}_{\mathrm{TD}}$ (2)

assuming that the stars have identical masses $m_{\star}$. The rest of this
section describes how we estimate $\dot{N}_{\mathrm{TD}}$ in a simulation.

###  2.2 Key Assumptions

In many observational studies, the stellar density profiles at the galactic
centers are typically fitted by a power law $\rho_{\star}(r)\propto
r^{-\gamma}$ with a power-law index $\gamma$ (e.g., Faber et al., 1997; Lauer
et al., 1998; Ferrarese et al., 2006a, b; Schödel et al., 2009). For most
galaxies, $\gamma$ is between 1.0 and 2.0. Theoretically, this power law can
be understood as a feature of a self-consistent stellar system (Binney &
Tremaine, 2011). For example, a self-gravitating system following the Maxwell-
Boltzmann distribution with a constant velocity dispersion gives us
$\gamma=2.0$. The outer region of a collapsing single-mass Plummer model
yields $\gamma=2.23$ (Takahashi, 1995). Also, a stellar system dominated by a
single MBH follows the power-law profile with $\gamma=1.75$ (Bahcall & Wolf,
1976).

Based on these observations, when estimating the TDR we assume that the NSC
follows a power-law stellar density profile. However, unfortunately, it is not
desirable to completely trust the density profile in a typical galaxy-scale
simulation near the MBH since the stellar distribution there is under-
resolved. Therefore, we use a sub-resolution, analytic approach to infer the
stellar profile at the NSC’s center needed to compute the TDR. A naive, simple
power-law profile such as the one assumed in Wang & Merritt (2004) may work
well for high-mass MBHs ($M_{\mathrm{BH}}\gtrsim 10^{5}\,\,{\rm M}_{\odot}$),
but may over-estimate the TDR for low-mass MBHs ($M_{\mathrm{BH}}\lesssim
10^{4}\,\,{\rm M}_{\odot}$).111The Wang & Merritt (2004) model predicts that,
for a fixed power-law stellar density profile with $\gamma>1.42$ (see Eq.(4)
and the description below), the TDR varies inversely with the MBH mass.
Applying their model to the lower BH mass range ($M_{\mathrm{BH}}\lesssim
10^{4}\,\,{\rm M}_{\odot}$) returns an exceedingly high TDR
($\dot{N}_{\mathrm{TD}}\gg 10^{-2}\,\mathrm{yr}^{-1}$), which is a sign that
the TDR could be over-estimated with a naive assumption for the stellar
profile. Thus, instead, while we use the simple power-law stellar density
profile for high-mass MBHs, we consider a “cored” power law for low-mass MBHs,
as motivated by Stone & Ostriker (2015):

$\rho_{\star}\left(r\right)=\rho_{\mathrm{c}}\left(1+{r^{2}}/{r_{\mathrm{c}}^{2}}\right)^{-{\gamma\over
2}},$ (3)

where $\rho_{\mathrm{c}}$ and $r_{\mathrm{c}}$ are the density and size of the
core, respectively.222In practice, our model entails that a low-mass MBH sits
in a uniform density ($\rho_{\mathrm{c}}$) stellar core, while a high-mass MBH
is at the center of a simple power-law profile. Two more assumptions simplify
our model: (i) the masses of stars are identically $m_{\star}$ (value to be
chosen by the user; see Eq.(14)), and (ii) the MBH stays at the center of a
NSC without “wandering” (for discussion, see Appendix A).

###  2.3 Tidal Disruption Rate (TDR) Estimation

We now discuss how the TDR is calculated using the stellar density profile
assumed in Section 2.2, and the various stellar distribution properties found
in a simulation. In our prescription that computes the TDR, we first estimate
the average stellar density $\rho_{0}$ and the average 1D stellar velocity
dispersion $\sigma_{0}$ in a simulation inside a sphere of a user-defined
radius $R_{\mathrm{s}}$ from the MBH ($R_{\mathrm{s}}$ is set to
$0.5\,\mathrm{pc}$, a value close to the minimum cell size $\Delta
x_{\mathrm{min}}$ in a simulation; see Section 3.2). $r_{0}$ is then defined
as the radius from the MBH at which the power-law density profile and the
velocity dispersion profile becomes $\rho_{0}$ and $\sigma_{0}$, respectively.
Also, the power-law index $\gamma$ and the size of the stellar core,
$r_{\mathrm{c}}$, are manually assumed by the user (see Eq.(3) and Table 2),
but not the core density $\rho_{\mathrm{c}}$ which will later be estimated by
our TDA model.

Figure 1: Tidal disruption rate (TDR) as a function of massive black hole
(MBH) mass. The solid lines indicate the TDR of our model (see Section 2.3)
while the dashed lines are the estimates by Wang & Merritt (2004, i.e., no
core; $r_{\mathrm{c}}=0$). The left panel shows how the TDR varies with
$\sigma_{0}$ near MBH with $\gamma=2.0$:
$\sigma_{\mathrm{0}}=\mathrm{112}\,\mathrm{km\,s^{-1}}$
($r_{\mathrm{c}}=0.1\,\mathrm{pc}$),
$\sigma_{\mathrm{0}}=34\,\mathrm{km\,s^{-1}}$
($r_{\mathrm{c}}=0.03\,\mathrm{pc}$), and
$\sigma_{\mathrm{0}}=11\,\mathrm{km\,s^{-1}}$
($r_{\mathrm{c}}=0.01\,\mathrm{pc}$). The $r_{\mathrm{c}}$ values are chosen
in tandem with $\sigma_{0}$ to match the estimates in Figure 8 of Stone et al.
(2017) shown in this panel as gray dashed lines. The middle panel depicts the
TDR for different $\gamma=$ 1.5, 1.75, 2.0 with $r_{\mathrm{c}}=$ 0.1 pc and
$\sigma_{\mathrm{0}}=112\,\mathrm{km\,s^{-1}}$. Lastly, the right panel shows
the TDR for different $r_{\mathrm{c}}=$ 0.01, 0.3, 0.1 pc with $\gamma=2.0$
and $\sigma_{\mathrm{0}}=112\,\mathrm{km\,s^{-1}}$.

Then, the 1D velocity dispersion profile is written for a spherically
symmetric system as

$\displaystyle\sigma^{2}\left(r\right)$
$\displaystyle=\frac{G}{\rho_{\star}\left(r\right)}\int^{\infty}_{r}\frac{M\left(r^{\prime}\right)\rho_{\star}\left(r^{\prime}\right)}{r^{\prime
2}}dr^{\prime}$ (4) $\displaystyle=\frac{2\pi
G\rho_{0}r^{2}_{0}}{\left(3-\gamma\right)\left(\gamma-1\right)}\left(\frac{r}{r_{0}}\right)^{2-\gamma},$
(5)

where we used
$\rho_{\star}\left(r\right)\simeq\rho_{0}\left(r/r_{0}\right)^{-\gamma}$, an
approximation of Eq.(3) to simplify the integration (Binney & Tremaine,
2011).333This approximation does not make a significant difference in the
obtained $\sigma(r)$ or $\sigma_{\mathrm{c}}$ because the integration is
dominated by $r>r_{\mathrm{c}}$. Plugging $r=r_{0}$ into Eq.(5) yields

$\sigma_{0}^{2}=\frac{2\pi
G\rho_{0}r^{2}_{0}}{\left(3-\gamma\right)\left(\gamma-1\right)},$ (6)

from which, as expected, $r_{0}$ is determined by $\rho_{0}$ and $\sigma_{0}$
as

$\displaystyle r_{0}$
$\displaystyle=\sqrt{\frac{\left(3-\gamma\right)\left(\gamma-1\right)}{2\pi}\frac{\sigma_{0}^{2}}{G\rho_{0}}}$
$\displaystyle\simeq
2.0\,\,\mathrm{pc}\times\left(\frac{\sigma_{0}}{100\,\,\mathrm{km}\,\mathrm{s}^{-1}}\right)\left(\frac{\rho_{0}}{10^{5}\,\mathrm{M_{\odot}}\,\mathrm{pc}^{-3}}\right)^{-{1\over
2}}.$ (7)

In addition, plugging $r=r_{\mathrm{c}}$ into Eq.(5) gives the 1D stellar
velocity dispersion in the core as

$\sigma_{\mathrm{c}}^{2}=\frac{2\pi
G\rho_{0}r_{0}^{2}}{\left(3-\gamma\right)\left(\gamma-1\right)}\left(\frac{r_{\mathrm{c}}}{r_{0}}\right)^{2-\gamma},$
(8)

whereas, from Eq.(3), the core stellar density is

$\displaystyle\rho_{\mathrm{c}}=\rho_{\mathrm{0}}(1+r_{0}^{2}/r_{\mathrm{c}}^{2})^{\gamma\over
2}\,\,\simeq\,\,\rho_{0}\left(r_{\mathrm{0}}/r_{\mathrm{c}}\right)^{\gamma},$
(9)

where $r_{\mathrm{0}}^{2}/r_{\mathrm{c}}^{2}\gg 1$ because $r_{0}\sim
1\,\mathrm{pc}$ and $r_{\mathrm{c}}\sim[0.03,0.3]\,\mathrm{pc}$ are assumed in
our estimation (see Table 2 and Appendix B).

Table 1: Structural characteristics of a model galaxy in our fiducial initial condition | Density profile | Structural properties††footnotemark:  | Data type
---|---|---|---
Dark matter halo | Navarro et al. (1997) |  $M_{200}=1.074\times 10^{12}\,{\rm M}_{\odot}$, $v_{\mathrm{c,\,200}}=150\,\,{\rm km}\,{\rm s}^{-1}$, |  $10^{7}$ particles
|  |  $r_{200}=205.5\,\,{\rm kpc}$, $c=10$, $\lambda=0.04$ | $(m_{\mathrm{DM}}=1.791\times 10^{5}\,{\rm M}_{\odot})$
Stellar disk | exponential |  $M_{\mathrm{d,\,star}}=3.292\times 10^{10}\,{\rm M}_{\odot}$, |  $10^{7}$ particles
|  |  $r_{\mathrm{d,\,star}}=3.432\,\,{\rm kpc}$, $\,\,{z_{\mathrm{\,d,\,star}}}=0.1r_{\mathrm{d,\,star}}=343\,\,{\rm pc}$ | $(m_{\mathrm{d}}=3.292\times 10^{3}\,{\rm M}_{\odot})$
Gas disk | exponential |  $M_{\mathrm{d,\,gas}}=8.593\times 10^{9}\,{\rm M}_{\odot}$, | adaptive mesh
|  |  $r_{\mathrm{d,\,gas}}=r_{\mathrm{d,\,star}}=3.432\,\,{\rm kpc}$, $\,\,z_{\mathrm{\,d,\,gas}}=z_{\mathrm{d,\,star}}=343\,\,{\rm pc}$ |
Stellar bulge | Dehnen (1993) |  $M_{\mathrm{b}}=4.115\times 10^{9}\,{\rm M}_{\odot}$, |  $1.25\times 10^{6}$ particles
|  |  $r_{\mathrm{\,b}}=0.1r_{\mathrm{d,\,star}}=343\,\,{\rm pc}$, $\gamma_{\mathrm{\,b}}=2.0$, $r_{\mathrm{cutoff,\,b}}=50\,{\rm pc}$ | $(m_{\mathrm{b}}=3.292\times 10^{3}\,{\rm M}_{\odot})$
Nuclear star cluster (NSC) | Dehnen (1993) |  $M_{\mathrm{NSC}}=10^{7}\,{\rm M}_{\odot}$, |  $2\times 10^{4}$ particles
|  |  $r_{\mathrm{\,NSC}}=3.0\,{\rm pc}$, $\gamma_{\mathrm{\,\,NSC}}=2.0$, $r_{\mathrm{cutoff,\,NSC}}=0.5\,{\rm pc}$ | $(m_{\mathrm{NSC}}=500\,\,{\rm M}_{\odot})$
Massive black hole (MBH) | N/A | $M_{\mathrm{BH,\,init}}=10^{3}\,{\rm M}_{\odot}\,\,\,\,\,\,({\rm or}\,\,\,8\times 10^{3}\,{\rm M}_{\odot})$ | a single particle

Now we consider two cases:

  * •

If $M_{\mathrm{BH}}$ is small (i.e., $M_{\mathrm{BH}}<M_{\mathrm{c}}$; see
Eq.(12)), the gravitational influence of the MBH is restricted to the uniform
stellar core. In such a system, we adopt the TDR estimate in Rees (1988): rCl
˙N_Rees & ≃ 10^-4 yr^-1 (MBH106M⊙)^4 3
×(nc105pc-3) (σc100 km s-1)^-1, where
$n_{\mathrm{c}}=\rho_{\mathrm{c}}/m_{\star}$ is the stellar number density at
the core (the value of $m_{\star}$ to be chosen by the user; see Eq.(14)).
Here, $\dot{N}_{\mathrm{Rees}}$ increases with $M_{\mathrm{BH}}$.

  * •

If $M_{\mathrm{BH}}$ is large (i.e., $M_{\mathrm{BH}}>M_{\mathrm{c}}$; see
Eq.(12)), the gravitational influence of the MBH reaches beyond the core
radius. As stated in Section 2.2 we assume that the stellar density follows a
simple power law. Then, according to Wang & Merritt (2004), the TDR becomes

$\displaystyle\dot{N}_{\mathrm{WM}}$
$\displaystyle\approx\frac{M_{\star}\left(r_{\mathrm{crit}}\right)}{m_{\star}\,t_{\mathrm{\,R}}\left(r_{\mathrm{crit}}\right)}$
(10) $\displaystyle\simeq\left(3-\gamma\right)(\ln{\Lambda})\,G^{1\over
2}\rho_{0}^{2}r_{0}^{9\over 2}M_{\mathrm{BH}}^{-{3\over
2}}\left(\frac{r_{\mathrm{crit}}}{r_{0}}\right)^{{9\over 2}-2\gamma},$ (11)

where $M_{\star}(r)$ is the stellar mass enclosed in $r$,
$t_{\mathrm{\,R}}(r)$ is the stellar relaxation time at radius $r$, and
$\ln\Lambda=\ln(0.4\,{M_{\mathrm{BH}}}/{m_{\star}})$ is the Coulomb logarithm
(Spitzer & Harm, 1958).444Here in Eq.(11), the deformation of stellar
distribution by the central MBH (Bahcall & Wolf, 1976) is not considered.
$r_{\mathrm{crit}}$ is the radius at which a star is scattered through an
angular size of the loss cone in its dynamical time and is related to
$\rho_{0}$ and $r_{0}$ as rCl (rcritr0)^4-γ & = η2 3π(r⋆3m⋆)^13m_⋆^-1
×ρ_0^-1r_0^-4(lnΛ)^-1M_BH^73 , for which we set $\eta=0.844$ and
$r_{\star}={\rm R}_{\odot}\left({m_{\star}}/{{\rm M}_{\odot}}\right)^{1/3}$
for our simulations.555Two angles appear in the loss cone theory. First, the
“angular size” of the loss cone at radius $r$ is defined as
$\theta_{\mathrm{\,lc}}(r)=({R_{\mathrm{T}}}/{r^{2}})(GM_{\mathrm{BH}}/{\sigma}(r))$.
Recall that $R_{\mathrm{T}}$ is the tidal radius defined in Eq.(1). Second,
$\theta_{\mathrm{d}}(r)$ is the angle through which a star is scattered in a
dynamical time $t_{\mathrm{d}}(r)=r/\sigma(r)$. If
$\theta_{\mathrm{\,lc}}>\theta_{\mathrm{d}}$, the loss cone is in the
diffusive regime and becomes empty. If
$\theta_{\mathrm{\,lc}}<\theta_{\mathrm{d}}$, the loss cone is always filled
with stars. The radius at which these two angles are equal is called the
critical radius, $r_{\mathrm{crit}}$.,666Our estimate weakly depends on the
choice of $r_{\star}$. For example, inserting $r_{\star}={\rm
R}_{\odot}\left({m_{\star}}/{{\rm M}_{\odot}}\right)^{0.8}$ gives a rate
decreased by $\lesssim 5$%. Here, $\dot{N}_{\mathrm{WM}}$ decreases with
$M_{\mathrm{BH}}$ for $\gamma>1.42$ because $t_{\mathrm{R}}$ increases with
$M_{\mathrm{BH}}$ (i.e., refilling the loss cone becomes harder as the MBH
grows).

From these considerations we build a model for the TDR by combining Eqs.(• ‣
2.3) and (11):

$\dot{N}_{\mathrm{TD,1}}\left(M_{\mathrm{BH}},\rho_{0},\sigma_{0},\gamma,r_{\rm
c}\right)=\begin{cases}\dot{N}_{\mathrm{Rees}}&\text{if}\
M_{\mathrm{BH}}<M_{\mathrm{c}},\\\ \dot{N}_{\mathrm{WM}}&\text{if}\
M_{\mathrm{BH}}>M_{\mathrm{c}},\end{cases}$ (12)

where the transition mass $M_{\mathrm{c}}\,\,(\sim
10^{4-5}\,\mathrm{M}_{\odot})$ is the mass of a MBH that makes
$r_{\mathrm{crit}}(M_{\mathrm{BH}})$ equal to $r_{\mathrm{c}}$.777Note that
Eq.(4) tells us that $r_{\mathrm{crit}}$ is a function of $M_{\mathrm{BH}}$.
Because the Eq.(11) assumes that the MBH is embedded in a power-law profile,
it requires $r_{\mathrm{crit}}\gtrsim r_{\mathrm{c}}$. Therefore, in
principle, the transition between the two estimates, Eqs.(• ‣ 2.3) and (11),
should occur at $r_{\mathrm{crit}}=r_{\mathrm{c}}$. However, since the two
estimates, Eqs.(• ‣ 2.3) and (11), do not exactly agree at
$M_{\mathrm{BH}}=M_{\mathrm{c}}$, we limit Eq.(12) as

$\dot{N}_{\mathrm{TD}}\,\,=\,\,{\tt
min}\left\\{\dot{N}_{\mathrm{TD,1}},\,\,\dot{N}_{\mathrm{Rees}}({\scriptstyle\,M_{\mathrm{BH}}=M_{\mathrm{c}}}),\,\,\dot{N}_{\mathrm{WM}}({\scriptstyle
M_{\mathrm{BH}}=M_{\mathrm{c}}})\right\\}$ (13)

to remove the discontinuity at $M_{\mathrm{c}}$. By taking the minimum, we can
also be conservative in our TDR estimates.

In Figure 1, we plot the resulting
$\dot{N}_{\mathrm{TD}}\left(M_{\mathrm{BH}}\right)$ in order to observe its
dependence on various parameters. First, for low-mass MBHs
($M_{\mathrm{BH}}\lesssim 10^{4}\,\,{\rm M}_{\odot}$), our model predicts that
$\dot{N}_{\mathrm{TD}}$ increases with $M_{\mathrm{BH}}$ (solid lines in all
panels), while the extrapolation of the Wang & Merritt (2004) model with no
stellar core expects a decreasing $\dot{N}_{\mathrm{TD}}$ with respect to
$M_{\mathrm{BH}}$ (dashed lines in the left and middle panels). For high-mass
MBHs ($M_{\mathrm{BH}}\gtrsim 10^{5}\,\,{\rm M}_{\odot}$) in the left panel,
the larger $\sigma_{0}$ is, the higher $\rho_{\mathrm{c}}$ becomes, and so
does $\dot{N}_{\mathrm{TD}}$. Overall, this panel shows that our model behaves
similarly as Stone et al. (2017, see their Figure 8) seen here as gray dashed
lines. In the middle panel, readers can notice that the slope index of the
stellar profile, $\gamma$, determines the slope of $\dot{N}_{\mathrm{TD}}$ for
high-mass MBHs ($M_{\mathrm{BH}}\gtrsim 10^{5}\,\,{\rm M}_{\odot}$). The
larger $\gamma$ is, the higher $\rho_{\mathrm{c}}$ becomes for a given
$\sigma_{0}$, and so does the peak value of $\dot{N}_{\mathrm{TD}}$. Lastly,
in the right panel of Figure 1, varying $r_{\mathrm{c}}$ changes the peak
value of $\dot{N}_{\mathrm{TD}}$ and the transition mass $M_{\mathrm{c}}$
between the two MBH mass regimes, Eqs.(• ‣ 2.3) and (11). But it does not
change the slope of $\dot{N}_{\mathrm{TD}}$ for the high-mass MBHs because,
for the chosen $\gamma=2.0$, Eq.(5) predicts a constant velocity dispersion.

###  2.4 Implementing the Tidal Disruption Accretion (TDA) Channel in the
Simulation Code

Once the stellar density $\rho_{0}$ and dispersion $\sigma_{0}$ around the MBH
are specified, we find the TDR with Eq.(13). Then, from Eq.(2) we increase the
MBH mass at each timestep by

$\Delta M_{\mathrm{BH,\,TDA}}=fm_{\star}\,\dot{N}_{\mathrm{TD}}\,\Delta\,t$
(14)

where $\Delta\,t$ is the size of the timestep. The mass of a disrupted star in
our model, $m_{\star}$, is set to $0.7\,\mathrm{M}_{\odot}$, the average value
of the Salpeter IMF (Salpeter, 1955). However, it should be noted that, due to
limited resolution, a newly-formed “star particle” in our simulation has a
mass of $m_{\rm star}\gtrsim 200\,M_{\odot}$. Therefore, in practice, we
subtract $\Delta M_{\mathrm{BH,\,TD}}$ uniformly from all the “star particles”
within $R_{\mathrm{s}}$ from the MBH ($R_{\mathrm{s}}=0.5\,\mathrm{pc}$; see
Section 2.3), and then add it to the MBH. In other words, each star particle
within $R_{\mathrm{s}}$ loses a fraction $\alpha$ of its mass as

$m_{\mathrm{star}}\,\,\rightarrow\,\,\,\,(1-\alpha)\,m_{\mathrm{star}}=\left(1-\frac{m_{\star}\,\dot{N}_{\mathrm{TD}}\,\Delta\,t}{M_{\star}(R_{\mathrm{s}})}\right)\,m_{\mathrm{star}}$
(15)

in which we set $f=1.0$ in our simulations for simplicity.

##  3 Simulations

In this section, we present the initial conditions, refinement criteria, and
baryonic physics employed in our suite of simulations. The high-resolution
Eulerian AMR code Enzo (Bryan et al., 2014; Brummel-Smith et al., 2019)
provides the key zoom-in simulation technology and the baseline physics
modules, and here we describe the ones that are closely related to the topic
of our interest.

###  3.1 Initial Condition

For a fiducial initial condition, we adopt an isolated Milky Way-mass galaxy
harboring a MBH and a NSC at its center. While we generate our initial
condition using the code Dice (Perret, 2016), many galactic properties are
motivated by one of the isolated galaxy initial conditions in the AGORA High-
resolution Galaxy Simulations Comparison Project (Kim et al., 2014, 2016;
Roca-Fàbrega et al., 2020, 2021). The galaxy in our initial condition includes
a dark matter halo, an exponential stellar/gas disk, a stellar bulge, a NSC,
and a MBH seed (summarized in Table 1).

The dark matter halo follows the Navarro-Frank-White profile (NFW; Navarro et
al., 1997) with the circular velocity
$v_{\mathrm{c,\,200}}=150\,\,\mathrm{km}\,\mathrm{s}^{-1}$ and the virial mass
$M_{200}=1.072\times 10^{12}\,{\rm M}_{\odot}$. The concentration parameter
$c$ and the spin parameter $\lambda$ are set to 10 and 0.04, respectively. The
dark matter halo is composed of $10^{7}$ particles with equal masses
$m_{\mathrm{DM}}=1.791\times 10^{5}\,{\rm M}_{\odot}$. The masses of the
stellar and gas disk are $M_{\mathrm{d,\,star}}=3.292\times 10^{10}\,{\rm
M}_{\odot}$ and $M_{\mathrm{d,\,gas}}=8.593\times 10^{9}\,{\rm M}_{\odot}$,
respectively. Each follows an exponential profile, e.g.,

$\rho_{\rm
d,\,star}(r,z)=\rho_{\mathrm{0,\,star}}\,\exp({-{r}/{r_{\mathrm{\,d,\,star}}}})\exp({-{\left|z\right|}/{z_{\mathrm{\,d,\,star}}}})$
(16)

with the scale radius $r_{\mathrm{d,\,star}}=3.432\,\mathrm{kpc}$, the scale
height $z_{\mathrm{d,\,star}}=0.1\,r_{\mathrm{d,\,star}}$, and
$\rho_{\mathrm{0,\,star}}={M_{\mathrm{d,\,star}}}/{(4\pi
r_{\mathrm{d,\,star}}^{2}\,z_{\mathrm{d,\,star}})}$. The stellar bulge of mass
$M_{\mathrm{b}}=4.297\times 10^{9}\,{\rm M}_{\odot}$ follows the Dehnen
profile (Dehnen, 1993; Tremaine et al., 1994),

$\rho_{\mathrm{b}}\left(r\right)=\frac{\left(3-\gamma_{\mathrm{\,b}}\right)M_{\mathrm{b}}}{4\pi}\frac{r_{\mathrm{\,b}}}{r^{\gamma_{\mathrm{\,b}}}\left(r+r_{\mathrm{\,b}}\right)^{4-\gamma_{\mathrm{\,b}}}}$
(17)

with a half-mass radius $r_{\mathrm{\,b}}=0.1\,r_{\mathrm{d,\,star}}$ and the
index $\gamma_{\mathrm{\,b}}=2.0$ (can be regarded as a generalized version of
Hernquist, 1990).888In contrast, the initial condition in the AGORA Project
(Kim et al., 2016) follows the Hernquist profile which gives an infinite
stellar density at $r=0$. In the initial condition, each “star particle” in
the disk and bulge has a mass of $3.292\times 10^{3}\,{\rm M}_{\odot}$.

To set up a NSC with its own structural properties, separately from the
stellar bulge, we cap the density of the central region of the bulge by
setting a “cutoff” radius, $r_{\mathrm{cutoff,\,b}}=50\,\,{\rm pc}$ — i.e.,
$\rho_{\mathrm{b}}(r<r_{\mathrm{cutoff,\,b}})=\rho_{\mathrm{b}}(r_{\mathrm{cutoff,\,b}})$.
Then, a NSC of mass $M_{\mathrm{NSC}}=10^{7}\mathrm{M}_{\odot}$ is inserted at
the galaxy’s center. It follows the same Dehnen profile but with a smaller
half-mass radius $r_{\mathrm{\,NSC}}=3.0\,\,{\rm pc}$, a smaller cutoff radius
$r_{\mathrm{cutoff,\,NSC}}=0.5\,\,{\rm pc}$, and
$\gamma_{\mathrm{\,\,NSC}}=2.0$.999It should noted that the power-law profile
with an index $\gamma_{\mathrm{\,\,NSC}}=2.0$ is used to merely initialize the
“star particles” within the NSC. The particles’ distribution later on may not
follow the initial profile, and is likely unrealistic anyways because the
central region of the NSC is under-resolved. In our sub-resolution TDR
estimates, we “assume” that the “individual stars” in the vicinity of the MBH
follow the power-law profile, Eq.(3), with $\gamma=1.5-2.0$ as listed in Table
2. See Section 5.2 for more discussion.  Each “star particle” in the NSC
represents $500\,\,{\rm M}_{\odot}$, a value smaller than the initial MBH
mass. Finally, we plant a MBH at the NSC’s center — with two choices of
initial masses, $M_{\rm BH,\,init}=10^{3}$ or $8\times 10^{3}\,{\rm
M}_{\odot}$ (for motivations for its value, see Appendix A). The MBH physics
implemented is described in Sections 2.4 and 3.4.

###  3.2 Refinement Strategy

In the simulation box of $(1.307\,\mathrm{Mpc})^{3}$, we apply AMR only in the
pre-defined, innermost $(200\,\,{\rm pc})^{3}$ box ($\sim$ size of the stellar
bulge), dubbed RefineRegion. On the outside, the RefineRegion is surrounded by
static nested volumes with successively coarser resolutions. The RefineRegion
is first uniformly resolved with cells of $5\,{\rm pc}$ width. Each of the
cells is then adaptively split into $2^{3}$ child cells if the cell’s stellar
or gas mass exceeds a certain threshold which depends on the cell size $\Delta
x$ (or the refinement level) as

$\begin{cases}\,\,M_{\mathrm{ref,\,gas}}^{\Delta x}=\sqrt{2^{3}\left(\Delta
x\over\,\,5\,{\rm pc}\,\,\right)}\times 5\times 10^{2}\,{\rm
M}_{\odot}\,,\\\\[5.69054pt] \,\,M_{\mathrm{ref,\,star}}^{\Delta
x}=\sqrt{2^{3}\left(\Delta x\over\,\,5\,{\rm pc}\,\,\right)}\times 5\times
10^{4}\,{\rm M}_{\odot}\,.\end{cases}$ (18)

Therefore, the mass threshold for refinement decreases from
$M_{\mathrm{ref,\,gas}}^{\Delta x=5\,{\rm pc}}=2^{3/2}\times 5\times
10^{2}\,{\rm M}_{\odot}$ at $\Delta x=5\,{\rm pc}$ down to
$M_{\mathrm{ref,\,gas}}^{\Delta x=0.625\,{\rm pc}}=5\times 10^{2}\,{\rm
M}_{\odot}$ at $\Delta x=0.625\,\,{\rm pc}=\Delta x_{\mathrm{min}}$ when the
refinement stops. This is a super-Lagrangian refinement scheme to adaptively
increase spatial resolution in the targeted region around the NSC at limited
computational cost.

Figure 2: The distribution of the masses of newly-formed “star particles” in our simulations. The black dashed line is the lowest possible mass of a “star particle”, $m_{\mathrm{star},\,\mathrm{min}}=235\,\,\mathrm{M}_{\odot}$ set by the star particle formation criteria (after it loses a fraction of its mass via stellar feedback). In a simulation without the tidal disruption accretion (TDA) of stars (blue), $m_{\mathrm{star},\,\mathrm{min}}$ is indeed the lowest mass in the distribution. In contrast, in a simulation with TDA (orange), star particle masses can be below $m_{\mathrm{star},\,\mathrm{min}}$ as they may lose a fraction of their masses via TDA. For detailed explanations, see Section 3.3. Table 2: A suite of idealized galaxy–NSC–MBH co-evolution simulations listed with their runtime parameters | Run name | MBH seed mass | Tidal disruption accretion (TDA)††footnotemark:  | NSC power-law index assumed | NSC core size assumed
---|---|---|---|---|---
|  | $M_{\mathrm{BH,\,init}}\,\,\left[\mathrm{M}_{\odot}\right]$ | [on/off] | $\gamma$ | $r_{\mathrm{c}}\,\,\left[\mathrm{pc}\right]$
Set A‡‡footnotemark:  | NoTDA1 | $8\times 10^{3}$ | off (i.e., only gas accretion) | - | -
gamma1.5 | $8\times 10^{3}$ | on (i.e., gas accretion$+$TDA) | 1.5 | 0.1
gamma1.75 | $8\times 10^{3}$ | on | 1.75 | 0.1
gamma2.0 | $8\times 10^{3}$ | on | 2.0 | 0.1
Set B | NoTDA2 | $10^{3}$ | off | - | -
rc0.03 | $10^{3}$ | on | 2.0 | 0.03
rc0.05 | $10^{3}$ | on | 2.0 | 0.05
rc0.1 | $10^{3}$ | on | 2.0 | 0.1
rc0.3 | $10^{3}$ | on | 2.0 | 0.3

###  3.3 Star Particle Formation and Feedback

A “star particle” forms if the density of a maximally refined cell exceeds a
certain threshold $\rho_{\mathrm{th}}$, and the cell cannot be reliably
treated in the hydrodynamics solver. Our star particle formation model is
based on Cen & Ostriker (1992) and Kim et al. (2011). Following the Jeans
argument, we choose the threshold $\rho_{\mathrm{th}}={\pi
c_{\mathrm{s}}^{2}}/{(G\lambda_{\mathrm{J}}^{2})}$ with the Jeans length set
to the finest spatial resolution, $\lambda_{\mathrm{J}}=\Delta
x_{\mathrm{min}}$, and with the sound speed $c_{\mathrm{s}}$ at
$100\,\mathrm{K}$. This gives
$\rho_{\mathrm{th}}=2571\,\,\mathrm{M}_{\odot}\mathrm{pc}^{-3}$ (or the
threshold mass $M_{\mathrm{J}}=628\,\,\mathrm{M}_{\odot}$ for the given
$\Delta x_{\mathrm{min}}=0.625\,\,{\rm pc}$). We turn 50% of the cell’s mass
into a newly-formed “star particle” with its dynamical timescale set to
$\tau_{\mathrm{dyn}}=1\,\mathrm{Myr}$. In the next $12\tau_{\mathrm{dyn}}$,
the star particle contributes to the thermal supernovae feedback while
returning 25 % of its mass back to gas phase. These considerations yield the
lowest possible star particle mass of
$m_{\mathrm{star,\,min}}=235\,\,\mathrm{M}_{\odot}$.

Figure 2 shows the mass distribution of newly-formed star particles. In the
baseline simulation without TDA (blue histogram for the NoTDA1 run; see Table
2 and Section 2.4), $m_{\mathrm{star},\,\mathrm{min}}$ is indeed the lowest
mass in the distribution. In contrast, in the run with TDA (orange histogram
for the rc0.03 run), the masses of star particles can be lower than
$m_{\mathrm{star},\,\mathrm{min}}$ because star particles near the MBH may
lose their masses to the MBH via the TDA channel (see Section
2.4).101010Therefore, the present stellar masses may not be suitable when
estimating the star formation rate (SFR). Instead, since the number of newly-
formed star particles is preserved, we estimate the SFR by counting all the
newly-formed star particles, and assuming that their masses at the time of
their formation are equally
$m_{\mathrm{star,\,min}}=235\,\,\mathrm{M}_{\odot}$ (see Section 4.4).

###  3.4 Massive Black Hole (MBH) Accretion and Feedback

The MBH particle inserted at the center of the simulation box (see Section
3.1) grows via two channels: GA and TDA. Then the total MBH accretion rate is
simply

$\dot{M}_{\mathrm{BH}}=\dot{M}_{\mathrm{BH,\,TDA}}+\dot{M}_{\mathrm{BH,\,GA}}$
(19)

where the first term (TDAR) is from Eq.(2) or (14) (see Section 2), and the
second term (GAR) is from the conventional Bondi-Hoyle-Lyttleton formalism
(Bondi & Hoyle, 1944; Bondi, 1952) with the Eddington limit,

$\displaystyle\dot{M}_{\mathrm{BH,\,GA}}$
$\displaystyle=\min\left(\dot{M}_{\mathrm{BH,\,Bondi}},\,\,\dot{M}_{\mathrm{BH,\,Edd}}\right)$
$\displaystyle=\min\left(\frac{4\pi
G^{2}M_{\mathrm{BH}}^{2}\rho_{\mathrm{B}}}{c_{\mathrm{s}}^{3}},\,\,\frac{4\pi
GM_{\mathrm{BH}}m_{\mathrm{p}}}{\epsilon_{\mathrm{r}}\sigma_{\mathrm{T}}c}\right),$
(20)

for a MBH residing in a cell with the sound speed $c_{\mathrm{s}}$, where
$m_{\mathrm{p}}$ is the proton mass, $\sigma_{\mathrm{T}}$ is the Thomson
scattering cross section, and $\epsilon_{\mathrm{r}}=0.1$ is the BH’s
radiative efficiency. $\rho_{\mathrm{B}}$ is the gas density at the Bondi
radius $R_{\mathrm{B}}={2GM_{\mathrm{BH}}}/{c_{\mathrm{s}}^{2}}$, and is
estimated from the density $\rho_{\rm gas}$ of the cell where the MBH particle
resides by $\rho_{\rm B}=\rho_{\rm gas}\,\cdot\,{\tt min}\\{(\Delta x/R_{\rm
B})^{1.5},1.0\\}$ (for details, see Kim et al., 2011, 2019). The MBH particle
returns thermal feedback energy to the cells around it at a rate of

$L_{\mathrm{BH}}=\epsilon_{\mathrm{r}}\,\dot{M}_{\mathrm{BH,\,GA}}\cdot c^{2}$
(21)

which is proportional only to the GAR, $\dot{M}_{\mathrm{BH,\,GA}}$.

##  4 Results

In this section, we analyze the effects of our TDA model on the growth of the
MBH using the simulations with varying TDA parameters. We also examine the
secondary effects arising from the TDA-boosted MBH accretion.

Figure 3: Snapshots of the NoTDA2 run at $t$ = 50 Myr. The right panel is the
face-on gas density of the central (5 kpc)2 sliced through the location of the
MBH. Among the three zoomed-in images on the left (250 pc width), the top two
panels are the edge-on and face-on sliced gas density centered on the MBH,
while the bottom panel is the projected stellar density. The black dot in each
image indicates the MBH’s position. The figure illustrates the high resolution
we retain near the MBH to apply the TDA model. For more information, see
Sections 3.1 and 4.1.

###  4.1 Simulation Suite and Initial Relaxation

We have performed a suite of simulations to investigate the galaxy–NSC–MBH co-
evolution and the effects of our TDA model on the growth of the MBH (as listed
in Table 2). As discussed in Section 2.3, the two prominent parameters that
the user needs to select for the sub-resolution prescription of the TDA
channel are: (i) the power-law index $\gamma$ of the NSC’s density profile,
and (ii) the size of the NSC’s stellar core, $r_{\mathrm{c}}$. The group of
runs with varying $\gamma$ (from the gamma1.5 run to the gamma2.0 run) and the
NoTDA1 run is called “Set A”. Another group of runs with varying
$r_{\mathrm{c}}$ (from the rc0.03 run to the rc0.3 run; for discussions on our
$r_{\mathrm{c}}$ choices, see Appendix B) and the NoTDA2 run is called “Set
B”.

We first discuss the initial relaxation of the simulation in the first few
Myrs using the fiducial NoTDA2 run (see Table 2). The gas density and the
stellar distribution in the NoTDA2 run at $t=50\,{\rm Myr}$ are shown in
Figure 3. Due to the relaxation of the initial density distribution, the
stellar density at the galactic center changes in the first $\sim 5\,{\rm
Myr}$. Figure 4 shows the stellar density (top) and the enclosed mass profile
(bottom) centered on the NSC in the NoTDA2 run at several epochs. An increase
in the stellar density between $t=0$ and 5 Myr is noticeable, while the inner
density slope changes from $\gamma_{\mathrm{\,\,NSC}}=2.0\,$ to $\,\sim 1.0$.
However, the profiles do not change significantly after 5 Myr. The changes in
the first few Myrs are not because of any astrophysical origin, but because
our idealized NSC and stellar bulge in the initial condition follow the
artificial analytic fits that are prone to further collapse. Considering that
there is no easy way to initialize a realistic, relaxed galaxy with all its
constituents in Table 1, we may regard the $t=5$ Myr galaxy as our de facto
initial condition.

After the initial relaxation, the central stellar density inside the NSC
($\lesssim 3$ pc from the MBH) is
$\sim$$10^{5}\,\mathrm{M}_{\odot}\mathrm{pc}^{-3}$, while the 1D velocity
dispersion is $\sim$$150\,\,{\rm km}\,{\rm s}^{-1}$. These values are
consistent with what was found in massive NSCs (see Figure 2 of Stone et al.,
2017). The average stellar density and velocity dispersion measured inside a
sphere of $R_{\mathrm{s}}=0.5\,\,{\rm pc}\simeq\Delta x_{\rm min}$ from the
MBH become $\rho_{0}$ and $\sigma_{0}$, respectively, that are used to
estimate the TDR (see Section 2.3).

Figure 4: The radial profiles of the stellar density (top panel) and the
enclosed mass (bottom panel) from the NSC’s gravitational center in the NoTDA2
run at $t=0$, 5, 10, 40, and 100 Myr. This initial relaxation by $t\sim 5$ Myr
is due to our idealized setup following an analytic fit. After the initial
relaxation, the profiles do not change by a significant amount. For more
information, see Section 4.1. Figure 5: The gas accretion rate (GAR
$=\dot{M}_{\mathrm{BH,\,GA}}$; blue lines) and the tidal disruption accretion
rate (TDAR $=\dot{M}_{\mathrm{BH,\,TDA}}$; red lines) onto the MBH in the “Set
A” simulations (see Table 2). TDAR is found to be comparable to GAR in the
first $\lesssim 100$ Myr of evolution in all of the runs. In the early stage
of the MBH evolution, the TDAR grows with $M_{\mathrm{BH}}$, but it begins to
saturate or even decline after $\sim$100 Myr. While the GAR is limited at all
times by the Eddington rate (marked by black dashed lines) in all simulations,
the TDAR may occasionally exceed the Eddington rate (bottom panel with
$\gamma=2.0$). For more information, see Section 4.2. Figure 6: The MBH growth
histories in the “Set A” simulations (see Table 2). The black dashed line
indicates a model MBH growing at the Eddington rate. Top: the MBH without TDA
(the NoTDA1 run) grows only to $\sim 2\times\,10^{5}\,{\rm M}_{\odot}$,
whereas the MBHs with TDA occasionally grow faster than the Eddington rate.
Middle: the ratio of TDAR to the total BH accretion rate,
$\dot{M}_{\mathrm{BH,\,TDA}}/\dot{M}_{\mathrm{BH}}$. The relative contribution
of TDA towards the MBH’s growth peaks at $\sim$50% around $t=50-100$ Myr.
Bottom: the cumulative gas mass consumed by the MBH. For more information, see
Section 4.2.

###  4.2 TDA’s Impact on the MBH’s Growth: Dependence on the Power-law Index
($\gamma$) of the NSC

We now investigate the impact of our TDA model on the evolution of the MBH. In
particular, by controlling the two key input parameters $\gamma$ and
$r_{\mathrm{c}}$, we monitor how the MBH grows during the 200 Myr of
evolution.

Figure 7: Same as Figure 5 but for the “Set B” simulations (see Table 2). The
TDAR reaches its peak earlier when a smaller $r_{\mathrm{c}}$ is assumed. For
more information, see Section 4.3. Figure 8: Same as Figure 6 but for the “Set
B” simulations (see Table 2). The TDAR may (briefly) dominate over the GAR
when a smaller $r_{\mathrm{c}}$ is assumed. For more information, see Section
4.3.

First, the runs in the “Set A” suite (see Table 2) explore how the TDAR is
affected by $\gamma$, with a MBH seed of mass $M_{\mathrm{BH,\,init}}=8\times
10^{3}\,\mathrm{M}_{\odot}$. We are specifically interested in how the TDAR
compares with the GAR. Figure 5 illustrates the GARs (blue) and TDARs (red)
measured in $t=0-200$ Myr for the “Set A” simulations. The black dashed line
refers to the Eddington rate, $\dot{M}_{\mathrm{BH,\,Edd}}$ in Eq.(20),
corresponding to the mass of the MBH at that moment. What is the most worth
noting is that the TDAR is comparable to GAR in the first $\lesssim 100$ Myr
of evolution in all of the runs we have tested. In the early stage of the MBH
evolution, the TDAR grows with $M_{\mathrm{BH}}$, but it begins to saturate or
even decline after $\sim$100 Myr, which is a sign that our TDAR prescription
is working as intended (see Figure 1). The TDAR becomes higher with a larger
$\gamma$, as is predicted from the middle panel of Figure 1. Whereas the GAR
is limited at all times by the Eddington rate, the TDAR may exceed the
Eddington rate, e.g., from 60 to 100 Myr when with $\gamma=2.0$ (bottom panel
of Figure 5). This implies that TDA can dominate the mass supply to the MBH at
a certain epoch ($M_{\rm BH}\lesssim 10^{5}\,\mathrm{M}_{\odot}$).

In Figure 6 we further compare the MBH growth histories in the “Set A”
simulation suite. By $200\,\mathrm{Myr}$, the MBH growing only via GA reaches
$\sim 2\times 10^{5}\,{\rm M}_{\odot}$ in mass (the NoTDA1 run; top panel). In
contrast, the MBHs in all of the tested runs with TDA grow faster than the
Eddington rate (denoted by the black dashed line). These MBHs reach $\gtrsim
10^{6}\,{\rm M}_{\odot}$ by 200 Myr, about an order of magnitude larger in
mass than in the NoTDA1 run. It is worth noting that even the MBH in the
gamma1.5 run — which predicts the least TDAR among “Set A” (see the middle
panel of Figure 1) — still grows at a rate slightly higher than the Eddington
rate. In the middle panel of Figure 6, the ratio of TDAR to the total BH
accretion rate, $\dot{M}_{\mathrm{BH,\,TDA}}/\dot{M}_{\mathrm{BH}}$, is shown.
We see that the relative contribution of TDA towards the MBH’s accretion
reaches $\sim$50% during $50-100$ Myr ($M_{\rm BH}\lesssim
10^{5}\,\mathrm{M}_{\odot}$). Indeed, for $\gamma=1.75$ and $2.0$, the TDARs
exceed the GARs in this period. However, the ratio
$\dot{M}_{\mathrm{BH,\,TDA}}/\dot{M}_{\mathrm{BH}}$ tends to decline after
$100\,\mathrm{Myr}$. In the gamma1.5 and gamma1.75 runs, the ratio decreases
to $\lesssim 0.1$. In gamma2.0, the ratio stays above 0.25 after 150 Myr, but
it is mainly attributed to the stalled GA (see the bottom panel in Figure 5
and Section 4.4). Lastly, the bottom panel of Figure 6 depicts the cumulative
gas consumption by the MBH. Although the GA and TDA models operate
independently, TDA may still affect GAR indirectly — that is, when the MBH
growth is boosted by TDA, the GAR in Eq.(20) is also enhanced because of the
larger MBH mass.

###  4.3 TDA’s Impact on the MBH’s Growth: Dependence on the Core Size
($r_{\mathrm{c}}$) of the NSC

We move to the results of the “Set B” simulations (see Table 2; with a MBH
seed of $M_{\mathrm{BH,\,init}}=10^{3}\,\mathrm{M}_{\odot}$) to explore how
the TDAR is affected by the assumed NSC core size $r_{\mathrm{c}}$. In Figure
7, one can see that the TDAR correlates inversely with $r_{\mathrm{c}}$, as
our TDR prescription predicts in the right panel of Figure 1. We also find
that the peak of TDAR appears earlier when a smaller $r_{\mathrm{c}}$ is
assumed, again as expected from Figure 1. These behaviors are because a
smaller $r_{\mathrm{c}}$ increases the TDR (see Eq.(9)) and at the same time
reduces the transition mass $M_{\mathrm{c}}$ (see Eq.(12) and Section 2.3).

Figure 8 is the same as Figure 6 but for the “Set B” suite. As in Figure 6,
the MBH growth is significantly enhanced by TDA, especially in its early
evolutionary stage. The relative contribution of TDA towards the MBH’s
accretion reaches $\gtrsim$ 50% when $M_{\rm BH}\lesssim
10^{5}\,\mathrm{M}_{\odot}$. The TDAR is greatly affected by the assumed value
of $r_{\mathrm{c}}$. In the rc0.03 run, the MBH grows to $\sim 10^{7}\,{\rm
M}_{\odot}$, nearly 3 orders of magnitude greater than in our control run,
NoTDA2. By contrast, TDA makes little difference in the rc0.3 run, and the
MBHs in the rc0.3 run and the NoTDA2 run reaches similar masses at $t=200$
Myr. This observation emphasizes the importance of the assumed characteristics
of the NSC in our TDA prescription, such as $r_{\mathrm{c}}$ (for discussion
on our $r_{\mathrm{c}}$ choices, see Appendix B).

To summarize our findings, we discover that TDA can significantly boost the
growth of the seed MBH in most of the runs we tested with reasonable parameter
choices for the NSC’s structure, $\gamma$ and $r_{\mathrm{c}}$. It is
especially true in the early phase of the MBH’s growth up to $M_{\rm
BH}\lesssim 10^{5}\,\mathrm{M}_{\odot}$. During this phase, the
$\dot{M}_{\mathrm{BH,\,TDA}}$ grows with $M_{\mathrm{BH}}$. After
$\dot{M}_{\mathrm{BH,\,TDA}}$ reaches its peak, however, it begins to decline
as $M_{\mathrm{BH}}$ grows, so does the contribution of TDA towards the MBH’s
growth (for all the run with TDA in Table 2 except rc0.1 and rc0.3). In this
later phase, GA becomes the main mass supplier to the MBH, in agreement with
the conventional MBH growth model in a galaxy-scale numerical experiment.

###  4.4 TDA’s Impact on the MBH Host Galaxy: Star Formation, Morphology

The TDA model we have tested describes an interaction between the MBH particle
and its neighboring star particles. Yet, TDA’s influence reaches beyond this
simplistic interaction. In Section 4.3, we discussed that TDA may indirectly
affect the GAR by boosting the MBH’s mass. In this section, we consider TDA’s
impact on other aspects of the MBH host galaxy, such as its star formation
history and morphology.

Figure 9: The star formation histories (SFHs) of the “Set B” simulations (see
Table 2). Top: star formation rates (SFRs). Bottom: Total masses of newly-
formed stars in the simulation. In the rc0.03 and rc0.05 run, star formation
halts after $t\sim 100\,\mathrm{Myr}$. In these two runs, the total masses of
the newly-formed stars by $t\sim 200\,\mathrm{Myr}$ are only half of those in
the other runs.

We first examine the star formation history (SFH) to study this effect. Figure
9 shows the star formation rates (SFRs) and the cumulative masses of newly-
formed star particles as functions of time for the “Set B” suite. In all the
runs, stars form at a rate of $\gtrsim
0.03\,\mathrm{M}_{\odot}\mathrm{yr}^{-1}$ before $t\sim 50\,\mathrm{Myr}$.
However, in the rc0.03 and rc0.05 run in which TDA is more active than others,
star formation is quenched after $\sim$100 Myr. As a result, the cumulative
masses of the newly-formed stars by $t\sim 200\,\mathrm{Myr}$ in the two runs
are $\sim 7\times 10^{6}\,{\rm M}_{\odot}$, while that of the other runs is
$\gtrsim 1.3\times 10^{6}\,{\rm M}_{\odot}$ (rc0.1, rc0.3, NoTDA2). This mass
gap is approximately the same as the difference in the cumulative gas mass
consumed by the MBH during the same period (see the bottom panel of Figure 8).
It implies that the suppressed star formation in the rc0.03 and rc0.05 run is
likely due to the increased gas consumption by the MBH with efficient TDA.

Figure 10: The projected gas densities of the “Set B” simulations (see Table
2) centered on the location of the MBH at T = 150 Myr. In the NoTDA2 run, the
concentrated gas around the MBH is noticeable. In contrast, once we include
TDA, the gas density around the black hole decreases. In the rc0.03 run, for
example, one can see a clear cavity (of size $\sim$ a few pc) in the vicinity
of the MBH.

When TDA is most active, the intensified gas consumption by the MBH may also
reduce the gas density near the MBH. Sometimes it creates a visible cavity in
gas in the vicinity of the MBH. Figure 10 shows the projected gas densities of
the “Set B” simulations at $t$ = 150 Myr. When compared with the NoTDA2 run,
cavities in gas density around the MBH (of sizes $\sim$ a few pc) are
pronounced in all other runs with TDA. One may also see a trend that the
cavity size increases with decreasing $r_{\mathrm{c}}$. This morphological
change is because the massive MBH — growing expeditiously with both GA and TDA
— could consume its neighboring gas faster than being replenished by gas
inflows. One could argue that this gas cavity has stalled the GAR in some
epochs; e.g., see the bottom panel of Figure 5 for the gamma2.0 run, and the
2nd/3rd panel of Figure 7 for the rc0.03/rc0.05 run. In these runs, the GAR is
significantly lower than the Eddington rate after $\sim 150\,\mathrm{Myr}$ due
to the lack of gas supply around the MBH.

##  5 Discussion

###  5.1 Comparison With Observations and Future Observational Prospects

In our simulations, we find that the TDR often reaches
$10^{-2}-10^{-3}\,\mathrm{yr}^{-1}$. In contrast, the estimates based on the
observed TDE samples give $\sim 10^{-5}\,{\rm yr}^{-1}\,\mathrm{galaxy}^{-1}$
(e.g., Donley et al., 2002; Esquej et al., 2008; Gezari et al., 2009; van
Velzen & Farrar, 2014; Holoien et al., 2016). We however note that the actual
TDR would likely be higher than $\sim 10^{-5}\,{\rm
yr}^{-1}\,\mathrm{galaxy}^{-1}$, given that the current survey selection
technique may be missing some population of TDEs and that there may be an
error in the assumed BH mass function in the lower-mass galaxies (Gezari,
2021). It should also be noted that the theorists have predicted the rates
$\sim 10^{-4}-10^{-3}\,{\rm yr}^{-1}\,\mathrm{galaxy}^{-1}$ in the local
universe, a value higher than the aforementioned observed values (e.g.,
Magorrian & Tremaine, 1999; Wang & Merritt, 2004, see also van Velzen et al.
2018). Furthermore, the TDR can be greater at high redshift than in the local
universe (Kochanek, 2016).

The growth of a MBH by TDA may be verified by the observations in the near
future. Recently, Baldassare et al. (2022) presented an observational test of
Stone et al. (2017)’s work. They found that the NSCs with velocity dispersions
above $40\,\,\mathrm{km}\,\mathrm{s}^{-1}$ (a threshold suggested by Miller &
Davies 2012) are twice as likely to contain a MBH. In addition, measurements
of the spin distribution of MBHs may determine the contribution of TDA on the
growth of a MBH (Zhang et al., 2019). TDE flares at $z\gtrsim 3$ could also be
observed (Padmanabhan & Loeb, 2021).

###  5.2 Limitations of Our Model and Future Work

Our model includes several simplifying assumptions to estimate the TDAR. Here
we discuss how lifting these assumption would change our findings, and present
future work.

  1. 1.

Mass function of stars in the NSC: We have assumed the masses of stars in the
NSC are identically $m_{\star}=0.7\,{\rm M}_{\odot}$ (see Sections 2.2 and 2.4
and Eq.(14)). In reality, a NSC consists of stars with a wide range of masses.
In such a stellar system, the equipartition principle entails that more
massive stars migrate inward. This implies that the average stellar mass at
the center of the NSC can be higher than $m_{\star}=0.7\,{\rm M}_{\odot}$.
Then, the relaxation time, $t_{\mathrm{\,R}}$ in Eq.(11), at the center of the
NSC would become shorter than what our model expects. Correspondingly, the TDR
would be higher.111111The change in $m_{\star}$ may also entail a steeper
power-law index, and thus a change in the TDAR (Bahcall & Wolf, 1977).

  2. 2.

Fraction of the tidal debris that accretes to the MBH: We have assumed that
after TDEs the fraction $f$ of stellar debris that eventually accretes to the
MBH is 1.0 (see Section 2.4 and Eq.(14)). However, it is likely that not all
the debris falls onto the MBH. In typical TDEs, approximately a half of the
disrupted mass is bound to the compact object (Rees, 1988; Evans & Kochanek,
1989). More materials may become unbound if an energetic shock is produced by
the infalling gas, making the actual accretion efficiency $f$ lower than unity
(e.g., Strubbe & Quataert, 2009).121212Another consideration may lower the
fraction $f$. Studies on super-critical accretion show that after a TDE, the
bound material eventually returns to the pericenter at a super-Eddington rate
(e.g., Evans & Kochanek, 1989; Strubbe & Quataert, 2009) and then feeds the BH
(e.g. Ohsuga et al., 2005). When super-Eddington accretion is allowed, the
typical timescale for this process is $\sim$days, and most of the bound
material can be consumed by the BH before the next TDE. However, if the
accretion is limited by the Eddington rate, the next TDE will occur before all
the bound debris returns to the BH. In this case, the actual accretion
efficiency could be even lower.

  3. 3.

MBH binaries: Our TDR estimates are based on the assumption that the NSC
follows a power-law density profile with an index $\gamma$ (see Sections 2.2).
In reality, the assumption may not be valid. For example, a second MBH
spiraling into the NSC’s center may perturb the stellar density significantly
(e.g., Merritt & Szell, 2006). In such a case, a different model would be
needed to compute TDRs of MBH binaries (e.g., Li et al., 2017, 2019). During
the close encounter of two MBHs the $\dot{M}_{\mathrm{BH,\,TDA}}$ may
significantly increase, albeit temporarily, for two reasons. First, strong
perturbation from the MBH companion enhances the loss cone feeding. Second,
more stars populate the loss cone via 3-body interaction in the triaxial
stellar distribution during this phase. These processes cannot be properly
captured by our present model. Therefore, more sophisticated TDA model will be
required in the future.

  4. 4.

MBH-BH mergers: Our work focuses on the MBH growth by tidal disruption of
stars. But, stellar mass BHs may also contribute to the (birth and) growth of
the MBH. Antonini et al. (2019) pointed out that BHs with masses $\sim
100\,\,{\rm M}_{\odot}$ can grow in a star cluster with high density ($\gtrsim
10^{5}\,{\rm M}_{\odot}\,{\rm pc}^{-3}$) and high escape speed ($\gtrsim
300\,{\rm km}\,{\rm s}^{-1}$). Hong et al. (2018, 2020) showed that IMBHs with
masses $\gtrsim 10^{4}\,{\rm M}_{\odot}$ can grow by repeated mergers of
stellar BHs within $12\,\mathrm{Gyr}$. However, this channel is likely to be
suppressed as MBH grows, as stellar mass BHs are removed by mergers, or as
they run away from the galactic center (Miller & Davies, 2012; Hong et al.,
2020).

  5. 5.

Tidal captures of stars: A tidal impulse during a close encounter between a
compact object and a star may help to form a binary system (Fabian et al.,
1975). This process is called a tidal capture (TC). Studies have suggested
that a stellar mass black hole in a dense stellar cluster may grow to an IMBH
by successively capturing nearby stars (e.g., Miller & Davies, 2012; Stone et
al., 2017). However, the TC is likely to be deactivated for BHs with masses
$\gtrsim 10^{3}\,{\rm M}_{\odot}$ (Stone et al., 2017). Then the TDA problem
can be dealt with using the classical loss cone theory adopted in our model.
Since the initial masses of MBH seeds in our simulations are above
$10^{3}\,{\rm M}_{\odot}$, we choose not to consider the TCs this time.

  6. 6.

Resolved dynamics of the NSC: When we determine the stellar profile near the
MBH, we take a sub-resolution approach by assuming a power-law profile with an
index $\gamma$ irrespective of the profile found in the simulation (see
Sections 2.2 and 2.3, and footnote 9). On the other hand, Pfister et al.
(2021) selects $\gamma$ by computing stellar masses within two spheres of
radii $2\Delta x_{\mathrm{min}}$ and $4\Delta x_{\mathrm{min}}$. Yet,
extrapolating the stellar density from the scale of $\sim\Delta
x_{\mathrm{min}}$ down to $\sim r_{\mathrm{c}}$ may yield an inaccurate
$\gamma$ value when the MBH’s gravitational influence is not well resolved.
Thus, one may say either approach has its own limitations. In our future work,
we plan to resolve the NSC with e.g., a direct $N$-body routine that subcycles
in a hydrodynamic calculation.

  7. 7.

Gaseous dynamical friction: A massive perturber in a gaseous medium loses its
angular momentum and sinks into the center due to dynamical friction.
Dynamical friction may bring massive stars far outside into the NSC center.
Further, Boco et al. (2020, 2021) argued that dynamical friction can drive
multiple mergers, thus making MBHs to grow to $\sim 10^{4-6}\,{\rm M}_{\odot}$
in $\sim 10$ Myr. Higher-resolution hydrodynamic calculations may thus provide
more accurate BH growth rates.

The prediction of the MBH growth history based on our simulations should be
interpreted with caution due to the limitations above. Nonetheless, our sub-
resolution prescription for TDA in a high-resolution galaxy-scale hydrodynamic
simulation is the first step towards understanding the possible contribution
of TDA to the MBH growth.

##  6 Conclusion

We have introduced a model of MBHs’ growth via TDA and quantified its effects
on the MBH evolution using high-resolution AMR simulations. Despite its
possible contribution, the role of TDA towards the MBH growth has been
overlooked in most galaxy simulations. GA may dominate the accretion to the
MBHs with masses $\gtrsim 10^{6}\,\mathrm{M}_{\odot}$, but TDA may
significantly boost the growth of the MBH seeds with masses $\lesssim
10^{5}\,\mathrm{M}_{\odot}$. Yet, it is nontrivial to attain sufficient
resolution to describe the interactions between stars and the MBH in a galaxy-
scale simulation, due to the limited computational resources (Section 1). To
tackle the challenge, we have built a sub-resolution prescription adopting a
statistical approach with a few simplifying assumptions (Section 2). With the
new model that considers both GA and TDA in a galaxy-scale hydrodynamic
simulation (Section 3), we have acquired three main results (Section 4).

  1. 1.

TDA significantly enhances the MBH’s growth. In some simulations, a MBH seed
grows rapidly from $10^{3}\,{\rm M}_{\odot}$ to $\gtrsim
10^{6}\,\mathrm{M}_{\odot}$ in 200 Myrs (Figures 5 and 7). The growth rate is
more than an order of magnitude higher than in the run where a MBH grows only
via Eddington-limited GA. In general, the MBH grows faster if its host NSC is
assumed to have a higher power-law index ($\gamma$) and a smaller core size
($r_{\mathrm{c}}$).

  2. 2.

TDA mainly contributes to the early growth of MBH, from $10^{3-4}\,{\rm
M}_{\odot}$ to $\lesssim 10^{5}\,\mathrm{M}_{\odot}$. In most tested runs, the
relative contribution of TDA towards the MBH’s accretion reaches $\gtrsim$ 50%
when $M_{\rm BH}\lesssim 10^{5}\,\mathrm{M}_{\odot}$ (Figures 6 and 8). As the
MBH grows, TDA becomes sub-dominant while the later evolution is driven by GA.

  3. 3.

We also find that the star formation around the MBH is suppressed when TDA is
most active, because the massive MBH growing by both GA and TDA could rapidly
consume its neighboring gas. A cavity in gas (of size $\sim$ a few pc) is
sometimes visible near the MBH.

While our experiments are ideal to see the relative contribution of TDA versus
GA in a simplified setting (e.g., the NSC with a power-law profile), the TDAR
in reality may be different from our estimates. We also note that our sub-
resolution TDA model in simulations behaves in the exact way that it is
designed (as a function of $M_{\mathrm{BH}},\rho_{0},\sigma_{0},\gamma,$ and
$r_{\rm c}$), and the model is not intended to probe how the TDAR changes as a
result of the detailed stellar dynamics in the NSC resolved in the simulation.
And yet, our study is sufficient to show the possible contribution of TDA to
the rapid growth of MBHs. Our calculations demonstrate the need to consider
different channels of MBH accretion that may provide clues for the existence
of supermassive black holes at high redshifts. We also emphasize that this is
the first step towards exploring TDA for MBHs. In our future work, we aim not
only to advance our TDA model, but also to test it in various environments.
Further improvements to our model will make our prediction more sophisticated
and reliable.

We are deeply grateful to Hyung Mok Lee and Jong-Hak Woo for their invaluable
advice during this work. We thank the anonymous referee for providing us with
insightful comments that improved this article. We also thank Vivienne
Baldassare, Lumen Boco, Jongsuk Hong, Hamsa Padmanabhan, Vivienne Baldassare,
Hugo Pfister, and Xiaoxia Zhang for their detailed comments on our manuscript,
and Yongseok Jo, Seoyoung Kim, and Ahram Lee for their advice during the early
phase of this work. Ji-hoon Kim acknowledges support by Samsung Science and
Technology Foundation under Project Number SSTF-BA1802-04, and by the POSCO
Science Fellowship of POSCO TJ Park Foundation. His work was also supported by
the National Institute of Supercomputing and Network/Korea Institute of
Science and Technology Information with supercomputing resources including
technical support, grants KSC-2020-CRE-0219 and KSC-2021-CRE-0442. The
publicly available Enzo and yt codes used in this work are the products of
collaborative efforts by many independent scientists from numerous
institutions around the world. Their commitment to open science has helped
make this work possible.

##  Appendix A The Mass of A MBH Seed in Our Simulations

Recall that, to simplify the TDR estimates, we have assumed that the MBH stays
at the center of the NSC without “wandering” (see Section 2.2). This
assumption must hold true during the simulation to make our TDA model self-
consistent. This consideration yields a requirement for the initial MBH seed
mass in our runs, $M_{\mathrm{BH,\,init}}$ (see Section 3.1).

To show this, we start by imposing a criterion

$r_{\mathrm{wan}}<r_{\mathrm{\,infl}}$ (A1)

where $r_{\mathrm{wan}}$ is the wandering radius of the MBH that can be
estimated from

$\frac{v_{\mathrm{BH}}^{2}}{r_{\mathrm{wan}}}=\frac{GM_{\star}\left(r_{\mathrm{wan}}\right)}{r_{\mathrm{wan}}^{2}},$
(A2)

while the MBH’s radius of influence (Peebles, 1972) is written as

$r_{\mathrm{\,infl}}=\frac{GM_{\mathrm{BH}}}{\sigma_{\rm NSC}^{2}}\,.$ (A3)

In both Eqs.(A2) and (A3), a nearly isothermal stellar distribution is assumed
(i.e., $\gamma=2.0$ in Eq.(3)). Therefore, the stellar mass enclosed in
$r_{\mathrm{wan}}$ is
$M_{\star}\left(r_{\mathrm{wan}}\right)=(4\pi/3)\rho_{\mathrm{NSC}}r_{\mathrm{wan}}^{3}$.
Because of the isothermal sphere assumption (with
$\rho_{\mathrm{c}}=\rho_{\mathrm{NSC}}$), we can then write

$r_{\rm
wan}^{2}=\frac{v_{\mathrm{BH}}^{2}}{(4\pi/3)G\rho_{\mathrm{NSC}}}\,\,\,\,\,\,\,\,{\rm
and}\,\,\,\,\,\,\,\,r_{\rm
c}^{2}=\frac{\sigma_{\mathrm{NSC}}^{2}}{(4\pi/3)G\rho_{\mathrm{NSC}}}$ (A4)

to which the equipartition principle,
$M_{\mathrm{BH}}v_{\mathrm{BH}}^{2}\simeq
3m_{\star}\sigma_{\mathrm{NSC}}^{2}$, is applied to get

$r_{\mathrm{wan}}\simeq
r_{\mathrm{c}}\left(\frac{m_{\star}}{M_{\mathrm{BH}}}\right)^{{1\over 2}}.$
(A5)

Finally, plugging Eqs.(A3) and (A5) into Eq.(A1), we obtain

$\displaystyle
M_{\mathrm{BH,\,init}}>\left(\frac{r_{\mathrm{c}}^{2}\sigma_{\rm
NSC}^{4}m_{\star}}{G^{2}}\right)^{\frac{1}{3}}$ $\displaystyle\approx$
$\displaystyle\,\,3\times 10^{3}\,{\rm
M}_{\odot}\times\left(\frac{r_{\mathrm{c}}}{0.1\,{\rm
pc}}\right)^{\frac{2}{3}}\left(\frac{\sigma_{\rm NSC}}{100\,{\rm km}\,{\rm
s}^{-1}}\right)^{\frac{4}{3}}\left(\frac{m_{\star}}{0.7\,{\rm
M}_{\odot}}\right)^{\frac{1}{3}}$ (A6)

which gives a lower limit of $M_{\mathrm{BH,\,init}}$ that makes the MBH
stationary and our TDA model self-consistent (see Section 2.2). This motivates
our choices of $M_{\mathrm{BH,\,init}}$ in Section 3.1.

##  Appendix B The Range of the NSC’s Core Size in Our Model

Our TDA prescription computes $\dot{M}_{\mathrm{BH,\,TDA}}$ based on the user-
defined parameters, such as the power-law index $\gamma$ of the NSC’s density
profile, and the size of its stellar core, $r_{\mathrm{c}}$ (see Section 2.3).
Therefore, it is crucial to provide the model with well-constrained parameters
in simulations with TDA, in order to avoid an unrealistic MBH growth scenario.

In particular, given that a small change in $r_{\mathrm{c}}$ results in a
sizable change in the TDAR, it is essential to constrain $r_{\mathrm{c}}$ to
be adopted for the TDA model. Stone et al. (2017, see their Figure 3) made a
2D Gaussian fit between $\overline{\sigma}_{\rm NSC}$ and $r_{\mathrm{c}}$
using observational data (e.g., Böker et al., 2004; Côté et al., 2006;
Georgiev & Böker, 2014) — where $\overline{\sigma}_{\rm
NSC}\equiv\sqrt{GM_{\rm NSC}/(3r_{\rm\,NSC})}$ is the average 1D velocity
dispersion in the NSC — and found a relation
$r_{\mathrm{c}}\propto(\overline{\sigma}_{\rm NSC})^{-1.3}$. It means that the
NSC with greater $\overline{\sigma}_{\rm NSC}$ tends to have a smaller
$r_{\mathrm{c}}$, thus a denser NSC core and more active TDA. From their
Figure 3 one can also observe that $r_{\mathrm{c}}$ is approximately within
$[0.1,\,1]\,{\rm pc}$ for $\overline{\sigma}_{\rm NSC}\sim 100\,\,{\rm
km}\,{\rm s}^{-1}$. This motivated our chosen $r_{\mathrm{c}}$ range for
simulations, $[0.03,\,0.3]\,{\rm pc}$ (see Table 2). However, because the
$\overline{\sigma}_{\rm NSC}-r_{\mathrm{c}}$ relation has a large scatter, it
is difficult to use the relation to constrain $r_{\mathrm{c}}$ exactly.

## References

  * Aarseth et al. (1974) Aarseth, S. J., Henon, M., & Wielen, R. 1974, A&A, 37, 183
  * Antonini et al. (2019) Antonini, F., Gieles, M., & Gualandris, A. 2019, MNRAS, 486, 5008, doi: 10.1093/mnras/stz1149
  * Bahcall & Wolf (1976) Bahcall, J. N., & Wolf, R. A. 1976, ApJ, 209, 214, doi: 10.1086/154711
  * Bahcall & Wolf (1977) —. 1977, ApJ, 216, 883, doi: 10.1086/155534
  * Baldassare et al. (2022) Baldassare, V. F., Stone, N. C., Foord, A., Gallo, E., & Ostriker, J. P. 2022, ApJ, 929, 84, doi: 10.3847/1538-4357/ac5f51
  * Begelman & Rees (1978) Begelman, M. C., & Rees, M. J. 1978, MNRAS, 185, 847, doi: 10.1093/mnras/185.4.847
  * Binney & Tremaine (2011) Binney, J., & Tremaine, S. 2011, Galactic dynamics, Vol. 13 (Princeton university press)
  * Boco et al. (2020) Boco, L., Lapi, A., & Danese, L. 2020, ApJ, 891, 94, doi: 10.3847/1538-4357/ab7446
  * Boco et al. (2021) Boco, L., Lapi, A., Sicilia, A., et al. 2021, J. Cosmology Astropart. Phys, 2021, 035, doi: 10.1088/1475-7516/2021/10/035
  * Böker et al. (2004) Böker, T., Sarzi, M., McLaughlin, D. E., et al. 2004, AJ, 127, 105, doi: 10.1086/380231
  * Bondi (1952) Bondi, H. 1952, MNRAS, 112, 195, doi: 10.1093/mnras/112.2.195
  * Bondi & Hoyle (1944) Bondi, H., & Hoyle, F. 1944, MNRAS, 104, 273, doi: 10.1093/mnras/104.5.273
  * Booth & Schaye (2009) Booth, C. M., & Schaye, J. 2009, MNRAS, 398, 53, doi: 10.1111/j.1365-2966.2009.15043.x
  * Brummel-Smith et al. (2019) Brummel-Smith, C., Bryan, G., Butsky, I., et al. 2019, The Journal of Open Source Software, 4, 1636, doi: 10.21105/joss.01636
  * Bryan et al. (2014) Bryan, G. L., Norman, M. L., O’Shea, B. W., et al. 2014, ApJS, 211, 19, doi: 10.1088/0067-0049/211/2/19
  * Cen & Ostriker (1992) Cen, R., & Ostriker, J. P. 1992, ApJ, 399, L113, doi: 10.1086/186620
  * Choi et al. (2014) Choi, E., Naab, T., Ostriker, J. P., Johansson, P. H., & Moster, B. P. 2014, MNRAS, 442, 440, doi: 10.1093/mnras/stu874
  * Côté et al. (2006) Côté, P., Piatek, S., Ferrarese, L., et al. 2006, ApJS, 165, 57, doi: 10.1086/504042
  * Dehnen (1993) Dehnen, W. 1993, MNRAS, 265, 250, doi: 10.1093/mnras/265.1.250
  * Di Matteo et al. (2005) Di Matteo, T., Springel, V., & Hernquist, L. 2005, Nature, 433, 604, doi: 10.1038/nature03335
  * Donley et al. (2002) Donley, J. L., Brandt, W. N., Eracleous, M., & Boller, T. 2002, AJ, 124, 1308, doi: 10.1086/342280
  * Dubois et al. (2010) Dubois, Y., Devriendt, J., Slyz, A., & Teyssier, R. 2010, MNRAS, 409, 985, doi: 10.1111/j.1365-2966.2010.17338.x
  * Ebisuzaki et al. (2001) Ebisuzaki, T., Makino, J., Tsuru, T. G., et al. 2001, ApJ, 562, L19, doi: 10.1086/338118
  * Emsellem et al. (2015) Emsellem, E., Renaud, F., Bournaud, F., et al. 2015, MNRAS, 446, 2468, doi: 10.1093/mnras/stu2209
  * Esquej et al. (2008) Esquej, P., Saxton, R. D., Komossa, S., et al. 2008, A&A, 489, 543, doi: 10.1051/0004-6361:200810110
  * Evans & Kochanek (1989) Evans, C. R., & Kochanek, C. S. 1989, ApJ, 346, L13, doi: 10.1086/185567
  * Faber et al. (1997) Faber, S. M., Tremaine, S., Ajhar, E. A., et al. 1997, AJ, 114, 1771, doi: 10.1086/118606
  * Fabian et al. (1975) Fabian, A. C., Pringle, J. E., & Rees, M. J. 1975, MNRAS, 172, 15, doi: 10.1093/mnras/172.1.15P
  * Ferrarese et al. (2006a) Ferrarese, L., Côté, P., Dalla Bontà, E., et al. 2006a, ApJ, 644, L21, doi: 10.1086/505388
  * Ferrarese et al. (2006b) Ferrarese, L., Côté, P., Jordán, A., et al. 2006b, ApJS, 164, 334, doi: 10.1086/501350
  * Frank & Rees (1976) Frank, J., & Rees, M. J. 1976, MNRAS, 176, 633, doi: 10.1093/mnras/176.3.633
  * Freitag et al. (2006a) Freitag, M., Gürkan, M. A., & Rasio, F. A. 2006a, MNRAS, 368, 141, doi: 10.1111/j.1365-2966.2006.10096.x
  * Freitag et al. (2006b) Freitag, M., Rasio, F. A., & Baumgardt, H. 2006b, MNRAS, 368, 121, doi: 10.1111/j.1365-2966.2006.10095.x
  * Georgiev & Böker (2014) Georgiev, I. Y., & Böker, T. 2014, MNRAS, 441, 3570, doi: 10.1093/mnras/stu797
  * Gezari (2021) Gezari, S. 2021, ARA&A, 59, 21, doi: 10.1146/annurev-astro-111720-030029
  * Gezari et al. (2009) Gezari, S., Heckman, T., Cenko, S. B., et al. 2009, ApJ, 698, 1367, doi: 10.1088/0004-637X/698/2/1367
  * Giersz & Spurzem (1994) Giersz, M., & Spurzem, R. 1994, MNRAS, 269, 241, doi: 10.1093/mnras/269.2.241
  * Graham & Spitler (2009) Graham, A. W., & Spitler, L. R. 2009, MNRAS, 397, 2148, doi: 10.1111/j.1365-2966.2009.15118.x
  * Gürkan et al. (2004) Gürkan, M. A., Freitag, M., & Rasio, F. A. 2004, ApJ, 604, 632, doi: 10.1086/381968
  * Hernquist (1990) Hernquist, L. 1990, ApJ, 356, 359, doi: 10.1086/168845
  * Holoien et al. (2016) Holoien, T. W. S., Kochanek, C. S., Prieto, J. L., et al. 2016, MNRAS, 463, 3813, doi: 10.1093/mnras/stw2272
  * Hong et al. (2020) Hong, J., Askar, A., Giersz, M., Hypki, A., & Yoon, S.-J. 2020, MNRAS, 498, 4287, doi: 10.1093/mnras/staa2677
  * Hong et al. (2018) Hong, J., Vesperini, E., Askar, A., et al. 2018, MNRAS, 480, 5645, doi: 10.1093/mnras/sty2211
  * Hopkins & Quataert (2010) Hopkins, P. F., & Quataert, E. 2010, MNRAS, 407, 1529, doi: 10.1111/j.1365-2966.2010.17064.x
  * Hopkins & Quataert (2011) —. 2011, MNRAS, 415, 1027, doi: 10.1111/j.1365-2966.2011.18542.x
  * Hopkins et al. (2006) Hopkins, P. F., Somerville, R. S., Hernquist, L., et al. 2006, ApJ, 652, 864, doi: 10.1086/508503
  * Inayoshi et al. (2020) Inayoshi, K., Visbal, E., & Haiman, Z. 2020, ARA&A, 58, 27, doi: 10.1146/annurev-astro-120419-014455
  * Johansson et al. (2009) Johansson, P. H., Burkert, A., & Naab, T. 2009, ApJ, 707, L184, doi: 10.1088/0004-637X/707/2/L184
  * Kim et al. (2019) Kim, J.-h., Wise, J. H., Abel, T., et al. 2019, ApJ, 887, 120, doi: 10.3847/1538-4357/ab510b
  * Kim et al. (2011) Kim, J.-h., Wise, J. H., Alvarez, M. A., & Abel, T. 2011, ApJ, 738, 54, doi: 10.1088/0004-637X/738/1/54
  * Kim et al. (2014) Kim, J.-h., Abel, T., Agertz, O., et al. 2014, ApJS, 210, 14, doi: 10.1088/0067-0049/210/1/14
  * Kim et al. (2016) Kim, J.-h., Agertz, O., Teyssier, R., et al. 2016, ApJ, 833, 202, doi: 10.3847/1538-4357/833/2/202
  * Kochanek (2016) Kochanek, C. S. 2016, MNRAS, 461, 371, doi: 10.1093/mnras/stw1290
  * Kormendy & Ho (2013) Kormendy, J., & Ho, L. C. 2013, ARA&A, 51, 511, doi: 10.1146/annurev-astro-082708-101811
  * Lauer et al. (1998) Lauer, T. R., Faber, S. M., Ajhar, E. A., Grillmair, C. J., & Scowen, P. A. 1998, AJ, 116, 2263, doi: 10.1086/300617
  * Li et al. (2019) Li, S., Berczik, P., Chen, X., et al. 2019, ApJ, 883, 132, doi: 10.3847/1538-4357/ab3e4a
  * Li et al. (2017) Li, S., Liu, F. K., Berczik, P., & Spurzem, R. 2017, ApJ, 834, 195, doi: 10.3847/1538-4357/834/2/195
  * Lynden-Bell (1969) Lynden-Bell, D. 1969, Nature, 223, 690, doi: 10.1038/223690a0
  * Lynden-Bell (1978) —. 1978, Phys. Scr, 17, 185, doi: 10.1088/0031-8949/17/3/009
  * Lynden-Bell & Rees (1971) Lynden-Bell, D., & Rees, M. J. 1971, MNRAS, 152, 461, doi: 10.1093/mnras/152.4.461
  * Magorrian & Tremaine (1999) Magorrian, J., & Tremaine, S. 1999, MNRAS, 309, 447, doi: 10.1046/j.1365-8711.1999.02853.x
  * Merritt & Szell (2006) Merritt, D., & Szell, A. 2006, ApJ, 648, 890, doi: 10.1086/506010
  * Miller & Davies (2012) Miller, M. C., & Davies, M. B. 2012, ApJ, 755, 81, doi: 10.1088/0004-637X/755/1/81
  * Navarro et al. (1997) Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490, 493, doi: 10.1086/304888
  * Neumayer et al. (2020) Neumayer, N., Seth, A., & Böker, T. 2020, A&A Rev., 28, 4, doi: 10.1007/s00159-020-00125-0
  * Nguyen et al. (2018) Nguyen, D. D., Seth, A. C., Neumayer, N., et al. 2018, ApJ, 858, 118, doi: 10.3847/1538-4357/aabe28
  * Oh et al. (2020) Oh, B. K., Smith, B. D., Peacock, J. A., & Khochfar, S. 2020, MNRAS, 497, 5203, doi: 10.1093/mnras/staa2318
  * Ohsuga et al. (2005) Ohsuga, K., Mori, M., Nakamoto, T., & Mineshige, S. 2005, ApJ, 628, 368, doi: 10.1086/430728
  * Padmanabhan & Loeb (2021) Padmanabhan, H., & Loeb, A. 2021, A&A, 656, A47, doi: 10.1051/0004-6361/202040273
  * Peebles (1972) Peebles, P. J. E. 1972, General Relativity and Gravitation, 3, 63, doi: 10.1007/BF00755923
  * Perret (2016) Perret, V. 2016, DICE: Disk Initial Conditions Environment, Astrophysics Source Code Library, record ascl:1607.002.  http://ascl.net/1607.002
  * Pfister et al. (2021) Pfister, H., Dai, J. L., Volonteri, M., et al. 2021, MNRAS, 500, 3944, doi: 10.1093/mnras/staa3471
  * Portegies Zwart et al. (2004) Portegies Zwart, S. F., Baumgardt, H., Hut, P., Makino, J., & McMillan, S. L. W. 2004, Nature, 428, 724, doi: 10.1038/nature02448
  * Portegies Zwart & McMillan (2002) Portegies Zwart, S. F., & McMillan, S. L. W. 2002, ApJ, 576, 899, doi: 10.1086/341798
  * Quinlan & Shapiro (1987) Quinlan, G. D., & Shapiro, S. L. 1987, ApJ, 321, 199, doi: 10.1086/165624
  * Rees (1988) Rees, M. J. 1988, Nature, 333, 523, doi: 10.1038/333523a0
  * Roca-Fàbrega et al. (2020) Roca-Fàbrega, S., Kim, J.-h., Primack, J. R., et al. 2020, arXiv e-prints, arXiv:2001.04354.  https://arxiv.org/abs/2001.04354
  * Roca-Fàbrega et al. (2021) Roca-Fàbrega, S., Kim, J.-H., Hausammann, L., et al. 2021, ApJ, 917, 64, doi: 10.3847/1538-4357/ac088a
  * Salpeter (1955) Salpeter, E. E. 1955, ApJ, 121, 161, doi: 10.1086/145971
  * Salpeter (1964) —. 1964, ApJ, 140, 796, doi: 10.1086/147973
  * Schödel et al. (2018) Schödel, R., Gallego-Cano, E., Dong, H., et al. 2018, A&A, 609, A27, doi: 10.1051/0004-6361/201730452
  * Schödel et al. (2009) Schödel, R., Merritt, D., & Eckart, A. 2009, A&A, 502, 91, doi: 10.1051/0004-6361/200810922
  * Sesana (2012) Sesana, A. 2012, Advances in Astronomy, 2012, 805402, doi: 10.1155/2012/805402
  * Seth et al. (2008) Seth, A., Agüeros, M., Lee, D., & Basu-Zych, A. 2008, ApJ, 678, 116, doi: 10.1086/528955
  * Sijacki et al. (2009) Sijacki, D., Springel, V., & Haehnelt, M. G. 2009, MNRAS, 400, 100, doi: 10.1111/j.1365-2966.2009.15452.x
  * Spitzer & Harm (1958) Spitzer, Lyman, J., & Harm, R. 1958, ApJ, 127, 544, doi: 10.1086/146486
  * Springel et al. (2005a) Springel, V., Di Matteo, T., & Hernquist, L. 2005a, MNRAS, 361, 776, doi: 10.1111/j.1365-2966.2005.09238.x
  * Springel et al. (2005b) —. 2005b, ApJ, 620, L79, doi: 10.1086/428772
  * Stone et al. (2017) Stone, N. C., Küpper, A. H. W., & Ostriker, J. P. 2017, MNRAS, 467, 4180, doi: 10.1093/mnras/stx097
  * Stone & Ostriker (2015) Stone, N. C., & Ostriker, J. P. 2015, ApJ, 806, L28, doi: 10.1088/2041-8205/806/2/L28
  * Strubbe & Quataert (2009) Strubbe, L. E., & Quataert, E. 2009, MNRAS, 400, 2070, doi: 10.1111/j.1365-2966.2009.15599.x
  * Syer & Ulmer (1999) Syer, D., & Ulmer, A. 1999, MNRAS, 306, 35, doi: 10.1046/j.1365-8711.1999.02445.x
  * Takahashi (1995) Takahashi, K. 1995, PASJ, 47, 561.  https://arxiv.org/abs/astro-ph/9507040
  * Tremaine et al. (1994) Tremaine, S., Richstone, D. O., Byun, Y.-I., et al. 1994, AJ, 107, 634, doi: 10.1086/116883
  * van Velzen et al. (2018) van Velzen, S., Bower, G. C., & Metzger, B. D. 2018, in Astronomical Society of the Pacific Conference Series, Vol. 517, Science with a Next Generation Very Large Array, ed. E. Murphy, 737.  https://arxiv.org/abs/1810.06677
  * van Velzen & Farrar (2014) van Velzen, S., & Farrar, G. R. 2014, ApJ, 792, 53, doi: 10.1088/0004-637X/792/1/53
  * Volonteri (2010) Volonteri, M. 2010, A&A Rev., 18, 279, doi: 10.1007/s00159-010-0029-x
  * Wang & Merritt (2004) Wang, J., & Merritt, D. 2004, ApJ, 600, 149, doi: 10.1086/379767
  * Wang et al. (2015) Wang, L., Spurzem, R., Aarseth, S., et al. 2015, MNRAS, 450, 4070, doi: 10.1093/mnras/stv817
  * Zel’dovich (1964) Zel’dovich, Y. B. 1964, Soviet Physics Doklady, 9, 195
  * Zhang et al. (2019) Zhang, X., Lu, Y., & Liu, Z. 2019, ApJ, 877, 143, doi: 10.3847/1538-4357/ab1d48