# Tracing the Evolution of SMBHs and Stellar Objects in Galaxy Mergers: An
Multi-mass Direct N-body Model

Shuo Li National Astronomical Observatories and Key Laboratory of
Computational Astrophysics,
Chinese Academy of Sciences, 20A Datun Rd., Chaoyang District, Beijing 100012,
China  Department of Astronomy, School of Physics, Peking University,
Yiheyuan Lu 5, Haidian Qu, Beijing 100871, China  Shiyan Zhong Yunnan
Observatories, Chinese Academy of Sciences,
396 Yang-Fang-Wang, Guandu District, 650216, Kunming, Yunnan, China  Peter
Berczik Astronomisches Rechen-Institut, Zentrum für Astronomie, University of
Heidelberg,
Möchhofstrasse 12-14, Heidelberg 69120, Germany  Konkoly Observatory, Research
Centre for Astronomy and Earth Sciences, Eötvös Loránd Research Network
(ELKH), MTA Centre of Excellence, Konkoly Thege Miklós út 15-17, 1121
Budapest, Hungary  Main Astronomical Observatory, National Academy of Sciences
of Ukraine,
27 Akademika Zabolotnoho St., 03680 Kyiv, Ukraine  Rainer Spurzem National
Astronomical Observatories and Key Laboratory of Computational Astrophysics,
Chinese Academy of Sciences, 20A Datun Rd., Chaoyang District, Beijing 100012,
China  Astronomisches Rechen-Institut, Zentrum für Astronomie, University of
Heidelberg,
Möchhofstrasse 12-14, Heidelberg 69120, Germany  Kavli Institute for Astronomy
and Astrophysics, Peking University,
Yiheyuan Lu 5, Haidian Qu, Beijing 100871, China  Xian Chen Department of
Astronomy, School of Physics, Peking University,
Yiheyuan Lu 5, Haidian Qu, Beijing 100871, China  Kavli Institute for
Astronomy and Astrophysics, Peking University,
Yiheyuan Lu 5, Haidian Qu, Beijing 100871, China  F.K. Liu Department of
Astronomy, School of Physics, Peking University,
Yiheyuan Lu 5, Haidian Qu, Beijing 100871, China  Kavli Institute for
Astronomy and Astrophysics, Peking University,
Yiheyuan Lu 5, Haidian Qu, Beijing 100871, China

###### Abstract

By using direct N-body numerical simulations, we model the dynamical co-
evolution of two supermassive black holes (SMBHs) and the surrounding stars in
merging galaxies. In order to investigate how different stellar components
evolve during the merger, we generate evolved stellar distributions with an
initial mass function. Special schemes have also been developed to deal with
some rare but interesting events, such as tidal disruption of main sequence
stars, the plunge of low mass stars, white dwarfs, neutron stars and stellar
mass black holes, and the partial tidal disruption of red giants or asymptotic
giant branch stars. Our results indicate that the formation of a bound
supermassive black hole binary (SMBHB) will enhance the capture rates of
stellar objects by the SMBHs. Compared to the equal stellar mass model, the
multi-mass model tends to result in a higher average mass of disrupted stars.
Instead of being tidally disrupted by the SMBH, roughly half of the captured
main sequence stars will directly plunge into the SMBH because of their small
stellar radius. Giant stars, on the other hand, can be stripped of their
envelopes if they are close enough to the SMBH. Though most remnants of the
giant stars can survive after the disruption, a small fraction still could
plunge into the SMBH quickly or after many orbital periods. Our results also
indicate significant mass segregation of compact stars at the beginning of the
merger, and then this effect is destroyed as the two SMBHs form a bound
binary.

Galaxies: evolution — Galaxies: interactions — Galaxies: kinematics and
dynamics — Galaxies: nuclei — Methods: numerical

##  1 Introduction

Supermassive black hole binaries (SMBHBs) are predicted as the descendents of
the hierarchical galaxy formation model (Begelman et al., 1980; Volonteri et
al., 2003). Over the past few decades, more and more observational evidence
has indicated that most massive galaxies, if not all, have a supermassive
black hole (SMBH) hidden in the center (Kormendy & Ho, 2013). Since massive
galaxies could undergo several mergers in their evolutionary history, it is
natural to predict the existence of SMBHBs in merging or merged galaxies.
Besides, many investigations find that there is a close connection between the
SMBH and its host galaxy (Magorrian et al., 1998; Ferrarese & Merritt, 2000;
Gebhardt et al., 2000; Tremaine et al., 2002; Kormendy & Ho, 2013). Therefore
it is essential to find out how SMBHBs evolve in galaxy mergers.

In a merging galaxy, two SMBHs will first approach each other mainly through
dynamical friction. However, as two SMBHs get closer, this effect gets weaker
and weaker. When two SMBHs are close enough to form a bound binary system, the
dynamical friction is not efficient enough to drive the SMBHB coalesce within
Hubble time in most of the massive galaxies (Begelman et al., 1980). A hard
SMBHB can eject surrounding stars to transfer their orbital energy and angular
momentum, which may be efficient to drive two SMBHs coalesce quickly (Saslaw
et al., 1974; Mikkola & Valtonen, 1992; Quinlan, 1996). But it needs enough
close encounter stars to be ejected by the SMBHB, which may be not always the
case, because stars scattered into the vicinity of the SMBHB through two-body
relaxation may be inefficient in spherical stellar distributions (Begelman et
al., 1980; Yu, 2002; Milosavljević & Merritt, 2003; Berczik et al., 2005).
Fortunately, both gas dynamics and more realistic stellar dynamics other than
spherical two-body relaxation can avoid this problem (Gould & Rix, 2000;
Chatterjee et al., 2003; Merritt & Poon, 2004; Berczik et al., 2006; Preto et
al., 2011; Khan et al., 2011).

There are many observational evidence to confirm the above scenario, such as
dual active galactic nuclei (AGNs), jet reorientation in X-shaped radio
galaxies, double-peaked emission lines, quasi-periodic outbursts in some
blazars and so on (Liu & Wu, 2002; Komossa et al., 2003; Liu, 2004; Shen et
al., 2013; Komossa et al., 2020; Tang et al., 2021). Most of these are
indirect evidence of SMBHBs and they are in gaseous environments. In gas poor
environments, hard SMBHBs are very difficult to detect. Besides emitting
electromagnetic waves, a close SMBHB has strong gravitational wave (GW)
emission, which could be detected by the ongoing Pulsar Timing Arrays (PTAs)
and the planned space borne GW detectors such as the Laser Interferometer
Space Antenna (LISA)111https://www.elisascience.org/ , Taiji Program and
Tianqin Project222http://tianqin.sysu.edu.cn/en/ (Foster & Backer, 1990;
Verbiest et al., 2016; Mei et al., 2021; Luo et al., 2021). However, after
years of effort, though PTAs may already get some clue on the stochastic GW
background, the specific detection is still missing (Shannon et al., 2015;
Babak et al., 2016; Arzoumanian et al., 2016, 2020). There are also other
processes connected to SMBHs and SMBHBs, which could generate GW emission for
LISA or similar instruments. For instance, compact stellar objects can
inspiral to SMBHs with detectable GW emissions for LISA, which is the so-
called extreme mass ratio inspirals (EMRIs) (Amaro-Seoane, 2018, and
references therein). Similar events should also exist in the SMBHB system. But
it has not been well studied. Nevertheless, due to limited spatial resolutions
of GW detectors, even if a GW signal from the SMBHB could be detected, it is
still difficult to locate its host galaxy. The inconsistency between
theoretical expectation on SMBHBs and GW observations may be due to many
processes that have not been clarified. It is important to design observations
that can be used to best constrain the evolution of SMBHBs.

A dormant SMBH can be temporarily illuminated by tidal disruption events
(TDEs). If a star closely approaches the SMBH by less than a critical
distance, it will be torn into debris by the SMBH, which may prompt temporary
flares with periods from days to years (Hills, 1975; Rees, 1988; Evans &
Kochanek, 1989; Guillochon & Ramirez-Ruiz, 2013). We use for the critical
distance tidal radius

$r_{\rm t}\backsimeq\mu\rm r_{*}(M_{\rm BH}/\rm m_{*})^{1/3},$ (1)

where $\mu$ is a dimensionless parameter of order unity which reflects the
stellar structure, $\rm r_{*}$, $\rm m_{*}$ and $M_{\rm BH}$ are the stellar
radius, the stellar mass and the mass of black hole (BH), respectively. Since
the first identified event in the 1990s, dozens of TDEs have been reported,
with emission range from $\gamma$-ray to radio bands (Komossa & Bade, 1999;
Komossa, 2015; Gezari, 2021, and references therein). Similar events can also
happen in SMBHB systems, but light curves may be different from normal single
SMBH TDEs. Due to the perturbation of the companion SMBH, a TDE in SMBHB
system could have repeated gaps or be truncated, which has been investigated
theoretically and numerically (Liu et al., 2009; Ricarte et al., 2016;
Coughlin et al., 2017; Vigneron et al., 2018). With the expansion of the TDE
sample in these years, a few SMBHB TDE candidates also have been founded in
observation (Liu et al., 2014; Shu et al., 2020; Huang et al., 2021). With
more and more optical and X-ray transient surveys, such as the Large Synoptic
Survey Telescope (LSST)333https://www.lsst.org, the All-Sky Automated Survey
for Supernovae (ASAS-SN)444http://www.astronomy.ohio-
state.edu/asassn/index.shtml, and the Einstein Probe
(EP)555http://ep.bao.ac.cn/, join the game, the sample of the SMBHB TDE could
significantly increase in the near future.

