# Exploring the origin of stars on bound and unbound orbits causing tidal
disruption events

Shiyan Zhong 11affiliation: Yunnan Observatories, Chinese Academy of Sciences,
396 Yang-Fang-Wang, Guandu District, 650216, Kunming, Yunnan, China , Kimitake
Hayasaki 22affiliation: Department of Astronomy and Space Science, Chungbuk
National University, Cheongju 361-763, Korea , Shuo Li 33affiliation: National
Astronomical Observatories of China, Chinese Academy of Sciences, 20A Datun
Rd., Chaoyang District, 100012, Beijing, China , Peter Berczik 33affiliation:
National Astronomical Observatories of China, Chinese Academy of Sciences, 20A
Datun Rd., Chaoyang District, 100012, Beijing, China 44affiliation: Main
Astronomical Observatory, National Academy of Sciences of Ukraine, 27
Akademika Zabolotnoho St., 03680, Kyiv, Ukraine 55affiliation: Astronomisches
Rechen-Institut, Zentrum für Astronomie, University of Heidelberg,
Mönchhofstrasse 12-14, 69120, Heidelberg, Germany , Rainer Spurzem
44affiliation: Main Astronomical Observatory, National Academy of Sciences of
Ukraine, 27 Akademika Zabolotnoho St., 03680, Kyiv, Ukraine 55affiliation:
Astronomisches Rechen-Institut, Zentrum für Astronomie, University of
Heidelberg, Mönchhofstrasse 12-14, 69120, Heidelberg, Germany 66affiliation:
Kavli Institute for Astronomy and Astrophysics, Peking University, Beijing,
China ††affiliationmark:

###### Abstract

Tidal disruption events (TDEs) probe properties of supermassive black holes
(SMBHs), their accretion disks, and the surrounding nuclear stellar cluster.
Light curves of TDEs are related to orbital properties of stars falling SMBHs.
We study the origin, density, and velocity distributions of bound and unbound
stars in the nuclear star cluster, which are causing TDEs as a function of
their orbital eccentricity $e$ and energy $E$. These quantities determine near
the SMBH the ratio of the orbit’s pericenter to tidal disruption radii
(denoted as penetration factor $\beta$). We develop an analytical model for
the density and velocity distribution of such stars in the cluster, which
agrees well with $N$-body experiments. Our model extends classical models of
angular momentum diffusion in the loss cone. We also derive an analytical
model for three characteristic eccentricities in the loss cone: the minimum
and maximum value for given $\beta$, respectively, and
$e_{\rm\leavevmode\nobreak\ lcb}$, which represents the orbital eccentricity
defining the boundary between empty and full loss cone regimes. With N-body
experiments, we show that stars causing TDEs are distributed between these
eccentricity limits on the $e-\beta$ plane. Moreover, we find most of the
bound stars between $e_{\rm\leavevmode\nobreak\ lcb}$ and $e=1$ (i.e., the
full loss cone regime), whereas the remaining bound stars are originating from
the empty loss cone regime. This is consistent with the loss cone theory. We
propose that the $e-\beta$ distribution of stars in a star cluster or galactic
nucleus can be a good tool to diagnose whether the stars can cause TDEs.

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

$\dagger$$\dagger$affiliationtext: Research Fellow at Frankfurt Institute for
Advanced Studies

##  1 Introduction

Most galaxies harbor supermassive black holes (SMBHs) with million to billion
solar masses at their center. Tidal disruption events (TDEs) provides a good
opportunity to identify dormant SMBHs in inactive galaxies. A star approaching
the SMBH is tidally destroyed, once it enters the tidal disruption radius,
where the star is torn apart by the black hole tidal force (Hills, 1975).
Mainly because the fallback material accretes onto the SMBH, a flaring event
occurs with large luminosity and a duration of months to years (Rees, 1988).

In the classical TDE theory, the star approaches the SMBH on a parabolic orbit
(i.e., $e=1$, where $e$ is the orbital eccentricity.) and is tidally disrupted
at $r_{\rm t}=(M_{\rm BH}/m_{*})^{1/3}r_{*}$, where $M_{\rm BH}$, $m_{*}$, and
$r_{*}$ are the black hole mass, stellar mass and radius, respectively. After
disruption, a half of the stellar debris is gravitationally bound and falls
back to the SMBH at a rate being proportional to $-5/3$ power of time because
of the Kepler’s third law, and thus the bolometric luminosity is proportional
to $t^{-5/3}$ (Evans & Kochanek, 1989; Phinney, 1989). However, recent
observations have revealed that some observed TDEs show light curves, which
deviate from the $t^{-5/3}$ decay rate (Gezari et al., 2012; Holoien et al.,
2014; Gezari et al., 2015; Miller et al., 2015; Holoien et al., 2016; van
Velzen et al., 2019). Dozens of X-ray TDEs have light curves shallower than
$t^{-5/3}$ (Auchettl et al., 2017), while many optical/UV TDEs are well fit by
$t^{-5/3}$ (e.g. Hung et al. 2017).

Lodato & Rossi (2011) have shown that the luminosity of the accretion disk
observed in different bands may decay with diverse power-law indexes. For TDEs
with the well-evolved optically thick disks, the observed X-ray light curves
should decay following the form of a power-law multiplied by an exponential,
which is caused by the Wien tail of the disk spectrum (Mummery & Balbus,
2020). Moreover, when mass falls back to the SMBH at super-Eddington rate, an
outflow could be launched from the accretion disk. The luminosity of the
outflow can decay with time more shallowly than $t^{-5/3}$ (Strubbe &
Quataert, 2009). The mass fallback rate can deviate from the $t^{-5/3}$ decay
rate because of the orbital eccentricity (Hayasaki et al., 2013, 2018; Park &
Hayasaki, 2020), the survived core through partial tidal disruption
(Guillochon & Ramirez-Ruiz, 2013), the stellar internal structure (Lodato et
al., 2009; MacLeod et al., 2012), and the penetration factor $\beta=r_{\rm
t}/r_{\rm p}$ (Park & Hayasaki, 2020), where $r_{\rm p}$ is the pericenter
distance radius.

Hayasaki et al. (2018) have examined the stellar distribution of a spherical
stellar cluster on the $e$-$\beta$ plane for a given number of particles and
tidal disruption radius by $N$-body experiments. They found that there are two
interesting features: first, there are two critical values of the orbital
eccentricity as a function of $M_{\rm BH}/m_{*}$ and $\beta$, indicating the
distribution of $e$ is coupled with the distribution of $\beta$. The light
curve characteristics of tidal disruption flares are divided by the two
critical eccentricities (see also Park & Hayasaki 2020). Second, the valid
region of the stars causing TDEs is constrained to a certain part of the
$e$-$\beta$ plane. Although the orbital eccentricity and the penetration
factor are important parameters to characterize the TDE flare, little has been
known how these orbital parameters are related to the stellar distribution and
dynamics in a spherical star cluster.

In this paper, we analytically derive the number densities of the bound and
unbound stars, which can cause TDEs, in a spherical cluster and test them by
performing $N$-body simulations. We also examine the distribution of the stars
on the $e$-$\beta$ plane. We construct our analytical models in Section 2, and
then describe the detail of $N$-body simulations in Section 3. We discuss our
results in Section 4 and draw our conclusions in Section 5.