There are some attempts to investigate the tidal disruption rate (TDR) of
SMBHBs. Theoretical analyses and numerical scattering experiments indicate
that, due to the perturbation of the companion SMBH and the triaxial stellar
distribution, TDR in galaxy merger remnants could be enhanced from a few times
to a few orders of magnitudes (Ivanov et al., 2005; Chen et al., 2009; Wegg &
Nate Bode, 2011; Liu & Chen, 2013). However, since the interaction of two
SMBHs and surrounding stars is chaotic during the formation of the SMBHB,
these results are limited and can not fully reveal the underlying physical
processes. In order to investigate the dynamical co-evolution of SMBHBs and
stars in galaxy mergers with more depth, we used GPU accelerated direct
N$-body~{}$ simulations to analyze a series of models on the TDR evolution of
SMBHBs in both major mergers Li et al. (2017, hereafter Paper I) and minor
mergers Li et al. (2019, hereafter Paper II). Both equal and unequal mass
models indicate significantly enhanced TDRs during two SMBHs forming bound
binary systems, which can be considered as a possible explanation for the high
detection rates of TDEs preferred in E+A galaxies (Arcavi et al., 2014).
However, for simplicity, it was assumed in both papers that all stars are
solar type stars with the same mass. In reality, for example, giant stars will
be partially tidally disrupted by the SMBH at relatively large distances.
Compact objects (neutron stars, black holes), on the other hand, could
directly plunge into the SMBH without disruption. Even main sequence stars
could have different fates after considering different masses. According to
Eq. 1, massive main sequence stars with larger stellar radius correspond to
larger tidal radius. Low mass main sequence stars usually correspond to
smaller tidal radius. Massive stars with large envelopes may as well have only
partial mass loss from TDE (see e.g.Zhong et al. (2022)). Finally, instead of
a tidal disruption with significant flare, many main sequence stars near the
low mass end actually will be entirely swallowed by the SMBH, because their
tidal radii are very close to the Schwarzschild radius of the SMBH.

In this work we improve our model by using a realistic distribution of stellar
properties (mass and radius) instead of equal mass stars as in previous papers
such as e.g. Li et al. (2019). In a dry major merger case, the typical
evolution time is $\sim 1{\rm\,Gyr}$ (Colpi, 2014). This period of time is
long enough for a stellar system to evolve from the zero-age main sequence to
different stellar components. For this reason, our stellar population is
initialized with an age of $\sim 1{\rm\,Gyr}$ to account for the lifetime of
the galaxies before their merger. We developed a special scheme to deal with
TDEs of different types of stars, which includes variation of the tidal
disruption radius as a function of stellar parameters, partial tidal
disruption of giant stars, and direct plunges of compact objects (see for
details next section and Fig. 1). With these modifications we are able to
trace the co-evolution of all kinds of stars and SMBHs in galaxy mergers. This
paper is organized as follows. We introduce our simulation models in Section
2. Our simulation results on different types of stars are presented in Section
3. In Section 4 we discuss how to extrapolate our results to more realistic
systems, and some observational implications are also provided. A short
summary is allocated to Section 5.

##  2 The Direct N-body Model With More Realistic Stellar Objects

Our numerical model base on a direct N-body code $\varphi$ -Grape/$\varphi$
-GPU, which is an accurate GPU accelerated code with fourth-order Hermite
integrator and efficient parallel scheme (Makino & Aarseth, 1992; Berczik et
al., 2005; Harfst et al., 2007). $\varphi$ -Grape/$\varphi$ -GPU has been
proved to be an efficient tool on investigating the dynamical co-evolution of
SMBH and stars in both single galaxy and galaxy mergers (Berczik et al., 2006;
Gualandris & Merritt, 2008; Khan et al., 2011; Preto et al., 2011). It can be
also adapted to investigate the tidal disruption evolution, with a simplified
tidal disruption scheme included (Zhong et al., 2014; Li et al., 2017;
Panamarev et al., 2018). Recently, Khan et al. (2018) have tried to study the
coalescence time of supermassive black holes in galaxy mergers with post-
Newtonian (PN) terms and stellar mass function included. In this work, we try
to introduce a more realistic model, which can self-consistently study how
different type stellar objects, such as main sequence(MS) star, red giant
(RG), asymptotic giant branch(AGB), white dwarf (WD), neutron star (NS),
stellar mass black hole (BH), interact with SMBHs in galaxy mergers.

###  2.1 Initial Mass Function and Stellar Evolution

Compared with previous equal stellar mass models, a more realistic multi-mass
model with initial mass function is adopted. We assume that the stellar mass
of a star cluster ranges from $0.1\,{\rm M_{\odot}}$ to $100\,{\rm
M_{\odot}}$, and the initial mass function follows a multiple-part power-law
as Kroupa (2001) suggested. This model corresponds to an average initial mass
$\overline{m}_{*}\sim 0.6\,{\rm M_{\odot}}$.

For simplicity, we are focusing on dry mergers with gas poor environments,
which are usually dominated by early type galaxies, and the starbursts induced
by mergers are not significant. It is reasonable to assume that all the stars
have evolved for a relatively long time as the initial condition. We evolve
all the stars for $1~{}{\rm\,Gyr}$ by using the stellar evolution package SSE
(Hurley et al., 2000). After the evolution, our model with largest particle
number $N=10^{6}$ has $4484$ giant branch stars (includes RGs, core He burning
stars, early AGBs and thermally pulsing AGBs), $29811$ WDs(includes C/O WDs
and O/Ne WDs), $5246$ NSs and 1884 BHs. And the average mass decreased to
$\overline{m}_{*}\sim 0.4\,{\rm M_{\odot}}$ due to the mass loss during the
stellar evolution. The same scheme as Panamarev et al. (2019) adopted is
involved to model the natal kick during the formation of NSs and BHs. The kick
amplitude of NSs is represented by a Maxwellian distribution with 1D velocity
dispersion $\sigma=265{\rm\,km\,s^{-1}}$ (Hobbs et al., 2005). For BHs, both
the mass and kick velocity sensitively depend on the ”fallback” of debris,
those materials failed to escape during the explosion of the progenitors
(Colgate, 1971; Chevalier, 1989; Zhang et al., 2008). We take this effect into
account by following the scheme suggested by Belczynski et al. (2002).

As a result, a group of stars with different stellar components including NSs
and BHs with initial kick velocities have been generated. In order to model a
dense star cluster, the next step is to assign position and velocity to every
star according to a proper stellar mass distribution.

###  2.2 Dense Star Cluster Model

The dynamical parameters of the initial condition are similar to those in
Paper I. The stellar distribution of the dense star cluster around the SMBH is
represented by a Dehnen model (Dehnen, 1993).

$\rho(r)=\frac{3-\gamma}{4\pi}\frac{Ma_{\rm D}}{r^{\gamma}(r+a_{\rm
D})^{4-\gamma}},$ (2)

where $a_{\rm D}$, $M$ and $\gamma$ denote the scaling radius, the total mass
of the galaxy/nucleus and the density profile index, respectively. Here we
adopt the units ${\rm G}=M=a_{\rm D}=1$, and assume $\gamma=1$ in the
following discussion for simplicity, because models with more steep central
density profiles are very time consuming in the integration. The influence of
$\gamma$ has been carefully discussed in Paper I, with equal stellar mass
models. The general results should be still informative here. According to
results in Paper I, steep density profiles usually correspond to significantly
higher TDRs, and the boosted TDRs in phase II are common for all models with
different $\gamma$. However, the magnitude of the enhancement of the averaged
TDRs in phase II only weakly depends on $\gamma$. Detailed discussions can be
found in Paper I. The relation between numerical and physical quantities can
be derived as

$\displaystyle[\rm{T}]$ $\displaystyle=$ $\displaystyle\left(\frac{{\rm
G}M}{a_{\rm D}^{3}}\right)^{-1/2}$ (3) $\displaystyle=$ $\displaystyle
1.491\times 10^{7}(2^{\frac{1}{3-\gamma}}-1)^{3/2}\left(\frac{M}{10^{9}\,{\rm
M_{\odot}}}\right)^{-1/2}\left(\frac{r_{1/2}}{1{\rm\,kpc}}\right)^{3/2}{\rm\,yr},$
$\displaystyle[\rm{V}]$ $\displaystyle=$ $\displaystyle\left(\frac{{\rm
G}M}{a_{\rm D}}\right)^{1/2}$ (4) $\displaystyle=$ $\displaystyle
65.58\times(2^{\frac{1}{3-\gamma}}-1)^{-1/2}\left(\frac{M}{10^{9}\,{\rm
M_{\odot}}}\right)^{1/2}\left(\frac{r_{1/2}}{1{\rm\,kpc}}\right)^{-1/2}{\rm\,km\,s^{-1}},$
$\displaystyle[\rm{R}]$ $\displaystyle=$ $\displaystyle a_{\rm
D}=(2^{\frac{1}{3-\gamma}}-1)\left(\frac{r_{1/2}}{1{\rm\,kpc}}\right){\rm\,kpc},$
(5) $\displaystyle[\dot{\rm{M}}]$ $\displaystyle=$ $\displaystyle M/[\rm T]$
(6) $\displaystyle=$ $\displaystyle
67.07\times(2^{\frac{1}{3-\gamma}}-1)^{-3/2}\left(\frac{M}{10^{9}\,{\rm
M_{\odot}}}\right)^{3/2}\left(\frac{r_{1/2}}{1{\rm\,kpc}}\right)^{-3/2}\,{\rm
M_{\odot}}/{\rm\,yr}.$

The SMBH is represented by a heavy particle with mass $M_{\rm BH}=0.01$ at the
center. The Dehnen model above does not considered the influence of SMBH and
multiple stellar mass distribution. Similar to Paper I, to relax the SMBH with
surrounding stars, we integrate the entire system for dozens of N-body unit
time before the simulation, which is roughly the two-body relaxation timescale
at the influence radius of the SMBH in the model. After this relaxation
procedure, the mass segregation of heavy stellar components in the central
region is significant. It has been confirmed by the result in the left panel
of Fig. 4 Based on this template model, two identical galaxies/nuclei are set
in a parabolic orbit with initial separation $d\sim 20$ and the first
pericenter $\sim 1$. The pericenter distance here is for the convenience of
comparison with Paper I. A closer encounter will lead to faster evolution, but
the general results should be similar. All the integrations are executed on
the $laohu$ GPU cluster in National Astronomical Observatories of China
(NAOC).

###  2.3 Scheme of Close Encounters with SMBHs

The interactions between merging SMBHs and the surrounding stars, especially
those stars close to SMBHs, are the focus of this work. In order to carefully
investigate these ”close encounters” we involve some special treatments. In
general, there are four kinds of close encounters in our model: normal stars
tidally disrupted by the SMBH, main sequence stars with relatively light mass
swallowed by the SMBH without tidal flares, giant branch stars partially
disrupted, and compact stars, such as NSs and BHs, forming EMRIs or the
extreme mass ratio bursts(EMRBs, more discussions can be found in Seciton
4.2).

Fig. 1 demonstrates what will happen if an ordinary star passes by a SMBH
within the tidal radius. Most of stars in our model are MS stars. If the
stellar mass of a star is massive enough, there will be a typical TDE.
Conversely, the tidal radius of a relatively light MS star or WD could be
smaller than the Schwarzschild radius. The star will be directly swallowed by
the SMBH without flare. The critical swallow radius of a BH, comparing to the
Schwarzschild radius $r_{\rm sch}$, can be enlarged when we take into account
the eccentricity of the intruder star. In this work we adopt the pericenter of
the marginally stable orbit $r_{\rm MSO}$ as the critical boundary. A star
will be marked as a plunge event when it’s tidal radius and separation to the
SMBH are both smaller than $r_{\rm MSO}$ (Cutler et al., 1994; Gair et al.,
2005)

$r_{\rm MSO}=\frac{3+e}{1+e}r_{\rm sch},$ (7)

where $r_{\rm sch}=2\rm{G}M_{\rm BH}/\rm{c}^{2}$ is the Schwarzschild radius,
$e$ is the eccentricity of the star. $r_{\rm MSO}$ determines the minimum
pericenter distance of a test particle with fixed eccentricity. A test
particle with pericenter distance smaller than $r_{\rm MSO}$ plunges directly
into the black hole. If $e=0$ the $r_{\rm MSO}$ is equivalent to the innermost
stable circular orbit(ISCO). For extreme hyperbolic situations with
$e\rightarrow\infty$, only orbits with pericenter distances larger than the
$r_{\rm sch}$ can survive. For simplicity, the mass and linear momentum of
both the disrupted and swallowed stars will directly add to the SMBH in our
model. Here we only take into account the linear momentum, and the angular
momentum is not considered because we can not trace the spin evolution of the
SMBH without the PN approximations included.

Figure 1: The fate of different types of stars after close encounters with
SMBHs.

In addition to main sequence stars, there are many giant branch stars (GBs),
such as RGs and AGBs. These stars, with their compact cores and diffuse
envelopes, usually have very large radii compared with MS stars, which
consequently corresponds to significantly large tidal radii. In most cases a
giant tidal disruption (GTD) will not cause a complete disruption. Inversely,
the stellar core, and even a large fraction of the envelope of the star can
survive (MacLeod et al., 2012). We assume that the envelope will be completely
striped and accreted by the SMBH immediately after the GTD, and the remnant
could be considered as a production of a fast evolved giant star. According to
the predicted evolution paths proposed by Hurley et al. (2000), we make a
simplified evolution scheme of GTDs in Fig. 1. Depending on the initial mass
of the disrupted star, GTDs in this scheme may result in different types of
remnants, such as naked helium main sequence stars, WDs and naked helium stars
in the Hertzsprung gap. More details can be found in the discussion of Hurley
et al. (2000).

Compact stars, on the other hand, are simply assumed to be directly swallowed
by the SMBH if their separations to the SMBH are smaller than the ISCO. In
that case, the linear momentum and mass of the compact stars will be added to
the SMBH. Therefore the capture criteria of compact stars is different from
other stars. Our integration does not include the post-Newtonian approximation
because it is very difficult to be involved in three body problems. The orbits
of those compact stars with heavier masses, especially when they are close to
ISCO, are not accurate enough. The eccentricity evolution due to GW emission
can not be well traced in this work. We prefer to leave this problem in future
works. For this reason we did not adopt the $r_{\rm MSO}$ as the criteria for
simplicity.

###  2.4 From numerical models to the reality

Direct N-body simulations are limited by particle resolution. Even taking into
account the acceleration of Heterogeneous Computing, the maximum particle
number that can be managed by direct N-body simulation with reasonable
integration time is only several million, which is obviously less than the
number of stars in a typical galaxy. A direct consequence of limited particle
resolution is that the TDEs are so rare in the simulation that we can not
collect enough events to make statistical analyses. For example, if there is a
galaxy with total mass $10^{9}\,{\rm M_{\odot}}$, and the half mass radius is
$r_{1/2}=1{\rm\,kpc}$, then a length scaling relation can be set up by Eq. 5.
For a solar type star, the tidal radius is $\sim 10^{-5}{\rm\,pc}$
corresponding to $\sim 10^{-8}$ in this simulation unit, which is a very tiny
scale in the simulation. The collected TDEs in such configuration will be only
a few in the entire simulation. In order to collect more TDEs in the
simulation, we have to adopt a larger ”tidal radius”, which is $r_{\rm t}\sim
10^{-4}$. And the Schwarzschild radius also adopts the same scaling to make
sure that $r_{\rm t}/r_{\rm sch}$ is the same both in simulation and reality.
By integrating several models with different tidal radii, we can try to
extrapolate the simulation result to reality. This scheme has been proved to
be feasible in Paper I. Since $r_{\rm t}\propto M_{\rm BH}^{1/3}$, while
$r_{\rm MSO}\propto M_{\rm BH}$, it is not recommended to make direct
extrapolations from simulation results to other real systems with different
$M_{\rm BH}$. Otherwise the results of the normal TDEs and those swallow
events could be inconsistent. For this reason, we have to fix the mass of the
SMBH in our simulations. In the following discussion, we set the total mass of
each galaxy/nucleus is $10^{9}\,{\rm M_{\odot}}$, and the mass of each SMBH is
$10^{7}\,{\rm M_{\odot}}$.

The limited particle resolution also induced other artificial effects. Some
dynamical processes are $N$ dependent. For instance, the two-body relaxation
timescale follows $t_{\rm r}\propto N/\ln N$, which means the two-body
relaxation in our model will be much faster than the case of a real galaxy. In
a spherical stellar system with a SMBH in the center, two-body relaxation is
the most important mechanism to scatter stars into the tidal disruption loss-
cone, the phase space region corresponding to the orbits of stars can be
tidally disrupted (Frank & Rees, 1976; Lightman & Shapiro, 1977). The process
that scatters stars into loss-cone is the so-called loss-cone refilling. The
situation in merging galaxies is quite different. In principle, as we did in
Paper I, we can divide the dynamical evolution of the SMBHB in a galaxy merger
into three phase. In phase I, two galaxies and their central SMBHs have so
large separation that the interaction between central SMBH and surrounding
stars is roughly the same as the case of a single SMBH in an isolated galaxy.
That means the loss-cone refilling will be dominated by two-body relaxation.
Thus the TDEs rate will be $N$ dependent in phase I. As two SMBHs get more and
more close to each other, the perturbation becomes stronger and stronger in
phase II. As a result, the perturbation will dominate the loss-cone refilling,
which will result in an approximate $N$ independent TDEs rate. In phase III,
two SMBHs form a compact binary and the loss-cone refilling is very
complicated. By carefully analyzing the different situations in different
phases, it is possible to make a credible extrapolation based on analytical
models and our numerical models with different $N$. Detailed discussions can
be found in Paper I and II. There are more discussions in Section 4.1.

In addition to the artificial effects mentioned above, the stellar evolution
model also needs to consider the influence of the limited particle resolution.
For instance, if we use $10^{6}$ particles to represent a dense star cluster,
the total mass of stars in the Kroupa model should be $\sim 6\times
10^{5}\,{\rm M_{\odot}}$. However, we have to adapt this cluster to represent
a galaxy with $10^{9}$ stars. That means every particle in our model actually
corresponds to a group of $1000$ stars in the real galaxy. A straightforward
discrepancy between our model and the real galaxy is that they have different
two-body relaxation timescale. For example, we can assume that all the stars
in the cluster have evolved for $1{\rm\,Gyr}$. In our model, due to the
limited particle number, the system could be well relaxed. While the same
thing should not happen in a real dense star cluster with $10^{9}$ stars.
Therefore we need to have a proper scaling relation between the stellar
evolution timescale in reality and the two-body relaxation timescale in our
model. A simplified solution is to make sure that the ratio between relaxation
timescale and stellar evolution time scale in the real dense star cluster
should be equal to the ratio in our model (Panamarev et al., 2019). This
scaling relation is essential for models with stellar evolution included (more
discussions can be found in Panamarev et al. (2019)). However, as we mentioned
before, we are focusing on the merging phase, which corresponds to a
relatively short period of time. Consequently, in dry mergers with all stars
that have already evolved for quite a long time, the stellar evolution during
the merger could be neglected. Instead of including stellar evolution, it is
more important to have an initial model with properly evolved stars before the
merger, which is adopted here.

##  3 Results

###  3.1 Dynamical evolution of SMBHs and surrounding stars

Compared with previous major merger models with equal mass stars, we have more
realistic models with different stellar components. Fig. 2 demonstrates the
dynamical evolution of SMBHs in merger galaxies, which are numerically
integrated by our models with $N=2\times 10^{6}$ particles. Unless otherwise
specified, all results in following discussions are based on this model. For
better comparison, we have two integrations with and without initial mass
function (IMF), which are represented by blue solid lines and red dashed
lines, respectively. The left panel represents the evolution of $d_{\rm{BH}}$,
the separation of two SMBHs, and the right panel demonstrates the evolution of
the semi-major axis $a$. Since bound SMBHBs are not formed in phase I, the
time begins from $t\sim 80$ in the figure. As the results indicated, though
the evolutions of two SMBHs in different models are roughly the same in phase
I, there is a significant difference in phase II. The SMBHB in the model with
IMF evolves slower than the case of the equal mass model, which looks
inconsistent with Khan et al. (2018). However, our integrations only continued
to $\sim 100$ N-body time unit, which is far less than Khan et al. (2018) did.
According to the Fig. 1 and Fig. 2 in their paper, the difference between
equal mass and multi-mass models is not obvious in the early stages of the
SMBHB formation.