##  2 Orbital parameter dependence of the stellar distribution

We first briefly review the theory of the stellar number density in a
spherical star cluster, which is proposed by Cohn & Kulsrud (1978) (hereafter
CK78). Then we derive the number densities of the stars on bound and unbound
orbits based on the CK78 solution. In what follows, we use subscript ‘b’ and
‘u’ to mark the quantities corresponding to bound and unbound cases,
respectively.

###  2.1 Stellar distributions in the energy-angular momentum phase space

When an SMBH resides in the center of a star cluster, a stellar cusp forms
within an influence radius of the SMBH, $r_{\rm h}$. A stellar orbit can be
determined by both the specific orbital energy of a star: $E=v^{2}/2-GM_{\rm
BH}/r$ and the normalized squared angular momentum of the star:

$\displaystyle\mathcal{R}=\left(\frac{J}{J_{\rm c}}\right)^{2},$ (1)

where $v$, $G$, $J$ and $J_{\rm c}$ are the velocity of the star, gravitation
constant, specific angular momentum of the star and the corresponding circular
angular momentum, respectively. In the spherically symmetric cusp with
isotropic velocity dispersion, the number density distribution of the stars
simply depends on the orbital energy. However, the isotropic velocity
distribution breaks down because the SMBH removes the stars with low angular
momentum through TDEs. In the classical loss cone theory, the number density
should also depend on $\mathcal{R}$, and thus should be a function of both $E$
and $\mathcal{R}$ as $n(\mathcal{R},E)$ (Frank & Rees, 1976). In phase space
the loss cone region is encompassed by $\mathcal{R}_{\rm lc}$, where
$\mathcal{R}_{\rm lc}$ is the square of the normalized loss cone angular
momentum (see equation 12). The stars inside the loss cone region can survive
for no more than one orbital period, unless they find a way out before being
disrupted. Thus the loss cone runs out of stars quickly and $n(\mathcal{R},E)$
vanishes to zero at $\mathcal{R}<\mathcal{R}_{\rm lc}$. On the other hand,
two-body encounters can replenish continuously new stars into the loss cone.
Consequently, $n_{\rm b}(\mathcal{R},E)$ is determined by the stellar
consumption and replenishment processes. Taken these two competing processes
into account, CK78 found an orbit-averaged number density by solving the
Fokker-Planck equation:

$n(\mathcal{R},E)\simeq A(E)\ln\left(\frac{\mathcal{R}}{\mathcal{R}_{\rm
0}}\right)\hskip 14.22636pt(\mathcal{R}>\mathcal{R}_{\rm 0}),$ (2)

where $A(E)$ is an energy-dependent normalization coefficient and
$\mathcal{R}_{\rm 0}$ is the square of the normalized angular momentum at the
zero-boundary below which the number density goes to zero.

Because the square of the normalized angular momentum should be in the range
of $\mathcal{R}_{\rm 0}\leq\mathcal{R}\leq\mathcal{R}_{\rm lc}$ for stars
which will enter inside the tidal disruption radius, the number density for a
given $E$ is given by

$n_{\rm TD}(E)=\int_{\mathcal{R}_{\rm 0}}^{\mathcal{R}_{\rm lc}}A_{\rm
TD}(E)\ln\left(\frac{\mathcal{R}}{\mathcal{R}_{\rm
0}}\right)\mathrm{d}\mathcal{R},$ (3)

where the normalization coefficient is estimated to be

$A_{\rm TD}(E)=\frac{n_{\rm TD}(E)}{\mathcal{R}_{\rm lc}\ln\mathcal{R}_{\rm
lc}-\mathcal{R}_{\rm lc}\ln\mathcal{R}_{\rm 0}-\mathcal{R}_{\rm
lc}+\mathcal{R}_{\rm 0}}.$ (4)

By introducing

$Q\equiv\frac{\Delta\mathcal{R}}{\mathcal{R}_{\rm lc}},$ (5)

where $\Delta\mathcal{R}$ is the cumulative change of $\mathcal{R}$ over one
orbital period of the star, Magorrian & Tremaine (1999) evaluated
$\mathcal{R}_{\rm 0}$ by

$\mathcal{R}_{\rm 0}=g(Q)\mathcal{R}_{\rm lc},$ (6)

where