Figure 2: The dynamical evolution of two SMBHs in merging galaxies. Red dashed
and blue solid lines correspond to equal stellar mass model and multi-mass
model, respectively. The left panel represents the separation evolution of two
SMBHs. The right panel represents the semi-major axis evolution of the SMBHB,
after a bound binary system is formed.

The evolution of Lagrangian radii is demonstrated in Fig. 3. A Lagrangian
radius is a sphere with the center of the stellar system. Stars inside this
sphere occupy a fixed fraction of the total mass of the system. By tracing the
evolution of several Lagrangian radii with different mass fractions, we can
trace the dynamical evolution of the star cluster. Here dashed lines and solid
lines denote equal mass and multi-mass models respectively. Different colors
correspond to $0.1\%$, $1\%$, $10\%$, $50\%$, and $90\%$ of the total mass
respectively, which have been marked in the legend. It should be noticed that
the $1\%$ curve is close to the influence radius, because the stellar mass
inside this radius is close to the mass of the central SMBH. The upper left
panel of the figure represents the Lagrangian radii relative to the center of
mass. Since two galaxies are far away from each other, there are only a few
stars around the center of mass, which results in the large Lagrangian radii
at the beginning. The upper right panel demonstrates the Lagrangian radii
relative to one of the SMBH. The inner region around the SMBH of the equal
mass model is slightly more compact than the multi-mass model, which is
consistent with a similar multi-mass model with single BH (Baumgardt et al.,
2004). The bottom panels are the corresponding average mass inside each
Lagrangian radii. In the equal mass model, all average masses are equal to
$10^{-6}$, the single particle mass in the simulation. In the multi-mass
model, the average masses inside larger radii, which correspond to all
Lagrangian radii equal and larger than $10\%$ in the figure, are very close to
$10^{-6}$. However, this is not the case in the inner region. According to
bottom panels, the average mass inside Lagrangian radii of mass fraction
$0.1\%$ and $1\%$ are significantly heavier than the average mass of the
entire system, which can be considered as the consequence of the mass
segregation at the inner region. In the bottom left panel, the average masses
of the inner region are close to $10^{-6}$ at the beginning of the
integration, because the center of mass is at the midpoint where the stellar
density is very low. The mass segregation effect is more significant in the
plot relative to one SMBH, which is demonstrated in the bottom right panel.
The inner region average masses are significantly higher than the equal mass
model at the beginning of the integration. However, during the formation of
the bound SMBHB, the central average masses decline sharply. The strong
interaction between the two SMBHs and the violent evolution of stars in the
central region suppress the mass segregation quickly. After the two SMBHs form
a compact binary, the average mass at the central region starts to recover
with gradually increasing.

Figure 3: The evolution of Lagrangian radii and corresponding average mass.
The upper left panel is the Lagrangian radii relative to the center of mass,
and the upper right panel corresponds to one of the SMBH. Dashed lines and
solid lines denote equal mass and multi-mass models respectively. The bottom
left and right panels demonstrate corresponding average mass inside each
Lagrangian radius. Different colors represent different mass fractions of the
total mass.

Fig. 4 demonstrates the evolution of Lagrangian radii and corresponding
average masses of different components relative to one of the SMBH. We
classified stars into three groups: normal stars (NORM, solid lines) including
MS star and WD, giant branch stars (GB, dotted lines) including RGs and AGBs
and compact stars (CPT, dashed lines) including NSs and BHs. Different colors
in the figure denote different mass fractions of corresponding component. Two
vertical gray lines divide the evolution into three phases. The criteria for
dividing three phases are based on the fluctuations of the TDE rate. Other
criteria, for instance, a criteria based on dynamical evolution, will only
slightly change the position and period of phase II, and will not lead to
significantly different results. More details can be found in Paper II. Since
NORM stars correspond to the largest fraction, the Lagrangian radii of NORM
stars are almost the same as the radii of all the stars. In the left panel,
CPT stars with heavier mass show significant mass segregation in the central
region, especially before the formation of SMBHB. Due to the heating of SMBHB,
all the Lagrangian radii in the central region expand after two SMBHs form a
bound system. Especially in the phase II, a large fraction of CPT stars in the
central region are kicked out or swallowed by the SMBHs, which result in a
significant jump of the Lagrangian radii inside the influence radius. This
result is consistent with Gualandris & Merritt (2012), who has numerically
investigated the evolution of SMBHBs in multi-component merging galaxies.
According to the results in the right panel, only CPT stars at central region
have significant evolution on average mass. NORM and GB stars roughly keep a
constant average mass during the integration, with latter corresponding to
heavier value. Similar to the Fig. 3, the average mass of CPT stars at inner
region has a sharp drop in phase II. Since the drops of other two components
are not significant, the sharp decline of the average mass in phase II should
be mainly contributed by CPT stars.

Figure 4: The evolution of Lagrangian radii and corresponding average masses
of different components relative to one of the SMBH. Solid lines, dotted lines
and dashed lines denote normal stars, GBs, and compact stars, respectively.
Different colors represent different mass fractions of the total mass of each
component. The entire evolution has been divided into three phases by two
vertical gray lines. The left panel is the evolution of Lagrangian radii, and
the right panel is the evolution of average mass inside each Lagrangian
radius.

###  3.2 Tidal disruptions in multi-mass models

In our multi-component simulation there are 13991 TDEs/swallow events within
the 200 N-body time unit, and $\sim 94\%$ are made by MS stars. The fraction
of GB, WD, NS and BH are, respectively, $\sim 2.3\%$, $\sim 2.7\%$, $\sim
0.7\%$ and $\sim 0.3\%$. Similar to the equal mass model, the tidal disruption
in multi-mass models also can be divided into three phases, which correspond
to before, during and after the formation of bound SMBHBs. Fig. 5 demonstrates
the total capture rate, which includes both disrupted stars and swallowed
normal stars or compact stars, in the multi-mass model and the equal mass
model. The red dashed and blue solid lines represent the result of equal mass
model and multi-mass model respectively. The left panel and middle panel are,
respectively, the evolution of the mass accretion rate and the event rate. In
the simulation, we assume that all the mass of the disrupted stars will be
accreted into the SMBH. That means, the evolution of the mass accretion rate
and the event rate in the equal mass models are the same thing. However, it is
different in multi-mass models. As a result, the peak mass accretion rates of
the equal mass and multi-mass model in phase II are at the same level. While
their event rates are different. Compared with the multi-mass model which has
peak event rate $\dot{N}=275/[\rm T]$, the equal mass model has significantly
higher peak event rate $\dot{N}=421.5/[\rm T]$. According to this result,
though the disrupted/swallowed stars are less than the case in the equal mass
model, the disrupted/swallowed stars in the multi-mass model are preferred at
the high mass end. It has been confirmed by the right panel of Fig. 5, which
represents the average mass of disrupted/swallowed stars. This effect is
mainly due to the larger tidal radii of heavier stars, and the mass
segregation may also have some contributions.

Figure 5: The evolution of the total tidal disruption/swallowed rate in the
multi-mass model and the equal mass model. The red dashed and blue solid lines
represent the result of equal mass model and multi-mass model respectively.
The left panel and middle panel are, respectively, the evolution of the mass
accretion rate and the event rate. The right panel is the average mass of
disrupted/swallowed stars.

####  3.2.1 Tidal disruptions of main sequence stars

Tidal disruption of MS stars dominate the TDEs in our models. As discussed in
Section 2.3, MS stars with relatively low mass generally correspond to small
radii. Stars with tidal radius smaller than $r_{\rm MSO}$ will be essentially
swallowed when they get close enough to the SMBH. According to our simulation
results with the largest particle number, there are $\sim 45\%$ MS stars, with
masses range from $\sim 0.1\,{\rm M_{\odot}}$ to $\sim 0.4\,{\rm M_{\odot}}$,
that have been swallowed by the SMBH. It should be noticed that all the WDs
will be swallowed by the SMBH too. But they only contribute less than $3\%$ of
the total disruption/swallow events. Other models with different particle
numbers give similar results. In general, nearly half of captured stars will
directly plunge into the SMBH without flares. Fig. 6 demonstrates the
disruption/swallow evolution of MS stars. The orange and green lines in the
figure denote stars with standard tidal disruption flare (STD) and stars with
plunge orbits without flare (PLG), respectively. Blue lines and gray lines
denote the evolution of all stars and the separation of two SMBHs,
respectively. The entire evolution is divided into three phases by two
vertical gray lines. The left y-axis represents the separation, and the right
y-axis in the left and right panels are, respectively, the mass accretion rate
and event number rate.

Fig. 6 indicates that, due to the strong perturbation and the rapid evolution
of the stellar distribution around two SMBHs, both swallowed stars and normal
disrupted stars have significantly enhanced rate in phase II. And there will
be a rate peak every time two SMBHs get close. This result is consistent with
the previous equal mass model in Paper I. From the right panel we can easily
see that the normal TDEs and the plunge cases have very close event rates.
However, it is obvious that the plunge stars only contribute a small fraction
of the mass accretion rate in the left panel, because they are dominated by
the stars at the low-mass end.

Figure 6: The tidal disruption/swallow evolution of MS stars. The left panel
and right panel are, respectively, the evolution of mass accretion rate and
the event rate. Blue, orange and green solid lines represent the result of all
of MS stars, standard TDE, and the MS stars which directly plunge into the
SMBH, respectively. Two vertical gray lines divide the evolution into three
phases.

####  3.2.2 Tidal disruptions of giant branch stars

GB stars usually have very large radii, which means a GB star has a very large
tidal radius. Usually, there is a dense core with significant mass
concentrated in a GB star. During the tidal disruption, instead of being
totally disrupted by the SMBH, only the envelope of the GB star will be
striped away, leading to a light accretion with only a fraction of the stellar
mass. As demonstrated in Fig. 7, the tidal disruption evolution of GB stars is
similar to MS stars. The evolution can be divided into three phases and there
is a significant enhanced rate in phase II. Here blue and orange lines
represent the evolution of CPT and GB stars respectively, the rest of the
legend and labels are the same as Fig. 6.

Figure 7: The tidal disruption evolution of CPT stars and GB stars. The left
panel and right panel are, respectively, the evolution of mass accretion rate
and the event rate. Blue and orange solid lines represent the result of CPT
and GB stars respectively. Two vertical gray lines divide the evolution into
three phases.

As mentioned in Section 2.3, after the disruption, the remnant stars usually
are WDs or naked helium MS stars. Most of them will fly by the SMBH and not be
swallowed. However, according to our result, there is $\sim 14\%$ remnants
will finally plunge into the SMBH. Some of them will quickly plunge into the
SMBH within one orbital period, and the rest may be captured by the SMBH first
and finally swallowed with a time delay $\Delta t$. Fig. 8 represents the
distribution of these accreted remnant stars. The x-axis and y-axis,
respectively, represent ”$t_{\rm{1stTD}}$”, the time that stars been tidally
disrupted by the SMBH and the time delay $\Delta t$ between the first tidal
disruption and the final plunge. The mass of remnants after the first tidal
disruption has been denoted in different colors, with a color bar in solar
mass. The size of the filled circles represents the mass of striped envelopes
during GTD with the range from $2.6\times 10^{-5}\,{\rm M_{\odot}}$ to
$1.8\,{\rm M_{\odot}}$. According to the figure, most of the plunge events of
the remnants are recorded soon after the first tidal disruption. There are
some exceptions which can survive for a relatively long time after the first
tidal disruption. And most of such large time delay cases are in early phases.
After a compact SMBHB formed in phase III, the remnants will be more easily
kicked by the companion SMBH and avoid the final plunge. Since our integration
terminated at $t=200$, there may be some plunge events at $t>200$ missed. In
addition, we also find that $\sim 17$ stars in Fig. 8 have been stripped off
more than half of their mass during the GTD.

Figure 8: The GB stars which have been finally swallowed by the SMBH after a
tidal disruption. The x-axis and y-axis represent the time that stars have
been tidally disrupted by the SMBH and the time delay between the tidal
disruption and the final plunge. Two vertical gray lines divide the evolution
into three phases. The color bar denotes the mass of the tidal disruption
remnants which finally plunge into the SMBH, in solar mass.

####  3.2.3 SMBHB and compact objects

Fig. 7 demonstrates the swallow rate evolution of CPT stars, which is similar
to GB stars. According to the heavier mass, though CPT plunge events are not
as frequent as GTD, their contributions on mass accretion is more significant
than GB stars. Fig. 9 represents the mass distribution of swallowed NSs and
BHs. Most swallowed NSs concentrate in the low mass end, while the BHs are the
opposite. It should be noted that, since general relativity effects are not
included in all of these integrations, the orbits of compact stars very close
to SMBHs are not accurate. In principle, by using post-Newtonian approximation
we can integrate the orbital evolution of a two-body system with relatively
high accuracy. However, such an approximation is not easy to be adopted in a
triple system which we are discussing here. Bonetti et al. (2016) has achieved
an improvement with corrections up to 2.5PN order. We prefer to solve this
problem with similar methods in future work. This simplification will not make
a significant influence on the TDR evolution. Because most of disrupted stars
are tend to be concentrated on the low mass end which was relatively less
affected. The rate may only be slightly underestimated, but the evolution
should be the same. While for CPT stars, specially those with orbits very
close to SMBHs, the absence of the PN corrections may lead to the suppressed
EMRIs formation.

Figure 9: The mass distribution of swallowed NSs and BHs. The left panel and
the right panel are, respectively, the distribution of swallowed NSs and BHs.
The masses of CPT stars are in solar mass and the y-axis represents the event
counts.

##  4 Discussion

###  4.1 Extrapolations to galaxies

As discussed in Section 2.4, the simulation results here can not be directly
used to estimate the tidal disruption rate evolution in a galaxy. Due to the
limited particle resolution, extrapolations based on the simulation results
are crucial. Here we follow the same scheme as Paper I. We choose $N=5\times
10^{5}$ for each galaxy and the average tidal radius $\overline{r}_{\rm
t}=5\times 10^{-4}$ as fiducial parameters. In order to investigate the
$r_{\rm t}$ and $N$ dependence, we vary $\overline{r}_{\rm t}$ with fixed
fiducial $N$, and vary $N$ with fiducial $\overline{r}_{\rm t}$ individually.
The $\overline{r}_{\rm t}$ varies in
$[10^{-5},10^{-4.5},10^{-4},10^{-3.5},10^{-3}]$, and the $N$ varies in
$[1.25\times 10^{5},2.5\times 10^{5},5\times 10^{5},10^{6}]$. In every model,
the average tidal disruption accretion rates and event rates in three phases
are calculated individually.

Fig. 10 demonstrates the corresponding simulation results of phase I and II
for different stellar components. The red, green and blue dots represent the
tidal disruption accretion rates of STD, PLG and CPT stars respectively.
Limited by the particle resolution, even the model with largest star particles
does not record enough events of GTD and CPT for statistical study in some
phases. Especially in phase II, due to the short time periods, only STD and
PLG stars have enough recorded events for statistical study. The left and
middle panels in the figure demonstrate the $\overline{r}_{\rm t}$ dependence
and $N$ dependence in phase I, respectively. And the right panel represents
the $\overline{r}_{\rm t}$ dependence in phase II.

Figure 10: The disruption rate extrapolations of phase I and II. The red,
green and blue dots represent the tidal disruption accretion rates of STD, PLG
and CPT events based on simulation results respectively. The red solid, the
green dashed and the blue dotted lines represent the corresponding fitting
results based on analytical estimations. The left and middle panels
demonstrate the $\overline{r}_{\rm t}$ dependence and $N$ dependence in phase
I, respectively. And the right panel denotes the $\overline{r}_{\rm t}$
dependence in phase II.

Based on the discussions in Paper I and II, the tidal disruption in phase I is
dominated by two-body relaxation, and the accretion rate can be well estimated
by

$\dot{M}\propto\left(\frac{N}{\ln\Lambda}\right)^{\alpha}\overline{r}_{\rm
t}^{\beta},$ (8)

where $\alpha$ and $\beta$ can be fitted through numerical simulation results,
and $\ln\Lambda$ is the Coulomb logarithm and can be estimated by

$\ln\Lambda\thickapprox\ln\left(\frac{M_{\rm BH}}{2\overline{m}_{*}}\right),$
(9)

where $\overline{m}_{*}$ is the average mass of stars (Preto et al., 2004).

However, unlike the equal stellar mass models, the mass accretion rates in
multi-mass models can not intuitively reflect the event rate. As a rough
estimation, we can simply assume $\dot{N}=\dot{M}/\overline{m}_{*}$. The TDE
rate in phase I can be estimated by

$\dot{N}\propto N\left(\frac{\ln\Lambda}{N}\right)^{\zeta}\overline{r}_{\rm
t}^{\xi},$ (10)

where $\zeta$ and $\xi$ can be fitted through numerical simulation results.

The tidal disruption loss cone refilling in phase II is almost full, which
means the tidal disruption mass accretion rate weakly depends on $N$.
Therefore the rate can be approximated as a power law of $\overline{r}_{\rm
t}$. But the event rate $\dot{N}$ should be divided by $\overline{m}_{*}$, or
equivalent to multiplying by $N$. The situation is very complicated in phase
III. Two-body relaxation, perturbation of companion SMBH and the triaxial
stellar distributions have significant contributions to the loss cone
refilling. It is very difficult to analytically estimate the disruption rate
in phase III. Therefore the extrapolation in phase III is not reliable.

With above analytical models, and assuming that different stellar components
roughly follow the same relations, we can make rough extrapolations for
different stellar components in phase I and II.

####  4.1.1 Phase I

Following the similar scheme in Paper I, combining with simulation results,
the tidal disruption mass accretion rate of STD stars in phase I can be
estimated as

$\dot{M}_{\rm STD}\sim
1.48\times\left(\frac{N}{\ln\Lambda}\right)^{-0.57}\overline{r}_{\rm
t}^{0.48},$ (11)

which can be extrapolated as

$\dot{M}_{\rm STD}\sim 2.7\times
10^{-2}\left(\frac{N}{\ln\Lambda}\right)^{-0.57}\left(\frac{r_{1/2}}{1{\rm\,kpc}}\right)^{-1.98}\left(\frac{\overline{r}_{\rm
t}}{10^{-6}\rm{pc}}\right)^{0.48}\rm M_{\odot}~{}yr^{-1}.$ (12)

The event rate is

$\dot{N}_{\rm STD}\sim
0.84N\left(\frac{N}{\ln\Lambda}\right)^{-0.57}\overline{r}_{\rm t}^{0.48},$
(13)

which can be extrapolated as

$\dot{N}_{\rm STD}\sim 1.5\times 10^{-2}\frac{\,{\rm
M_{\odot}}}{\overline{m}_{*}}\left(\frac{N}{\ln\Lambda}\right)^{-0.57}\left(\frac{r_{1/2}}{1{\rm\,kpc}}\right)^{-1.98}\left(\frac{\overline{r}_{\rm
t}}{10^{-6}\rm{pc}}\right)^{0.48}\rm yr^{-1}.$ (14)

For our fiducial model with $M=10^{9}\,{\rm M_{\odot}}$, $M_{\rm
BH}=10^{7}\,{\rm M_{\odot}}$, $r_{1/2}=1{\rm\,kpc}$, and according to our
simulation results with $\overline{m}_{*}=0.43\,{\rm M_{\odot}}$ corresponding
to stellar radius $r_{*}=0.40\rm R_{\odot}$, we can derive the mass accretion
rate of the disrupted stars should be $\sim 1.0\times 10^{-6}\rm
M_{\odot}~{}yr^{-1}$, and the event rate is $\sim 1.3\times
10^{-6}\rm{yr^{-1}}$.

Similarly, the mass accretion rate of PLG stars can be estimated as

$\dot{M}_{\rm PLG}\sim
0.16\times\left(\frac{N}{\ln\Lambda}\right)^{-0.47}\overline{r}_{\rm
t}^{0.51},$ (15)

with corresponding event rate

$\dot{N}_{\rm PLG}\sim
0.32N\left(\frac{N}{\ln\Lambda}\right)^{-0.47}\overline{r}_{\rm t}^{0.51}.$
(16)