$g(Q)=\left\\{\begin{array}[]{cc}\mathrm{exp}(-Q),&(Q>1)\\\
\mathrm{exp}(-0.186Q-0.824\sqrt{Q})&(Q<1)\end{array}\right.$ (7)

(see also Merritt 2013). Since $\mathcal{R}_{\rm 0}$ is very close to
$\mathcal{R}_{\rm lc}$ in the case of $Q<1$, the number density almost goes to
zero in the loss cone region. We call the $Q<1$ case as the empty loss cone
regime. On the other hand, $\mathcal{R}_{\rm 0}$ is smaller than
$\mathcal{R}_{\rm lc}$ in the $Q>1$ case and therefore $R_{0}\approx 0$ in the
limit of $\mathcal{R}_{0}\ll{R}_{\rm lc}$. We call this limit the full loss
cone regime. A significant part of the loss cone region is thus populated with
the stars, which are originating from a transition region between the empty
and full loss cone regimes. We will explain how to compute $Q$ in Section 3.

###  2.2 Number density of bound stars

The CK78 solution is solved for the stars moving inside the stellar cusp
surrounding an SMBH. Therefore, all of the stars are bound to the SMBH. In
this subsection, we first obtain the relation between the orbital eccentricity
and the square of the normalized angular momentum, and then we derive the
number density for the bound stars in the range of
$\mathcal{R}\leq\mathcal{R}_{\rm lc}$.

The specific angular momentum of the bound stars, $J$, is given by

$J=\sqrt{GM_{\rm BH}a(1-e^{2})},$ (8)

where $e<1$ and $a$ is the semi-major axis of the star. Putting $e=0$ into the
above equation, we get the circular angular momentum $J_{\rm c}=\sqrt{GM_{\rm
BH}a}$. At the $e=1$ limit,

$J\approx\sqrt{2GM_{\rm BH}r_{\rm p}},$ (9)

where $r_{\rm p}=a(1-e)$. Equating $r_{\rm p}$ with $r_{\rm t}$, we get the
loss cone angular momentum:

$J_{\rm lc}=\sqrt{2GM_{\rm BH}r_{\rm t}}.$ (10)

Substituting equations (9) and (10) into equation (1), we obtain

$\mathcal{R}=1-e^{2}$ (11)

and

$\mathcal{R}_{\rm lc}=\frac{2r_{\rm t}}{a}.$ (12)

Since $n_{\rm TD,b}(\mathcal{R},E)\mathrm{d}\mathcal{R}=n_{\rm
TD,b}(e,E)\mathrm{d}e$ due to the mass conservation law, equations (2), (6),
and (11) provide the number density of the bound stars in the range of $e_{\rm
ll}\leq e\leq e_{\rm ul}$:

$n_{\rm TD,b}(e,E)=2A_{\rm
TD,b}(E)e\ln\left(\frac{1-e^{2}}{g(Q)\mathcal{R}_{\rm lc}}\right),$ (13)

where

$e_{\rm ll}=\sqrt{1-2\frac{r_{\rm t}}{a}}$ (14)

comes from the $\mathcal{R}=\mathcal{R}_{\rm lc}$ limit and

$e_{\rm ul}=\sqrt{1-2g(Q)r_{\rm t}/a}$ (15)

is obtained by the $\mathcal{R}=\mathcal{R}_{\rm 0}$ limit. Note that $n_{\rm
TD,b}(e,E)$ goes to zero outside of this eccentricity range.

Next, let us obtain the number density of the bound stars as a function of
$\beta$ and $E$. From equations (9) and (10), $r_{\rm t}$ and $r_{\rm p}$ can
be expressed as

$r_{\rm p}=\frac{J^{2}}{2GM_{\rm BH}},$ (16)

and

$r_{\rm t}=\frac{J_{\rm lc}^{2}}{2GM_{\rm BH}},$ (17)

respectively. Then, we get the relation between $\beta$ and $\mathcal{R}$
through equation (1)

$\beta=\frac{\mathcal{R}_{\rm lc}}{\mathcal{R}}$ (18)

because of $\beta=r_{\rm t}/r_{\rm p}$. Since $n_{\rm
TD,b}(\mathcal{R},E)\mathrm{d}\mathcal{R}=n_{\rm
TD,b}(\beta,E)\mathrm{d}\beta$ for a given $E$, we obtain

$n_{\rm TD,b}(\beta,E)=\frac{A_{\rm TD,b}(E)\mathcal{R}_{\rm
lc}}{\beta^{2}}(-\ln g(Q)-\ln\beta)$ (19)

in the range of $1\leq\beta\leq 1/g(Q)$, where we used equations (2) and (6)
for the derivation. Note that the number density vanishes to zero at
$\beta=1/g(Q)$, which corresponds to that $\mathcal{R}$ equals
$\mathcal{R}_{\rm 0}$ in the original number density (see equation 2).

###  2.3 Number density of unbound stars

Most of the theoretical works have so far focused on bound stars within the
cusp (Cohn & Kulsrud, 1978; Shapiro & Marchant, 1978; Magorrian & Tremaine,
1999; Wang & Merritt, 2004; Stone & Metzger, 2016). In this section, we derive
the number density of the unbound stars outside of the cusp. The loss cone
effect should be negligible at the outside of the influence radius of the SMBH
because of $r_{\rm t}\ll{r}_{\rm h}$, so that the stars follow the Maxwellian
velocity distribution (we refer this case as the Maxwellian regime in the rest
of the paper). In the Maxwellian regime, the number density of the unbound
stars is proportional to $\beta^{-2}$ (Rees, 1988):

$n_{\rm TD,u}(\beta,E)=A_{\rm TD,u}(E)\frac{1}{\beta^{2}},$ (20)

in the range of $\beta\geq 1$.

Since the pericenter distance is given by $r_{\rm p}=a(e-1)$ for the stars on
hyperbolic orbits, we obtain

$e=1+\frac{E}{|E_{\rm t}|}\frac{1}{\beta},$ (21)

through $E=GM/(2a)$. This equation gives us the number density as a function
of only $E$ through $n_{\rm TD,u}(\beta,E)\mathrm{d}\beta=n_{\rm
TD,u}(e,E)\mathrm{d}e$ as

$n_{\rm TD,u}(e,E)=A_{\rm TD,u}(E)\frac{|E_{\rm t}|}{E}$ (22)

in the range of $1<e\leq 1+r_{\rm t}/a$. Note that $n_{\rm TD,u}(e,E)=0$ if
$e>1+r_{\rm t}/a$.

##  3 $N$-body experiments

In this section, we test the number density distributions that we have
analytically derived by $N$-body simulations. For comparison purposes, we
first compute the simulated number density as a function of $e$ or $\beta$ for
a given $E$ because of the low resolution. To do so, we integrate equations
(13) and (19) over the whole energy range as

$n_{\rm TD,b}(e)=\int_{E_{\rm t}}^{0}n_{\rm TD,b}(e,E)\mathrm{d}E,$ (23)
$n_{\rm TD,b}(\beta)=\int_{E_{\rm t}}^{0}n_{\rm TD,b}(\beta,E)\mathrm{d}E,$
(24)

where $E_{\rm t}=-GM_{\rm BH}/(2r_{\rm t})$ is the specific energy estimated
at the tidal disruption radius. For the number densities of the unbound stars,
we obtain $n_{\rm TD,u}(e)$ and $n_{\rm TD,u}(\beta)$ in the same way but
integrate from $E=0$ to infinity.

There are the following procedures for us to evaluate $Q$ for a given $E$:
First we estimate $E$ at the apocenter distance radius, $r_{\rm a}$. Since the
number densities of unbound stars are independent of $Q$ (see equations 20 and
22), we do not need to define $r_{\rm a}$ for the unbound stars. Because
$r_{\rm a}=(1+e)a\approx 2a$ for the bound stars on very highly eccentric
orbits, we obtain $r_{\rm a}\approx-GM_{\rm BH}/E$ for a given $E$. In the
case that $r_{a}$ is larger than the influence radius $r_{\rm h}$, we set
$r_{\rm a}\approx{r}_{\rm h}$ because the bound stars should be, by
definition, inside the $r_{\rm h}$.

Next, we compute the averaged change-per-orbit of the angular momentum due to
the two-body encounters at for a given $r_{a}$ as $\langle\Delta J\rangle^{2}$
because most of the angular momentum of the stars are likely to change at the
apocenter passage. This assumption can be justified because the orbiting star
passes its apocenter so slowly that has more time to interact with the
surrounding stars, and also, the perturbing forces may exert a non-negligible
torque to the passing star (Touma & Tremaine, 1997; Zhong et al., 2015). From
the relaxation theory (Frank & Rees, 1976; Merritt, 2013), we get

$\langle\Delta J\rangle^{2}\approx J_{\rm K}(r_{\rm a})^{2}\left(\frac{t_{\rm
dyn}}{t_{\rm relax}}\right),$ (25)

where $J_{\rm K}(r)=\sqrt{GMr}$ is the Keplerian specific angular momentum at
$r$, $t_{\rm dyn}=r_{\rm a}/\sigma(r_{\rm a})$ is the dynamical timescale,

$t_{\rm relax}=\frac{0.065\sigma^{3}(r_{\rm a})}{G^{2}m\rho(r_{\rm
a})\ln\Lambda}.$ (26)

is the local relaxation time (Spitzer, 1987), $\sigma(r)$ is the velocity
dispersion of the stars, $\rho(r)$ is the radial density profile of the star
cluster, $m=M_{\rm c}/N$ is the mass of the star, and $\ln\Lambda=\ln(0.11N)$
is the Coulomb logarithm (Giersz & Spurzem, 1994a). Final we can evaluate $Q$
from equations (5), (12) and (25) as

$Q=\frac{\Delta{R}}{R_{\rm lc}}\approx\left(\frac{t_{\rm dyn}}{t_{\rm
relax}}\right)\left(\frac{a}{r_{\rm t}}\right).$ (27)

We find from the above equation that $Q$ is a function of only energy, for the
given density profile and velocity dispersion, so that $g(Q)$ becomes thus a
function of only $E$ because of equation (7).

###  3.1 Basic model

We compare the analytically derived number densities with those obtained by
$N$-body simulations. The spherical star cluster consisting of $N$ equal mass
stars is initialized with a Plummer model, with an SMBH fixed at the center.
The quantities are expressed with the dimensionless Hénon unit, in which
$G=M_{\rm c}=1$ and the total energy of the system is $E=-1/4$ (Heggie, 2014a,
b). In the simulations, we take $r_{\rm t}$ as a fixed accretion radius, in
which all the stars are regarded as being tidally disrupted and removed from
the simulations. More details can be also seen in Hayasaki et al. (2018).
There are two simulation models with $N=512$K and $r_{\rm t}=10^{-5}$. In the
two models, only the black hole mass is different: While one model has $M_{\rm
BH}=0.01$, another one has $M_{\rm BH}=0.05$. These models are the most-close-
to-reality star cluster model in the $N$-body simulation performed by Hayasaki
et al. (2018).

Figure 1 shows the dependence of the number densities on the orbital energy
for the stars on both the bound and unbound orbits around the black hole with
$M_{\rm BH}=0.01$. In each panel, the red histogram shows the simulated number
density, while the red line is a fitted curve by the double broken power-law
function. The stars are distributed between $E_{\rm min}$ ($<0$) and $E_{\rm
max}$ ($>0$). Note that $E_{\rm min}$ (roughly $-2$ in model unit) is much
larger than $E_{\rm t}$ ($-500$ in model unit). This is consistent with the
loss cone theory that the stars are originating far from the tidal disruption
radius. We also employ a double broken power-law function to manually fit the
simulated number density.

Figure 2 depicts the density profile and the velocity dispersion profile of
the star cluster of our models. We use the broken power-law functions to fit
these profiles. These functions are used to compute $Q$ of the number
densities of the bound stars so that we can extrapolate the density profile
and velocity dispersion to the influence radius of the SMBH. It is estimated
to be roughly $0.2$ in our simulations.

Figure 1:  Energy dependence of the number densities for both the bound (left
panel) and unbound (right panel) stars with $M_{\rm BH}=0.01$ cases. In each
panel, the red histogram shows the simulated number density, while the black
dotted line is a fitted curve by the double broken power-law function.

Figure 2:  Radial profiles of the density and of the square of the stellar
velocity dispersion. The black dotted line represents the simulated density
profile in left panel, while it shows the square of the simulated stellar
velocity dispersion in the right panel. In both panels, the red line
represents the fitted curve by the broken power-law function. Note that the
red line extends to the influence radius of the SMBH, which corresponds to
$\sim 0.2$, in each panel.

###  3.2 Number densities of the bound stars

Figure 3:  Dependence of the number density of the bound stars on the orbital
eccentricity. The red histogram and black dotted line represent the simulated
and theoretical number densities. The left panel corresponds to $M_{\rm
BH}=0.01$ case, while the right panel does to $M_{\rm BH}=0.05$ case.

We compare the theoretical number densities with the simulated ones. Figure 3
shows the dependence of the number density of the bound stars on the orbital
eccentricity. The red histogram and black dotted line represent the simulated
and theoretical number densities. In both panels, the simulated number
densities are in good agreement with the theoretical ones, except for around
$0.998$ or the lower eccentricity regime. This deviation of the lower
eccentricity regime is mainly caused by a very low particle resolution of the
simulations.

The number density of the $M_{\rm BH}=0.01$ case is wider for the orbital
eccentricity than that of the $M_{\rm BH}=0.05$ case. This trend can be
interpreted as follows: from equation (14), the lowest eccentricity of the
bound stars can be estimated to be $e_{\rm ll}(M_{\rm BH},E_{\rm
min})=\sqrt{1-4r_{\rm t}|E_{\rm min}|/GM_{\rm BH}}$. We find that $E_{\rm
min}=-2$ in the $M_{\rm BH}=0.01$ case, whereas $E_{\rm min}=-4$ (ignoring the
isolated bins) in the $M_{\rm BH}=0.05$ case. Substituting each quantity into
the above equation, we find $e_{\rm ll}(0.01,-2)=0.9960$ and $e_{\rm
ll}(0.05,-4)=0.9984$. These evaluations are consistent with the number density
distributions shown in Figure 3.

Figure 4:  $\beta$-dependence of the number density of the bound stars. The
red histogram and black dotted line represent the simulated and theoretical
number densities, respectively. The left panel is the plot of $M_{\rm
BH}=0.01$ case, while the right panel is that of $M_{\rm BH}=0.05$ case.

Figure 4 shows the dependence of the penetration factor on the number density
of the bound stars. The figure format of the two panels is the same as Figure
3. We find from the panels that the analytical number densities are good
agreement with the simulated ones within the range of $\beta\lesssim 10$. For
$\beta\gtrsim 10$, the deviation between the analytical and simulated number
densities gets larger because of the poor resolution regime.

###  3.3 Number density of the unbound stars

Figure 5: Same format as Figure 3 but for the unbound stars.

Figure 6: Same format as Figure 4 but for the unbound stars.

Figures 5 and 6 compare the analytical number densities of $n_{\rm TD,u}(e)$
and $n_{\rm TD,u}(\beta)$ with those of the simulated number densities. Our
theoretical predictions match well with the simulated number densities in both
figures. As in the bound star case, the number densities of $M_{\rm BH}=0.01$
case are more widely distributed over the eccentricity than the $M_{\rm
BH}=0.05$ case. Substituting $a_{\rm max}=GM_{\rm BH}/(2E_{\rm max})$ into
equation (21) with $\beta=1$, we obtain $e_{\rm max}(M_{\rm BH},E_{\rm
max})=1+2r_{\rm t}E_{\rm max}/(GM_{\rm BH})$. We find from the simulated value
of $E_{\rm max}$ that $e_{\rm max}(0.01,1.259)=1.0025$ and $e_{\rm
max}(0.05,1.585)=1.0006$. These suggest that the number density is more widely
distributed over the orbital eccentricity in the star cluster with the less
massive black hole.

###  3.4 Distribution of the stars on the $e$-$\beta$ plane

Hayasaki et al. (2018) proposed that there are two critical eccentricities to
classify TDEs into five different types by the stellar orbit: eccentric
($e<e_{\rm crit,1}$), marginally eccentric ($e_{\rm crit,1}\lesssim{e}<1$),
purely parabolic ($e=1$), marginally hyperbolic ($1<e<e_{\rm crit,2}$), and
hyperbolic ($e\gtrsim{e_{\rm crit,2}}$) TDEs, respectively, where $e_{\rm
crit,1}=1-2(M_{\rm bh}/m_{*})^{-1/3}\beta^{-1}$ and $e_{\rm crit,2}=1+2(M_{\rm
bh}/m_{*})^{-1/3}\beta^{-1}$. They found that stars on marginally eccentric
and hyperbolic orbits are the main TDE source in a spherical star cluster by
N-body experiments. Moreover, these two critical eccentricities have recently
modified by taking account of the effect of the penetration factor on the
tidal potential as $e_{\rm crit,1}=1-2(M_{\rm bh}/m_{*})^{-1/3}\beta^{k-1}$
and $e_{\rm crit,2}=1+2(M_{\rm bh}/m_{*})^{-1/3}\beta^{k-1}$, where
$0\leq{k}\lesssim 2$ (Park & Hayasaki, 2020). In our simulation case, we find
that they are reduced to be $e_{\rm ctit,1}=1-2(NM_{\rm
bh})^{-1/3}\beta^{k-1}$ and $e_{\rm ctit,2}=1+2(NM_{\rm
bh})^{-1/3}\beta^{k-1}$, respectively. In this section we give new, tighter
constraints on the $e$-$\beta$ distribution of the stars.

For a given energy $E=-GM_{\rm BH}/2a$ at the $\mathcal{R}=\mathcal{R}_{0}$
limit, we can define the minimum orbital eccentricity of the bound stars by
using $r_{\rm t}/\beta=a(1-e)$ as

$e_{\rm min}=1-\frac{2}{\beta}\frac{r_{\rm t}|E|}{GM_{\rm BH}},$ (28)

where $\beta=1/g(Q)$ is obtained through equations (6) and (18). Adopting
$\beta=1$, we find equation (28) corresponds to equation (14) at $r_{\rm
t}/a\ll 1$: $e_{\rm ll}=\sqrt{1-4r_{\rm t}|E|/(GM_{\rm BH}})\approx 1-2r_{\rm
t}|E|/(GM_{\rm BH})$.

For the unbound stars we use another condition that all the unbound stars have
$E\leq E_{\rm max}$. Substituting $a_{\rm max}=GM_{\rm BH}/(2E_{\rm max})$
into equation (21), we obtain

$e_{\rm max}=1+\frac{2}{\beta}\frac{r_{\rm t}E_{\rm max}}{GM_{\rm BH}}.$ (29)

Adopting $Q=1$, we obtain the orbital eccentricity on the boundary between the
empty loss cone and the full loss cone regimes:

$e_{\rm lcb}=1-\frac{r_{\rm t}}{a_{\rm lcb}\beta},$ (30)

where the corresponding semi-major axis, $a_{\rm lcb}$, can be obtained by
solving equation (27) with $Q=1$ semi-analytically. Here, we analytically
estimate it as

$a_{\rm lcb}=\left(\frac{t_{\rm relax}}{t_{\rm dyn}}\right)r_{\rm
t}=\left(\frac{0.13\pi}{\ln\Lambda}\frac{M_{\rm BH}}{m}\right)r_{\rm t}\approx
4.1\times 10^{4}r_{\rm
t}\,\left(\frac{\ln\Lambda}{10}\right)^{-1}\left(\frac{M_{\rm
BH}/m}{10^{6}}\right)$ (31)

simply by assuming that $\rho(r)=\sigma^{2}/(2\pi{G}r^{2})$ and
$\sigma=\sqrt{GM/r}$ are the density of the singular isothermal sphere and the
Keplerian velocity, respectively (e.g., Merritt 2013). Note also that we adopt
$r_{\rm a}=2a$ here. In this case, we find $e_{\rm lcb}\approx 1$ for
$\beta\geq 1$.

Figure 7:  Distribution of the stars, which can cause TDEs, on the $e$-$\beta$
plane. The solid blue, solid black, dashed blue lines denote $e_{\rm min}$,
$e_{\rm max}$, and $e_{\rm lcb}$, respectively (see also equations 28-30).
Note that $e_{\rm crit,1}\approx 0.885$ and $e_{\rm crit,2}\approx 1.12$ for
$M_{\rm bh}=0.01$, $\beta=1$, and $k=0$ case and $e_{\rm crit,1}\approx 0.933$
and $e_{\rm crit,2}\approx 1.067$ for $M_{\rm bh}=0.05$, $\beta=1$, and $k=0$
case. Therefore, the two critical eccentricities are much larger or smaller
than the current eccentricity scales.

Figure 7 shows the distribution of the bound and unbound stars, which can
cause TDEs, on the $e$-$\beta$ plane. The solid blue, solid black, dashed blue
lines denote $e_{\rm min}$, $e_{\rm max}$, and $e_{\rm lcb}$, respectively. We
evaluate these three characteristic eccentricities by using $E_{\rm min}$,
$E_{\rm max}$, and $a_{\rm lcb}$, which are obtained from our $N$-body
simulations. While the space between $e_{\rm lcb}$ and $e_{\rm min}$
corresponds to the empty loss cone regime, the space between $e_{\rm lcb}$ and
$e=1$ corresponds to the full loss cone regime. The stars located above the
$e=1$ axis are supplied to the black hole from the Maxwellian distribution
regime. We find that only a small fraction of the stars are originated from
the empty loss cone regime and the corresponding values of $\beta$ are
distributed around unity, whereas most of the stars are originally supplied
from the region between $e_{\rm lcb}$ and $e_{\rm max}$ over the much wider
range of $\beta$. These results are consistent with the loss cone theory: due
to the diffusion nature in the empty loss cone regime, the loss cone flux is
much smaller than the full loss cone regime, and the stars inside the empty
loss cone cannot penetrate the tidal radius too much.

We also find that there are some stars outside $e_{\rm min}$ or $e_{\rm max}$.
These outliers seem to violate the loss cone theory. There are two possible
reasons for the outliers that unexpectedly behave: first, the more energetic
two-body encounter can occur in the $N$-body simulations, leading to the
enhancement of the angular momentum exchange so that the stars have
$\mathcal{R}<\mathcal{R}_{0}$. Second, the quantity $Q$ depends on the density
and the velocity dispersion of the star cluster, which can fluctuate with
radius and evolve with time. Therefore we need to reevaluate $e_{\rm min}$ and
$e_{\rm max}$ by taking account of these effects, but they are beyond the
scope of the CK78 solution. Both $e_{\rm min}$ and $e_{\rm max}$ are thus
likely to deviate from the correct value, producing some outliers on the
$e$-$\beta$ plane.

##  4 Discussion

Our analytical models have been tested by $N$-body simulations in the previous
section. Our model can be applied to the observed spherical star cluster
systems if we obtain the density and velocity dispersion profile of the star
cluster from the observations in addition to the energy distribution of the
stars. In this paper, we assumed the energy distribution of the stars because
they can be obtained by solving the Fokker-Planck equation in the framework of
loss cone theory. We note that we need to consider the effects of the local
self-gravitational potential of stars and other potentials such as the
galactic potential or the non-dominant SMBH potential in order to handle the
unbound stars more precisely.

Let us qualitatively discuss how the number density varies with the black hole
mass in the star cluster with the assumption that the density profile of the
cusp is $\rho(r)=\rho_{h}(r/r_{\rm h})^{-s}$, where $\rho_{\rm h}$ is the
density at the influence radius. Note that the total stellar mass inside the
influence radius equals to the black hole mass. The velocity dispersion inside
the influence radius follows $\sigma^{2}(r)\approx GM_{\rm BH}/[(1+s)r]$ (Syer
& Ulmer, 1999). Following Frank & Rees (1976), the critical radius, where the
star consumption is balanced with its repopulation, is calculated to be

$r_{\rm crit}=\Bigl{[}\frac{2M^{2}_{\rm BH}r_{\rm
t}N}{15.38(1+s)^{2}\ln(0.11N)M_{\rm c}\rho_{\rm h}r^{s}_{\rm
h}}\Bigr{]}^{\frac{1}{4-s}}\propto M_{\rm BH}^{\frac{17-3s}{6(4-s)}},$ (32)

where we use the scaling procedure of Zhong et al. (2014): $r_{\rm t}\propto
M^{1/3}_{\rm BH}$, $N\propto M_{\rm BH}$, $M_{c}\propto M_{\rm BH}$, $r_{\rm
h}\propto M^{1/2}_{\rm BH}$, and $\rho_{\rm h}\propto M_{\rm BH}/r^{3}_{\rm
h}\propto M^{-1/2}_{\rm BH}$, and we simply adopt $1/2$ for the scaling of
$r_{\rm h}$ instead of $0.54$ which is obtained from the $M_{\rm BH}-\sigma$
relation (Schulze & Gebhardt, 2011). Since the ratio of critical to influence
radii, $r_{\rm crit}/r_{\rm h}$, is proportional to $M_{\rm BH}^{5/[6(4-s)]}$,
it is smaller as the black hole is smaller as long as $s<4$. In other words,
it is expected that the number of the bound stars, which can cause TDEs,
increases as the black hole mass decreases. Figure 8 depicts $r_{\rm
crit}/r_{\rm h}$ as a function of $M_{\rm BH}$ (see also equation 32) in the
range of $10^{3}\,M_{\odot}\leq{M_{\rm BH}}\leq 10^{8}\,M_{\odot}$, where we
adopt $s=1.75$ for the Bahcall-Wolf cusp (Bahcall & Wolf, 1976) and $s=1$ for
the cusp obtained from the $N$-body simulations.

Assuming that the semi-major axis $a_{\rm min}$ corresponding to $E_{\rm min}$
is given by a fixed fraction of $r_{\rm crit}$ as $a_{\rm min}=fr_{\rm crit}$,
where $f$ is a parameter determined by $N$-body simulations, we have $E_{\rm
min}\propto-M_{\rm BH}^{\frac{7-3s}{6(4-s)}}$. When $s<7/3$, $E_{\rm min}$
decreases as $M_{\rm BH}$ decreases. Substituting $E_{\rm min}=-GM_{\rm
BH}/(2fr_{\rm crit})$ into equation (28), we obtain

$e_{\rm min}=1-\frac{1}{\beta}\frac{1}{f}\frac{r_{\rm t}}{r_{\rm crit}}.$ (33)

Since the ratio of $r_{\rm t}$ to $r_{\rm crit}$ is less than $3\times
10^{-4}$ over the whole range of the black hole mass, $e_{\rm min}$ is always
larger $0.994$ for $f=0.05$ and $\beta=1$. Because of $1-e_{\rm
min}\propto{r_{\rm t}/r_{\rm crit}}\propto{M_{\rm BH}^{(s-9)/[6(4-s)]}}$,
$e_{\rm min}$ is closer to 1 for $s<4$ as the black hole mass is larger. We
also confirm that the black hole mass dependence of the pericenter radius is
the same as that of the tidal disruption radius, i.e., $r_{p,\rm{min}}=a_{\rm
min}(1-e_{\rm min})\propto{M_{\rm BH}^{1/3}}\propto{r_{\rm t}}$.

Next, let us see how $e_{\rm min}$ and $e_{\rm lcb}$ (see equations 30 and 33)
depends on the black hole mass the $e$-$\beta$ plane. The left panel of Figure
9 shows it for the $M_{\rm BH}=10^{3}$ and $10^{4}M_{\odot}$ cases. The
$e_{\rm lcb}$ curve of $M_{\rm BH}=10^{3}M_{\odot}$ gets larger than $M_{\rm
BH}=10^{4}M_{\odot}$ case (thick curve). In addition, as mentioned in Section
3.4 (see also Figure 7), most of the bound stars are distributed around the
$e_{\rm lcb}$ curve on the plane. These suggest that the star is distributed
closer to $e=1$ over the whole range of $\beta$ as the black hole mass
increases. The right panel of Fig. 9 depicts how the $e_{\rm min}$ curve
depends on the black hole mass in the $e$-$\beta$ plane. It is clear from the
panel that $e_{\rm min}$ gets closer to 1 as the black hole mass is larger.
This is consistent with that $e_{\rm min}\propto 1-M_{\rm
BH}^{(s-9)/[6(4-s)]}\,(s<4)$ as we estimated in the previous paragraph. In
summary, Figure 9 suggests that the stars are supplied into the black hole on
extremely marginally eccentric to parabolic orbits for a spherical cluster
with $M_{\rm BH}>10^{7}M_{\odot}$ black hole.

In our $N$-body simulations, we have adopted a fixed position and mass of the
central black hole. For a star cluster with an IMBH, however, the Brownian
motion of the black hole modifies the energy distribution of stars so that the
number density and $e$-$\beta$ distributions can be significantly affected.
This is because $M_{\rm BH}/m_{*}$ of the IMBH system is not so large enough
to suppress the Brownian motion. In a forthcoming paper, we will examine how
the Brownian motion affects the number density and $e$-$\beta$ distributions.

Our analysis of the $e$-$\beta$ distributions could provide a good tool to
probe the dynamical status of the stars causing TDEs in a star cluster.
Hayasaki et al. (2013, 2016) and Bonnerot et al. (2016) studied the accretion
disk formation by performing the hydrodynamic simulations, where the authors
have adopted $(\beta,e)=(5,0.8)$ as initial values, although they also used
the other combinations of $(e,\beta)$. This parameter set is clearly ruled out
in our model. However, this does not mean that such a very tightly bound TDE
cannot occur. The tightly bound stars are likely to supply to the loss cone
not by two-body encounters but by the other mechanisms: the tidal separation
of stellar or compact binaries approaching the SMBH (Fragione & Sari, 2018),
accretion disk mediated TDEs (Kennedy et al., 2016), TDEs produced by a
recoiling SMBH (Gualandris & Merritt, 2008; Li et al., 2012) or by a merging
SMBH binary (Hayasaki & Loeb, 2016; Li et al., 2017).

Figure 8:  Black hole mass dependence of $r_{\rm crit}/r_{\rm h}$. The two
slopes for the density profiles: $s=1.75$ (Bahcall-Wolf cusp) and $s=1$
($N$-body simulations) are adopted.

Figure 9:  Black hole mass dependence of two characteristic eccentricities:
$e_{\rm min}$ and $e_{\rm lcb}$ in the $e$-$\beta$ plane. We adopt the
Bahcall-Wolf density cusp for all the models. The left panel shows
$\beta$-dependence of $e_{\rm min}$ and $e_{\rm lcb}$ for the $M_{\rm
BH}=10^{3}M_{\odot}$ and $10^{4}M_{\odot}$ cases. The dashed lines denote the
eccentricity, $e_{\rm lcb}$, between the empty and full loss cone regimes (see
equation 30), whereas the solid lines denote the eccentricity, $e_{\rm min}$,
of most tightly bound stars (see equation 33). The different color shows the
different black hole mass. In the right panel, all the lines represent $e_{\rm
min}$. We adopt the different line style for the different black hole mass for
$10^{3}\,M_{\odot}\leq{M}_{\rm BH}\leq 10^{8}{M}_{\odot}$.

##  5 Conclusion

Understanding tidal disruption events (TDEs) and their light curves provide a
clue to the properties of the central supermassive black hole (SMBH), the
accretion disk around it, and to the stellar density and kinematic
distributions in the nuclear star cluster surrounding the SMBH. The link
between TDEs and central star clusters in galactic nuclei has been a classical
subject of seminal papers, such as Frank & Rees (1976); Bahcall & Wolf (1976)
finding the classical density distribution near central SMBHs; Dokuchaev &
Ozernoi (1977a, b) first noted that accretion and tidal disruption of stars
with low angular momentum causes the density profile to flatten out towards
the black hole (the energy distribution function $f(E)$ drops towards zero as
they showed, today we would call this the empty loss cone region, see also
Ozernoi & Reinhardt (1978) for a summary of the topic at the time). CK78 put
this on a more quantitative footing by using the technique of solving the
orbit-averaged Fokker-Planck equation. Among the classical work in this field
also Rees’s conjecture about the fate of tidal debris (Rees, 1988) is most
noteworthy.

Our study generalizes and expands this by computing approximate distributions
of bound and unbound stars in energy-angular momentum space and by
transforming the results into density and velocity dispersion profiles. At the
same time, we extend the classical Rees’ picture by considering different
orbital eccentricities into the tidal disruption radius, which may result - in
extreme cases - in 100% mass loss (i.e., hyperbolic TDEs) or 100% accretion of
the tidal debris (i.e., eccentric TDEs). Here we follow the work of Hayasaki
et al. (2013, 2018) and Park & Hayasaki (2020). By their works, we can improve
the understanding of the link between light curves and other TDE
characteristics and the stellar distribution of stars in the surrounding
nuclear star clusters.

Amaro-Seoane & Spurzem (2001) and Amaro-Seoane et al. (2004) proposed an
another model for the loss cone in a spherical star cluster. It is interesting
to note that they have also derived the density and velocity dispersion of
bound and unbound loss-cone stars by using moment equations of the basic
Fokker-Planck equation and very similar principles to CK78 and this paper. Due
to the use of moment equations (so-called gaseous model of star clusters, see
e.g. Giersz & Spurzem 1994b) their analysis is completely based on density and
velocity profiles rather than orbits with energy and angular momentum. Their
model takes into account an anisotropic velocity distribution also for the
unbound stars. In the future, a more quantitative comparison of the two models
could be done.

We have tested our theoretical models by performing direct $N$-body
simulations; comparison with our analytical model helps to understand the
scaling behavior of the $N$-body simulations since we can still not yet do
realistic particle numbers for them. Our primary conclusions are summarized as
follows:

  1. 1.

The number density of the bound stars has been derived as a function of an
orbital eccentricity and energy $n(e,E)$ and of a penetration factor and
energy $n(\beta,E)$. These two-number densities are given by equations (13)
and (19), respectively. For the unbound stars, we have found that
$n(\beta,E)\propto 1/\beta^{2}$ and $n(e,E)$ is a function of only energy,
which are given by equations (20) and (22). We have confirmed that the energy
integrals of the number densities of both the bound and unbound stars are good
agreement with those obtained by performing $N$-body simulations.

  2. 2.

We have also analytically derived the three characteristic orbital
eccentricities: $e_{\rm min}$, $e_{\rm max}$, and $e_{\rm lcb}$ in the loss
cone region, where $e_{\rm min}$ and $e_{\rm max}$ take the minimum and
maximum values for a given $\beta$, respectively, whereas $e_{\rm lcb}$
represents the orbital eccentricity which gives the boundary between the empty
and full loss cone regimes. These eccentricities are given by equations (28),
(29), and (30), respectively. We have confirmed that the stars causing TDEs
are distributed between $e_{\rm min}$ and $e_{\rm max}$ on the $e-\beta$ plane
by N-body experiments. Moreover, we find most of the bound stars are focused
between $e_{\rm lcb}$ and $e=1$, i.e., in the full loss cone regime, whereas
the remaining bound stars are originating from the empty loss cone regime.
This result is consistent with the loss cone theory.

## Acknowledgements

S.Z. has been supported by the National Natural Science Foundation of China
(NSFC 11603067) and acknowledges the support from Yunnan Astronomical
Observatories, Chinese Academic of Sciences. K.H. has been supported by the
Korea Astronomy and Space Science Institute (KASI) under the R&D program
supervised by the Ministry of Science, ICT and Future Planning, and by the
Basic Science Research Program through the National Research Foundation of
Korea (NRF) funded by the Ministry of Education (2016R1A5A1013277,
2017R1D1A1B03028580, and 2020R1A2C1007219 (K.H.)). K.H. has been also
supported by the National Supercomputing Center with supercomputing resources
including technical support (KSC-2019-CRE-0082 (K.H.)). The authors
acknowledge the Yukawa Institute for Theoretical Physics (YITP) at Kyoto
University. Discussions during the YITP workshop YITP-T-19-07 on International
Molecule-type Workshop ”Tidal Disruption Events: General Relativistic
Transients” were useful to complete this work. The authors also acknowledge
support by the Chinese Academy of Sciences (CAS) through the Silk Road Project
at NAOC. We are grateful for the support from the Sino-German Center
(DFG/NSFC) under grant no. GZ1289. S.L., P.B. and R.S. acknowledge the
Strategic Priority Research Program (Pilot B) Multi-wavelength gravitational
wave universe of the Chinese Academy of Sciences (No. XDB23040100). P.B.
acknowledges the President International Fellowship for Visiting Scientists
program of CAS. P.B. and R.S. acknowledge the support of the Volkswagen
Foundation under the Trilateral Partnerships grant No. 97778. 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.

## References

  * Amaro-Seoane et al. (2004) Amaro-Seoane, P., Freitag, M., & Spurzem, R. 2004, MNRAS, 352, 655
  * Amaro-Seoane & Spurzem (2001) Amaro-Seoane, P. & Spurzem, R. 2001, MNRAS, 327, 995
  * Auchettl et al. (2017) Auchettl, K., Guillochon, J., & Ramirez-Ruiz, E. 2017, ApJ, 838, 149
  * Bahcall & Wolf (1976) Bahcall, J. N. & Wolf, R. A. 1976, ApJ, 209, 214
  * Bonnerot et al. (2016) Bonnerot, C., Rossi, E. M., Lodato, G., & Price, D. J. 2016, MNRAS, 455, 2253
  * Cohn & Kulsrud (1978) Cohn, H. & Kulsrud, R. M. 1978, ApJ, 226, 1087
  * Dokuchaev & Ozernoi (1977a) Dokuchaev, V. I. & Ozernoi, L. M. 1977a, Zhurnal Eksperimentalnoi i Teoreticheskoi Fiziki, 73, 1587
  * Dokuchaev & Ozernoi (1977b) —. 1977b, Pisma v Astronomicheskii Zhurnal, 3, 295
  * Evans & Kochanek (1989) Evans, C. R. & Kochanek, C. S. 1989, ApJ, 346, L13
  * Fragione & Sari (2018) Fragione, G. & Sari, R. 2018, ApJ, 852, 51
  * Frank & Rees (1976) Frank, J. & Rees, M. J. 1976, MNRAS, 176, 633
  * Gezari et al. (2015) Gezari, S., Chornock, R., Lawrence, A., Rest, A., Jones, D. O., Berger, E., Challis, P. M., & Narayan, G. 2015, ApJ, 815, L5
  * Gezari et al. (2012) Gezari, S., Chornock, R., Rest, A., Huber, M. E., Forster, K., Berger, E., Challis, P. J., Neill, J. D., Martin, D. C., Heckman, T., Lawrence, A., Norman, C., Narayan, G., Foley, R. J., Marion, G. H., Scolnic, D., Chomiuk, L., Soderberg, A., Smith, K., Kirshner, R. P., Riess, A. G., Smartt, S. J., Stubbs, C. W., Tonry, J. L., Wood-Vasey, W. M., Burgett, W. S., Chambers, K. C., Grav, T., Heasley, J. N., Kaiser, N., Kudritzki, R.-P., Magnier, E. A., Morgan, J. S., & Price, P. A. 2012, Nature, 485, 217
  * Giersz & Spurzem (1994a) Giersz, M. & Spurzem, R. 1994a, MNRAS, 269, 241
  * Giersz & Spurzem (1994b) —. 1994b, MNRAS, 269, 241
  * Gualandris & Merritt (2008) Gualandris, A. & Merritt, D. 2008, ApJ, 678, 780
  * Guillochon & Ramirez-Ruiz (2013) Guillochon, J. & Ramirez-Ruiz, E. 2013, ApJ, 767, 25
  * Hayasaki & Loeb (2016) Hayasaki, K. & Loeb, A. 2016, Scientific Reports, 6, 35629
  * Hayasaki et al. (2013) Hayasaki, K., Stone, N., & Loeb, A. 2013, MNRAS, 434, 909
  * Hayasaki et al. (2016) —. 2016, MNRAS, 461, 3760
  * Hayasaki et al. (2018) Hayasaki, K., Zhong, S., Li, S., Berczik, P., & Spurzem, R. 2018, ApJ, 855, 129
  * Heggie (2014a) Heggie, D. C. 2014a, ArXiv e-prints
  * Heggie (2014b) —. 2014b, MNRAS, 445, 3435
  * Hills (1975) Hills, J. G. 1975, Nature, 254, 295
  * Holoien et al. (2016) Holoien, T. W. S., Kochanek, C. S., Prieto, J. L., Stanek, K. Z., Dong, S., Shappee, B. J., Grupe, D., Brown, J. S., Basu, U., Beacom, J. F., Bersier, D., Brimacombe, J., Danilet, A. B., Falco, E., Guo, Z., Jose, J., Herczeg, G. J., Long, F., Pojmanski, G., Simonian, G. V., Szczygieł, D. M., Thompson, T. A., Thorstensen, J. R., Wagner, R. M., & Woźniak, P. R. 2016, MNRAS, 455, 2918
  * Holoien et al. (2014) Holoien, T. W.-S., Prieto, J. L., Bersier, D., Kochanek, C. S., Stanek, K. Z., Shappee, B. J., Grupe, D., Basu, U., Beacom, J. F., Brimacombe, J., Brown, J. S., Davis, A. B., Jencson, J., Pojmanski, G., & Szczygieł, D. M. 2014, MNRAS, 445, 3263
  * Hung et al. (2017) Hung, T., Gezari, S., Blagorodnova, N., Roth, N., Cenko, S. B., Kulkarni, S. R., Horesh, A., Arcavi, I., McCully, C., Yan, L., Lunnan, R., Fremling, C., Cao, Y., Nugent, P. E., & Wozniak, P. 2017, ApJ, 842, 29
  * Kennedy et al. (2016) Kennedy, G. F., Meiron, Y., Shukirgaliyev, B., Panamarev, T., Berczik, P., Just, A., & Spurzem, R. 2016, MNRAS, 460, 240
  * Li et al. (2012) Li, S., Liu, F. K., Berczik, P., Chen, X., & Spurzem, R. 2012, ApJ, 748, 65
  * Li et al. (2017) Li, S., Liu, F. K., Berczik, P., & Spurzem, R. 2017, ApJ, 834, 195
  * Lodato et al. (2009) Lodato, G., King, A. R., & Pringle, J. E. 2009, MNRAS, 392, 332
  * Lodato & Rossi (2011) Lodato, G. & Rossi, E. M. 2011, MNRAS, 410, 359
  * MacLeod et al. (2012) MacLeod, M., Guillochon, J., & Ramirez-Ruiz, E. 2012, ApJ, 757, 134
  * Magorrian & Tremaine (1999) Magorrian, J. & Tremaine, S. 1999, MNRAS, 309, 447
  * Merritt (2013) Merritt, D. 2013, Classical and Quantum Gravity, 30, 244005
  * Miller et al. (2015) Miller, J. M., Kaastra, J. S., Miller, M. C., Reynolds, M. T., Brown, G., Cenko, S. B., Drake, J. J., Gezari, S., Guillochon, J., Gultekin, K., Irwin, J., Levan, A., Maitra, D., Maksym, W. P., Mushotzky, R., O’Brien, P., Paerels, F., de Plaa, J., Ramirez-Ruiz, E., Strohmayer, T., & Tanvir, N. 2015, Nature, 526, 542
  * Mummery & Balbus (2020) Mummery, A. & Balbus, S. A. 2020, MNRAS, 492, 5655
  * Ozernoi & Reinhardt (1978) Ozernoi, L. M. & Reinhardt, M. 1978, Ap&SS, 59, 171
  * Park & Hayasaki (2020) Park, G. & Hayasaki, K. 2020, ApJ, 900, 3
  * Phinney (1989) Phinney, E. S. 1989, in IAU Symposium, Vol. 136, The Center of the Galaxy, ed. M. Morris, 543
  * Rees (1988) Rees, M. J. 1988, Nature, 333, 523
  * Schulze & Gebhardt (2011) Schulze, A. & Gebhardt, K. 2011, ApJ, 729, 21
  * Shapiro & Marchant (1978) Shapiro, S. L. & Marchant, A. B. 1978, ApJ, 225, 603
  * Spitzer (1987) Spitzer, L. 1987, Dynamical evolution of globular clusters
  * Stone & Metzger (2016) Stone, N. C. & Metzger, B. D. 2016, MNRAS, 455, 859
  * Strubbe & Quataert (2009) Strubbe, L. E. & Quataert, E. 2009, MNRAS, 400, 2070
  * Syer & Ulmer (1999) Syer, D. & Ulmer, A. 1999, MNRAS, 306, 35
  * Touma & Tremaine (1997) Touma, J. & Tremaine, S. 1997, MNRAS, 292, 905
  * van Velzen et al. (2019) van Velzen, S., Stone, N. C., Metzger, B. D., Gezari, S., Brown, T. M., & Fruchter, A. S. 2019, ApJ, 878, 82
  * Wang & Merritt (2004) Wang, J. & Merritt, D. 2004, ApJ, 600, 149
  * Zhong et al. (2014) Zhong, S., Berczik, P., & Spurzem, R. 2014, ApJ, 792, 137
  * Zhong et al. (2015) —. 2015, ApJ, 811, 22