The extrapolation of the mass accretion rate and the event rate are

$\dot{M}_{\rm PLG}\sim 1.6\times
10^{-3}\left(\frac{N}{\ln\Lambda}\right)^{-0.47}\left(\frac{r_{1/2}}{1{\rm\,kpc}}\right)^{-2.01}\left(\frac{\overline{r}_{\rm
t}}{10^{-6}\rm{pc}}\right)^{0.51}\rm M_{\odot}~{}yr^{-1},$ (17) $\dot{N}_{\rm
PLG}\sim 3.2\times 10^{-3}\frac{\,{\rm
M_{\odot}}}{\overline{m}_{*}}\left(\frac{N}{\ln\Lambda}\right)^{-0.47}\left(\frac{r_{1/2}}{1{\rm\,kpc}}\right)^{-2.01}\left(\frac{\overline{r}_{\rm
t}}{10^{-6}\rm{pc}}\right)^{0.51}\rm yr^{-1}.$ (18)

In our fiducial model, the mass accretion rate should be $\sim 4.0\times
10^{-7}\rm M_{\odot}~{}yr^{-1}$, and the event rate is $\sim 1.8\times
10^{-6}\rm{yr^{-1}}$.

GTD and CPT events are not popular in all simulations. GTD events are quite
rare. However, we can still find some CPT stars swallowed by SMBHs in phase I.
Although the data is not very good for statistical study, we can still try
some extrapolations, which gives the mass accretion rate

$\dot{M}_{\rm CPT}\sim
62.74\times\left(\frac{N}{\ln\Lambda}\right)^{-1.04}\overline{r}_{\rm
t}^{0.55},$ (19)

and the event rate

$\dot{N}_{\rm CPT}\sim
7.18N\left(\frac{N}{\ln\Lambda}\right)^{-1.04}\overline{r}_{\rm t}^{0.55}.$
(20)

The extrapolation relations of the mass accretion rate and the event rate are

$\dot{M}_{\rm CPT}\sim
0.29\left(\frac{N}{\ln\Lambda}\right)^{-1.04}\left(\frac{r_{1/2}}{1{\rm\,kpc}}\right)^{-2.05}\left(\frac{\overline{r}_{\rm
t}}{10^{-6}\rm{pc}}\right)^{0.55}\rm M_{\odot}~{}yr^{-1},$ (21) $\dot{N}_{\rm
CPT}\sim 3.3\times 10^{-2}\frac{\,{\rm
M_{\odot}}}{\overline{m}_{*}}\left(\frac{N}{\ln\Lambda}\right)^{-1.04}\left(\frac{r_{1/2}}{1{\rm\,kpc}}\right)^{-2.05}\left(\frac{\overline{r}_{\rm
t}}{10^{-6}\rm{pc}}\right)^{0.55}\rm yr^{-1}.$ (22)

And the corresponding accretion rate and the event rate of the fiducial model
are, respectively, $\dot{M}\sim 1.6\times 10^{-9}\rm M_{\odot}~{}yr^{-1}$ and
$\dot{N}\sim 4.4\times 10^{-10}\rm{yr^{-1}}$. It should be noticed that our
results here are purely Newtonian. Since the post-Newtonian approach is not
included in our model, the results here can only be considered as a crude
approximation.

Figure 11: The disruption event rate extrapolations of phase I and II. The
legends are similar to the Fig. 10. The first and the second rows demonstrate
the results of phase I and phase II, respectively. The left and right panels
in each row demonstrate the $\overline{r}_{\rm t}$ dependence and $N$
dependence, respectively.

Fig. 10 and 11 demonstrate the extrapolation fitting results of the mass
accretion rate and event rate respectively. Red, green and blue dots in two
figures represent simulation results of STD, PLG and CPT stars respectively.
Red solid, green dashed and blue dotted lines are, respectively, fitting
results of STD, PLG and CPT stars based on the derived extrapolation formulas
in Section 4.1.1 and 4.1.2. It should be aware that, due to very limited event
records, the extrapolation fitting results of swallowed CPT stars in PI are
not as good as STD or PLG stars. Compare the two figures, the STD and PLG have
similar event rates in phase I. But the STD stars correspond to significantly
higher mass accretion rates. Our extrapolation results of the fiducial model
also indicate that more than half events are contributed by PLG stars, which
may not induce any observational effects. But more than half accreted mass is
contributed by STD stars.

####  4.1.2 Phase II

Both analytical estimations and our numerical simulations indicate that the
tidal disruption rate in phase II does not significantly depend on $N$. The
accretion rate of STD stars can be estimated by

$\dot{M}_{\rm STD}\sim 0.11\overline{r}_{\rm t}^{0.82},$ (23)

with extrapolation relation

$\dot{M}_{\rm STD}\sim 2.2\times
10^{-6}\left(\frac{r_{1/2}}{1{\rm\,kpc}}\right)^{-2.32}\left(\frac{\overline{r}_{\rm
t}}{10^{-6}\rm{pc}}\right)^{0.82}\rm M_{\odot}~{}yr^{-1},$ (24)

which corresponds to $\dot{M}\sim 5.0\times 10^{-6}\rm M_{\odot}~{}yr^{-1}$
for our fiducial model.

The event rate is

$\dot{N}_{\rm STD}\sim 5.97\times 10^{-2}N\overline{r}_{\rm t}^{0.82},$ (25)

and the extrapolation can be wrote as

$\dot{N}_{\rm STD}\sim 1.2\times 10^{-6}\frac{\,{\rm
M_{\odot}}}{\overline{m}_{*}}\left(\frac{r_{1/2}}{1{\rm\,kpc}}\right)^{-2.32}\left(\frac{\overline{r}_{\rm
t}}{10^{-6}\rm{pc}}\right)^{0.82}\rm yr^{-1},$ (26)

which corresponds to $\dot{N}\sim 6.3\times 10^{-6}\rm{yr^{-1}}$ for the
fiducial model.

Similarly, the mass accretion rate of PLG stars can be estimated by

$\dot{M}_{\rm PLG}\sim 2.69\times 10^{-2}\overline{r}_{\rm t}^{0.85},$ (27)

with extrapolation relation

$\dot{M}_{\rm PLG}\sim 3.3\times
10^{-7}\left(\frac{r_{1/2}}{1{\rm\,kpc}}\right)^{-2.35}\left(\frac{\overline{r}_{\rm
t}}{10^{-6}\rm{pc}}\right)^{0.85}\rm M_{\odot}~{}yr^{-1}.$ (28)

For our fiducial model, that corresponds to $\dot{M}\sim 7.6\times 10^{-7}\rm
M_{\odot}~{}yr^{-1}$.

The event rate of PLG stars is

$\dot{N}_{\rm PLG}\sim 5.36\times 10^{-2}N\overline{r}_{\rm t}^{0.85}.$ (29)

and the extrapolation can be wrote as

$\dot{N}_{\rm PLG}\sim 6.5\times 10^{-7}\frac{\,{\rm
M_{\odot}}}{\overline{m}_{*}}\left(\frac{r_{1/2}}{1{\rm\,kpc}}\right)^{-2.35}\left(\frac{\overline{r}_{\rm
t}}{10^{-6}\rm{pc}}\right)^{0.85}\rm yr^{-1},$ (30)

which corresponding to $\dot{N}\sim 3.5\times 10^{-6}\rm{yr^{-1}}$ for the
fiducial model. Obviously, both the mass accretion rate and the event rate in
phase II have significant increases. The rate of STD has increased several
times compared with PI, which is consistent with the result in Paper I and II.

Since the period of phase II in all of our models is around 20 N-body time
units, there are not many recorded events contributed by CPT and GTD stars.
The extrapolation can not be well managed for them.

Table. 1 summarizes the extrapolation results of averaged mass accretion rates
and event rates in phase I and II, for disrupted/swallowed STD, PLG and CPT
stars respectively. Since the swallowed CPT stars in phase II are relatively
rare, we do not estimate their rates in the table.

Table 1: Averaged mass accretion rates and event rates in different phases.

Stage | $\dot{M}_{\rm STD}$ | $\dot{N}_{\rm STD}$ | $\dot{M}_{\rm PLG}$ | $\dot{N}_{\rm PLG}$ | $\dot{M}_{\rm CPT}$ | $\dot{N}_{\rm CPT}$
---|---|---|---|---|---|---
| $\times 10^{-6}\rm M_{\odot}/yr$ | $\times 10^{-6}\rm{/yr}$ | $\times 10^{-6}\rm M_{\odot}/yr$ | $\times 10^{-6}\rm{/yr}$ | $\times 10^{-9}\rm M_{\odot}/yr$ | $\times 10^{-9}\rm{/yr}$
(1) | (2) | (3) | (4) | (5) | (6) | (7)
PI | 1.0 | 1.3 | 0.4 | 1.8 | 1.6 | 0.4
PII | 5.0 | 6.3 | 0.8 | 3.5 | - | -

Note. — Col.(1): Stage of the evolution. Col.(2): Mass accretion rate of
disrupted STD stars. Col.(3): Event rate of disrupted STD stars. Col.(4): Mass
accretion rate of swallowed PLG stars. Col.(5): Event rate of swallowed PLG
stars. Col.(6): Mass accretion rate of swallowed CPT stars. Col.(7): Event
rate of swallowed CPT stars.

###  4.2 Delectability by space based GW instruments

As strong GW sources, SMBHBs could be directly detected by PTAs or space born
detectors such as LISA, TaiJi and TianQin in the future. There are many
discussions on this topic, which is out of the scope of this paper. Here we
want to discuss EMRIs, another kind of typical GW source produced by SMBHs
which could be detected by LISA/TaiJi/TianQin. Such events could also be
prompted around a SMBHB. However, a typical EMRI should spiral the SMBH for
many orbits. It needs accurate integrations with general relativity effects
considered, which is not included in our integrations. Therefore, it is not
reliable to estimate the rate of EMRIs based on our models. Another
interesting event are EMRBs, which are GW bursts from stellar objects passing
by a SMBH with very small pericenter distances. It has been considered as the
precursor to EMRIs, because many EMRBs will lose their energy and angular
momentum through GW radiation and finally evolve into EMRIs (Rubbo et al.,
2006). Hopman et al. (2007) demonstrates that the burst rates for stellar BHs
and MSs/WDs are, respectively, $1{\rm\,yr}^{-1}$ and $0.1{\rm\,yr}^{-1}$ in
the Milky Way. If extragalactic sources could be included, a detector like
LISA could manage the detection out to $\sim 100{\rm\,Mpc}$ for a $10\,{\rm
M_{\odot}}$ BH, with event rate of $\sim 0.2{\rm\,yr}^{-1}$ (Berry & Gair,
2013b, c). Recently, Han et al. (2020) calculates the event rates of very
extreme mass ratio bursts with a mass ratio about $10^{-8}$. They especially
considered the contribution of plunge stellar objects such as brown dwarfs
with unbound orbits. Their estimation indicates that, for small stellar
objects with mass $\sim 0.1\,{\rm M_{\odot}}$, the space based facilities
could detect the bursts inside $10{\rm\,Mpc}$, with corresponding event rate
$4-8{\rm\,yr}^{-1}$.

In our model, as discussed in Section 4.1, the rates of compact stars such as
NSs and BHs getting swallowed by SMBHs are quite low. However, there are many
low mass main sequence stars that plunge into SMBHs. According to the result
of our largest simulation, the average mass of swallowed stars is $\sim
0.2\,{\rm M_{\odot}}$, which have similar mass as Han et al. (2020) discussed.
In principal, the signal-to-noise ratio (SNR) of such kind of EMRBs is
proportional to $\rm m_{*}R^{-1}M_{\rm BH}^{2/3}$, with $R$ is the distance
from the source to observers (Berry & Gair, 2013b). With carefully numerical
estimations, Han et al. (2020) finds that the SNRs of plunging events with
$0.1\,{\rm M_{\odot}}$ stars in the Galactic Center can be up to ten thousands
for LISA. If we consider a $0.2\,{\rm M_{\odot}}$ star swallowed by an
$10^{7}\,{\rm M_{\odot}}$ SMBH, the SNR can be still as large as $\sim 8$ at
$\sim 50{\rm\,Mpc}$. However, according to estimations in Section 4.1, the
event rates of the PLG are only around $10^{-6}{\rm\,yr}^{-1}$ in phase I and
II. Since there are not too many galaxies inside $\sim 50{\rm\,Mpc}$, it is
unlikely to detect such kind of EMRBs with space borne GW detectors in the
near future.

##  5 Summary

We investigated full and partial tidal disruption of stars and direct plunges
into supermassive black holes (MSBH) in nuclear star clusters during and after
a galaxy merger. For that a full direct N-body simulation has been used with
stars obtained from a realistic stellar mass distribution, evolved for $\sim
1{\rm\,Gyr}$ (to account for the age of the galaxies before the merger),
surrounding two SMBH situated in the centres of the two merging galaxies.

With the stellar evolution included, there are different stellar components,
from low mass main sequence stars to white dwarfs, neutron stars and stellar
mass black holes. Different stellar objects, according to their properties
like mass and radius, have different fates after the close encounter with
central SMBHs. Compared to the equal mass model, the multi-mass model with
similar parameters tends to have similar or even slightly higher mass
accretion rate. However, although their event rates in phase I and III (before
and after the galaxy merger) are similar, the multi-mass model has lower event
rate compared to the equal mass model in phase II (the short time while the
merger is dynamically ongoing), which indicates that the disrupted/swallowed
stars in multi-mass model prefer high mass end.

In the multi-mass model, if a main sequence star is heavy enough, its close
encounter with a SMBH may lead to a standard TDE. Otherwise a light main
sequence star may correspond to a tiny tidal radius which makes it inside the
marginally bound orbits. That will result in a plunge event. During the galaxy
merger, due to the perturbation of the companion SMBH and stars around it,
both STD and PLG event rates will be enhanced in phase II. STD and PLG have
similar event rates. But STD events correspond to significantly higher mass
accretion rate, because the PLG events are preferring to low mass stars. Since
PLG events may not have any observational signatures, probably nearly half of
the SMBH tidal capture events are invisible.

Post sequence stars, such as RGs or AGB stars, could be partially disrupted.
There might be a core that survived after a GTD which stripped away the
envelope. The remnant could escape or plunge into the SMBH. Our largest
numerical simulation gets $\sim 10\%$ GTD stars finally plunge into SMBHs
after their disruption. Some of them can have bound orbits around the SMBH and
survive for many orbits.

CPT stars show significant mass segregation in the central region during phase
I. Due to the heating of the newly formed SMBHB, the Lagrangian radii in the
central region expand quickly during phase II. Since our integrations are
totally Newtonian, the orbits of compact stars close to the SMBH are not
accurate. According to our limited results, the rate of compact stars getting
swallowed by SMBHs in phase II could be significantly enhanced. And most
swallowed NSs concentrate on the low mass end, while BHs are the opposite.

We are grateful to the support of the National Natural Science Foundation of
China (NSFC11988101,NSFC11303039), the Key International Partnership Program
of the Chinese Academy of Sciences (CAS) (No.114A11KYSB20170015), and the
Strategic Priority Research Program (Pilot B) Multiwavelength gravitational
wave universe of CAS (No.XDB23040100). We (SL, PB, RS) acknowledge support by
CAS through the Silk Road Project at National Astronomical Observatories
(NAOC) of China, and the support by Key Laboratory of Computational
Astrophysics. The computations have been done on the Laohu supercomputer at
the Center of Information and Computing at NAOC, CAS, funded by the Ministry
of Finance of People’s Republic of China under the grant $ZDYZ2008-2$. LS and
RS acknowledge the support of Yunnan Academician Workstation of Wang Jingxiu
(No. 202005AF150025). LS acknowledges support from the K.C.Wong Education
Foundation. PB acknowledges the special support by the CAS President’s
International Fellowship for Visiting Scientists (PIFI) program during his
stay in NAOC, CAS. XC acknowledges the support of the National Natural Science
Foundation of China (No. 11873022). The work of PB was supported by the
Volkswagen Foundation under the special stipend No. 9B870 and the grant No.
97778. PB acknowledge the support within the grant No. AP14870501 of the
Science Committee of the Ministry of Science and Higher Education of
Kazakhstan. The work of PB was supported under the special program of the NRF
of Ukraine Leading and Young Scientists Research Support - ”Astrophysical
Relativistic Galactic Objects (ARGO): life cycle of active nucleus”, No.
2020.02/0346. PB thanks the support from the ACIISI, Consejería de Economía,
Conocimiento y Empleo del Gobierno de Canarias and the European Regional
Development Fund (ERDF) under grant with reference PROID2021010044.

## References

  * Amaro-Seoane (2018) Amaro-Seoane, P. 2018, Living Reviews in Relativity, 21, 4. doi:10.1007/s41114-018-0013-8
  * Arcavi et al. (2014) Arcavi, I., Gal-Yam, A., Sullivan, M., et al. 2014, ApJ, 793, 38. doi:10.1088/0004-637X/793/1/38
  * Arzoumanian et al. (2016) Arzoumanian, Z., Brazier, A., Burke-Spolaor, S., et al. 2016, ApJ, 821, 13. doi:10.3847/0004-637X/821/1/13
  * Arzoumanian et al. (2020) Arzoumanian, Z., Baker, P. T., Blumer, H., et al. 2020, ApJ, 905, L34. doi:10.3847/2041-8213/abd401
  * Babak et al. (2016) Babak, S., Petiteau, A., Sesana, A., et al. 2016, MNRAS, 455, 1665. doi:10.1093/mnras/stv2092
  * Baumgardt et al. (2004) Baumgardt, H., Makino, J., & Ebisuzaki, T. 2004, ApJ, 613, 1143
  * Begelman et al. (1980) Begelman, M. C., Blandford, R. D., & Rees, M. J. 1980, Nature, 287, 307
  * Belczynski et al. (2002) Belczynski, K., Kalogera, V., & Bulik, T. 2002, ApJ, 572, 407. doi:10.1086/340304
  * Berczik et al. (2005) Berczik, P., Merritt, D., & Spurzem, R. 2005, ApJ, 633, 680. doi:10.1086/491598
  * Berczik et al. (2006) Berczik, P., Merritt, D., Spurzem, R., & Bischof, H.-P. 2006, ApJ, 642, L21. doi:10.1086/504426
  * Berry & Gair (2013b) Berry, C. P. L. & Gair, J. R. 2013b, MNRAS, 433, 3572. doi:10.1093/mnras/stt990
  * Berry & Gair (2013c) Berry, C. P. L. & Gair, J. R. 2013c, MNRAS, 435, 3521. doi:10.1093/mnras/stt1543
  * Bonetti et al. (2016) Bonetti, M., Haardt, F., Sesana, A., et al. 2016, MNRAS, 461, 4419. doi:10.1093/mnras/stw1590
  * Chatterjee et al. (2003) Chatterjee, P., Hernquist, L., & Loeb, A. 2003, ApJ, 592, 32
  * Chen et al. (2009) Chen, X., Madau, P., Sesana, A., & Liu, F. K. 2009, ApJ, 697, L149
  * Chevalier (1989) Chevalier, R. A. 1989, ApJ, 346, 847. doi:10.1086/168066
  * Colgate (1971) Colgate, S. A. 1971, ApJ, 163, 221. doi:10.1086/150760
  * Colpi (2014) Colpi, M. 2014, Space Sci. Rev., 183, 189. doi:10.1007/s11214-014-0067-1
  * Coughlin et al. (2017) Coughlin, E. R., Armitage, P. J., Nixon, C., & Begelman, M. C. 2017, MNRAS, 465, 3840. doi:10.1093/mnras/stw2913
  * Cutler et al. (1994) Cutler, C., Kennefick, D., & Poisson, E. 1994, Phys. Rev. D, 50, 3816. doi:10.1103/PhysRevD.50.3816
  * Dehnen (1993) Dehnen, W. 1993, MNRAS, 265, 250
  * Evans & Kochanek (1989) Evans, C. R., & Kochanek, C. S. 1989, ApJ, 346, L13
  * Ferrarese & Merritt (2000) Ferrarese, L., & Merritt, D. 2000, ApJ, 539, L9
  * Foster & Backer (1990) Foster, R. S., & Backer, D. C. 1990, ApJ, 361, 300. doi:10.1086/169195
  * Frank & Rees (1976) Frank, J., & Rees, M. J. 1976, MNRAS, 176, 633. doi:10.1093/mnras/176.3.633
  * Gair et al. (2005) Gair, J. R., Kennefick, D. J., & Larson, S. L. 2005, Phys. Rev. D, 72, 084009. doi:10.1103/PhysRevD.72.084009
  * Gebhardt et al. (2000) Gebhardt, K., Bender, R., Bower, G., et al. 2000, ApJ, 539, L13
  * Gezari (2021) Gezari, S. 2021, ARA&A, 59, 21, doi:10.1146/annurev-astro-111720-030029
  * Gould & Rix (2000) Gould, A., & Rix, H.-W. 2000, ApJ, 532, L29. doi:10.1086/312562
  * Gualandris & Merritt (2008) Gualandris, A., & Merritt, D. 2008, ApJ, 678, 780. doi:10.1086/586877
  * Gualandris & Merritt (2012) Gualandris, A. & Merritt, D. 2012, ApJ, 744, 74. doi:10.1088/0004-637X/744/1/74
  * Guillochon & Ramirez-Ruiz (2013) Guillochon, J., & Ramirez-Ruiz, E. 2013, ApJ, 767, 25
  * Han et al. (2020) Han, W.-B., Zhong, X.-Y., Chen, X., et al. 2020, MNRAS, 498, L61. doi:10.1093/mnrasl/slaa115
  * Harfst et al. (2007) Harfst, S., Gualandris, A., Merritt, D., Spurzem, R., Portegies Zwart, S., & Berczik, P. 2007, New Astronomy, 12, 357. doi:10.1016/j.newast.2006.11.003
  * Hills (1975) Hills, J. G. 1975, Nature, 254, 295. doi:10.1038/254295a0
  * Hobbs et al. (2005) Hobbs, G., Lorimer, D. R., Lyne, A. G., et al. 2005, MNRAS, 360, 974. doi:10.1111/j.1365-2966.2005.09087.x
  * Hopman et al. (2007) Hopman, C., Freitag, M., & Larson, S. L. 2007, MNRAS, 378, 129. doi:10.1111/j.1365-2966.2007.11758.x
  * Huang et al. (2021) Huang, S., Hu, S., Yin, H., et al. 2021, ApJ, 920, 12. doi:10.3847/1538-4357/ac0eff
  * Hurley et al. (2000) Hurley, J. R., Pols, O. R., & Tout, C. A. 2000, MNRAS, 315, 543. doi:10.1046/j.1365-8711.2000.03426.x
  * Ivanov et al. (2005) Ivanov, P. B., Polnarev, A. G., & Saha, P. 2005, MNRAS, 358, 1361. doi:10.1111/j.1365-2966.2005.08843.x
  * Khan et al. (2018) Khan, F. M., Berczik, P., & Just, A. 2018, A&A, 615, A71. doi:10.1051/0004-6361/201730489
  * Khan et al. (2011) Khan, F. M., Just, A., & Merritt, D. 2011, ApJ, 732, 89. doi:10.1088/0004-637X/732/2/89
  * Komossa (2015) Komossa, S. 2015, Journal of High Energy Astrophysics, 7, 148. doi:10.1016/j.jheap.2015.04.006
  * Komossa & Bade (1999) Komossa, S., & Bade, N. 1999, A&A, 343, 775.
  * Komossa et al. (2003) Komossa, S., Burwitz, V., Hasinger, G., et al. 2003, ApJ, 582, L15. doi:10.1086/346145
  * Komossa et al. (2020) Komossa, S., Grupe, D., Parker, M. L., et al. 2020, MNRAS, 498, L35. doi:10.1093/mnrasl/slaa125
  * Kormendy & Ho (2013) Kormendy, J., & Ho, L. C. 2013, ARA&A, 51, 511
  * Kroupa (2001) Kroupa, P. 2001, MNRAS, 322, 231. doi:10.1046/j.1365-8711.2001.04022.x
  * Li et al. (2017) Li, S., Liu, F. K., Berczik, P., & Spurzem, R. 2017, ApJ, 834, 195
  * Li et al. (2019) Li, S., Berczik, P., Chen, X., Liu, F. K., Spurzem, R., & Qiu, Y. 2019, ApJ, 883, 132
  * Lightman & Shapiro (1977) Lightman, A. P., & Shapiro, S. L. 1977, ApJ, 211, 244. doi:10.1086/154925
  * Liu (2004) Liu, F. K. 2004, MNRAS, 347, 1357. doi:10.1111/j.1365-2966.2004.07325.x
  * Liu & Chen (2013) Liu, F. K., & Chen, X. 2013, ApJ, 767, 18. doi:10.1088/0004-637X/767/1/18
  * Liu & Wu (2002) Liu, F. K., & Wu, X.-B. 2002, A&A, 388, L48. doi:10.1051/0004-6361:20020566
  * Liu et al. (2009) Liu, F. K., Li, S., & Chen, X. 2009, ApJ, 706, L133
  * Liu et al. (2014) Liu, F. K., Li, S., & Komossa, K. 2014, ApJ, 786, 103
  * Luo et al. (2021) Luo, Z., Wang, Y., Wu, Y., et al. 2021, Progress of Theoretical and Experimental Physics, 2021, 05A108. doi:10.1093/ptep/ptaa083
  * MacLeod et al. (2012) MacLeod, M., Guillochon, J., & Ramirez-Ruiz, E. 2012, ApJ, 757, 134. doi:10.1088/0004-637X/757/2/134
  * Makino & Aarseth (1992) Makino, J., & Aarseth, S. J. 1992, PASJ, 44, 141
  * Magorrian et al. (1998) Magorrian, J., Tremaine, S., Richstone, D., et al. 1998, AJ, 115, 2285
  * Mei et al. (2021) Mei, J., Bai, Y.-Z., Bao, J., et al. 2021, Progress of Theoretical and Experimental Physics, 2021, 05A107. doi:10.1093/ptep/ptaa114
  * Merritt & Poon (2004) Merritt, D., & Poon, M. Y. 2004, ApJ, 606, 788
  * Mikkola & Valtonen (1992) Mikkola, S., & Valtonen, M. J. 1992, MNRAS, 259, 115
  * Milosavljević & Merritt (2003) Milosavljević, M., & Merritt, D. 2003, in American Institute of Physics Conference Series, Vol. 686, The Astrophysics of Gravitational Wave Sources, ed. J. M. Centrella, 201–210, doi:10.1063/1.1629432
  * Panamarev et al. (2018) Panamarev, T., Shukirgaliyev, B., Meiron, Y., et al. 2018, MNRAS, 476, 4224. doi:10.1093/mnras/sty459
  * Panamarev et al. (2019) Panamarev, T., Just, A., Spurzem, R., et al. 2019, MNRAS, 484, 3279. doi:10.1093/mnras/stz208
  * Preto et al. (2004) Preto, M., Merritt, D., & Spurzem, R. 2004, ApJ, 613, L109
  * Preto et al. (2011) Preto, M., Berentzen, I., Berczik, P., & Spurzem, R. 2011, ApJ, 732, L26. doi:10.1088/2041-8205/732/2/L26
  * Quinlan (1996) Quinlan, G. D. 1996, New A, 1, 35
  * Rees (1988) Rees, M. J. 1988, Nature, 333, 523
  * Ricarte et al. (2016) Ricarte, A., Natarajan, P., Dai, L., et al. 2016, MNRAS, 458, 1712. doi:10.1093/mnras/stw355
  * Rubbo et al. (2006) Rubbo, L. J., Holley-Bockelmann, K., & Finn, L. S. 2006, ApJ, 649, L25. doi:10.1086/508326
  * Saslaw et al. (1974) Saslaw, W. C., Valtonen, M. J., & Aarseth, S. J. 1974, ApJ, 190, 253
  * Shannon et al. (2015) Shannon, R. M., Ravi, V., Lentati, L. T., et al. 2015, Science, 349, 1522. doi:10.1126/science.aab1910
  * Shen et al. (2013) Shen, Y., Liu, X., Loeb, A., & Tremaine, S. 2013, ApJ, 775, 49. doi:10.1088/0004-637X/775/1/49
  * Shu et al. (2020) Shu, X., Zhang, W., Li, S., et al. 2020, Nature Communications, 11, 5876. doi:10.1038/s41467-020-19675-z
  * Tang et al. (2021) Tang, S., Silverman, J. D., Ding, X., et al. 2021, ApJ, 922, 83. doi:10.3847/1538-4357/ac1ff0
  * Tremaine et al. (2002) Tremaine, S., Gebhardt, K., Bender, R., et al. 2002, ApJ, 574, 740
  * Verbiest et al. (2016) Verbiest, J. P. W., Lentati, L., Hobbs, G., et al. 2016, MNRAS, 458, 1267. doi:10.1093/mnras/stw347
  * Vigneron et al. (2018) Vigneron, Q., Lodato, G., & Guidarelli, A. 2018, MNRAS, 476, 5312
  * Volonteri et al. (2003) Volonteri, M., Haardt, F., & Madau, P. 2003, ApJ, 582, 559
  * Wegg & Nate Bode (2011) Wegg, C., & Nate Bode, J. 2011, ApJ, 738, L8. doi:10.1088/2041-8205/738/1/L8
  * Yu (2002) Yu, Q. 2002, MNRAS, 331, 935
  * Zhang et al. (2008) Zhang, W., Woosley, S. E., & Heger, A. 2008, ApJ, 679, 639. doi:10.1086/526404
  * Zhong et al. (2014) Zhong, S., Berczik, P., & Spurzem, R. 2014, ApJ, 792, 137
  * Zhong et al. (2022) Zhong, S., Li, S., Berczik, P., et al. 2022, ApJ, 933, 96