# Crime and COVID-19 in Rio de Janeiro: How does organized crime shape the
disease evolution?

Nuno Crokidakis and Lucas Sigaud  Instituto de Física, Universidade Federal
Fluminense
Niterói, Rio de Janeiro, Brazil
<EMAIL_ADDRESS>
(Day Month Year; Day Month Year)

###### Abstract

The city of Rio de Janeiro is one of the biggest cities in Brazil. Drug gangs
and paramilitary groups called milícias control some regions of the city where
the government is not present, specially in the slums. Due to the
characteristics of such two distinct groups, it was observed that the
evolution of COVID-19 is different in those two regions, in comparison with
the regions controlled by the government. In order to understand qualitatively
those observations, we divided the city in three regions controlled by the
government, by the drug gangs and by the milícias, respectively, and we
consider a SIRD-like epidemic model where the three regions are coupled.
Considering different levels of exposure, the model is capable to reproduce
qualitatively the distinct evolution of the COVID-19 disease in the three
regions, suggesting that the organized crime shapes the COVID-19 evolution in
the city of Rio de Janeiro. This case study suggests that the model can be
used in general for any metropolitan region with groups of people that can be
categorized by their level of exposure.

###### keywords:

Dynamics of social systems, Compartmental model, Collective phenomena

PACS Nos.: 87.23.Ge, 89.20.-a, 89.65.-s, 89.75.Fb

##  1 Introduction

The evolution of epidemics is one of the most dangerous problems for a society
and its spreading mechanisms have been the source of interest to the exact
sciences’ scientific community for over a century [1]. The humanity faced
severe pandemics during its evolution, like the Spanish flu in 1917, the Honk
Kong flu (H3N2) of 1968 and the swine flu (H1N1) in 2009. Several efforts were
done since the 70’s in order to understand the mathematical evolution and
spreading of diseases [2, 3].

The recent outbreak of the COVID-19 pandemic and its effect on society in
different countries and its many diverse aspects are now the object of a long
list of different studies on the subject of epidemic modelling. An immense and
not always reliable amount of data has been produced over the past year
regarding infection rates, contagion levels, population-exposure-related
contamination, age-related fatality rates, statistical spreading, quarantine
efficiency, death tolls and mortality rates by population; these have been
exhaustively given worldwide coverage by the press, health organizations and
data banks [4, 5] and we will not delve further here in details regarding
these parameters.

For many different diseases, compartmental models focusing on a population
divided into Susceptible (S), Infected (I) and Recovered (R) individuals,
namely SIR, have been used extensively. It often employs mean-field
approximation for the individuals network which, albeit unrealistic, was able
to model many different aspects of epidemic behaviour, including the necessity
of artificial immunity through vaccination in order to contain the spreading
of the disease [6]. Taking into account the pre-symptomatic, symptomatic and
asymptomatic possibilities for COVID-19 infection and spreading, SIR models
fail because of two important issues: i) the possibility of pre-symptomatic
and asymptomatic contagion through individuals which are not known as patients
(i.e. belonging to the I compartment) led to social isolation and quarantine
measures throughout the world; and, in consequence, ii) different levels of
social isolation among groups coexisting in the same geographical area make
the mean-field approximation inadequate for COVID-19 contagion behaviour [7].

In order to tackle the first of these two issues, many groups proposed
different compartmental modelling of the COVID-19 pandemic spreading using
SIQR, with the introduction of a fourth class of individuals, namely the
Quarantined ones (Q) [7, 8, 9, 10, 11, 12]. The introduction of the Q
population is nothing new, however, and quarantine effects on epidemic models
have been studied in this manner for a few decades [13]. Nonetheless, SIQR
models have been useful in modelling the COVID-19 spreading [9] and, for
example, predicting phenomena such as the so-called second wave incidences
[7]. Yah and Li have taken a further step and included age as a novel
parameter that influences the connections between the different compartments
in a SIQR model [14]. El Fatini et al. [10] have included vaccination
strategies employed in previous models [6] in the very early stages of the
COVID-19 onset, while Hackbart [11] has proposed models including mass testing
and selected quarantined individuals based on these testing parameters, in
order to target infected individuals to reduce contagion rise.

On the other hand, the second issue listed above needs further addressing.
Odagaki included the fluctuation of social isolation policies in Japan and how
it was obeyed by citizens in order to explain the crescent three-wave pattern
that has arisen there over the course of the pandemic [8], but still kept the
compartmental model with mean-field approximation. Weisbuch attempts to
circumvent the limitations of the mean-field approximation by including a
mobile initial population at both city and countryside in order to simulate
urban exodus and its influence on the COVID-19 spreading [15], and then
including crossover parameters linking interchange between city and
countryside populations. Although it is an oversimplification of the dynamics
of evolving networks [16], it is an efficient move towards including different
types of population in the same model. Kraemer et al. have studied the
consequences of human mobility and social isolation policies in China during
the COVID-19 first months and found strong correlation patterns [17].

In this work, we aim to address the question of different patterns of social
isolation and contact in a population confined to a single geographical
location, in order to model groups of individuals either with different
degrees of compliance to the same social isolation policy or subjected to
different social isolation policies. In order to have numerical data to which
to compare the proposed model, we take the particular situation of Rio de
Janeiro as a case study, which has been the object of previous works [18]. We
chose Rio de Janeiro city because of its very particular subdivision into
different areas, as is more thoroughly explained elsewhere [19, 20, 21]. In
short, although in principle the entire city is subjected to the regular
official municipal, state and federal policies, two other groups, with
divergent governance systems, retain political and enforcement control over
many different sub-regions of Rio de Janeiro’s metropolitan area, namely drug
trafficking gangs and the paramilitary groups called milícias. While drug
trafficking gangs exert stay-at-home regulations in their controlled areas,
since their financial activities rely on underground commerce and minimal
governmental interference, milícias, operated by police and ex-police officers
as a parallel power, force local populations to not stay at home, in order to
maintain their commercial tax and protection monetary charges. This behavior
in different slums has been monitored using, for example, credit card and
mobile cellphones usage data [22] and an illustrative map of the different
regions is shown by Bruce et al. [21], but it is not reproduced here because
we do not take the geographical subdivisions into account in our model.
Therefore, in the same metropolitan area there is an interchange of
individuals subdivided into three different but connected populations, each
with its own exposure parameter, which is related to the level of social
distancing that the individuals for each group adopt. In other words, the
exposure coefficients are related to the degree of exposure that each
individual both is subjected to other individuals and subject other
individuals to him/herself. We aim to verify if the distinct contagion rates
observed by the authors [21] are a direct consequence of different exposure
patterns.

In order to try to accomplish this, we refrain from using the SIQR model,
since, in general, our goal is to address the effects of exposure patterns on
the spread of COVID-19, and, in particular for our case study, quarantined
individuals data are very hard to obtain in both drug trafficking gangs and
milícias controlled regions. We opt instead for the SIRD model, which includes
a compartment for Dead (D) individuals and is very useful when death rates
data are more readily available, which is the case for COVID-19 both in
general [23, 24] and in our particular case study, where death rates due to
COVID-19 are reported to be higher for milícia-controlled areas and lower for
drug-trafficking-gang-controlled areas, with respect to the rates observed for
regions under the direct influence of neither type [21].

##  2 Model

Let us recall the Susceptible-Infected-Recovered-Dead (SIRD) model [1]. The
population is divided in four compartments, and the equations for the
evolution of such compartments are

$\displaystyle\frac{dS}{dt}$ $\displaystyle=$ $\displaystyle-\beta\,S\,I$ (1)
$\displaystyle\frac{dI}{dt}$ $\displaystyle=$
$\displaystyle\beta\,S\,I-\gamma\,I-\delta\,I$ (2)
$\displaystyle\frac{dR}{dt}$ $\displaystyle=$ $\displaystyle\gamma\,I$ (3)
$\displaystyle\frac{dD}{dt}$ $\displaystyle=$ $\displaystyle\delta\,I$ (4)

where $S,I,R$ and $D$ denote the fractions of Susceptible, Infected, Recovered
and Dead individuals, respectively. The parameters $\beta$, $\gamma$ and
$\delta$ denote the infection, recovery and death probabilities, respectively.

As discussed in the Introduction, we are interested to verify if the
hypothesis pointed by the authors in Ref. [21] is plausible. In other words,
if the restrictions to stay at home imposed by the drug trafficking gangs in
some regions of the city and the pressure to not stay at home in areas
dominated by the milícias influenced the COVID-19 spreading in distinct areas
of Rio de Janeiro. For this purpose, we divided the Rio de Janeiro city in
three regions, namely:

  * •

1: areas controlled by the government

  * •

2: areas controlled by drug trafficking gangs

  * •

3: areas controlled by milícias

In each of such regions, we consider a Susceptible-Infected-Recovered-Dead
(SIRD) model. In addition, we consider in the models’ equations the exposure
restrictions (in the case of drug gangs’ areas) or the absence of such
restrictions (in the case of milícias). In such case, we consider three
additional parameters to the SIRD model: the coefficients related to the
population’s exposure in each region, namely $C_{1},C_{2}$ and $C_{3}$ for
regions $1,2$ and $3$, respectively. Those coefficients take into account the
couplings of the population of a given region with the population of the other
regions.

Thus, the model is defined through the equations

$\displaystyle\frac{dS_{u}}{dt}$ $\displaystyle=$
$\displaystyle-\beta\,S_{u}\,C_{u}\sum_{v=1}^{3}C_{v}\,I_{v}$ (5)
$\displaystyle\frac{dI_{u}}{dt}$ $\displaystyle=$
$\displaystyle\beta\,S_{u}\,C_{u}\sum_{v=1}^{3}C_{v}\,I_{v}-\gamma\,I_{u}-\delta\,I_{u}$
(6) $\displaystyle\frac{dR_{u}}{dt}$ $\displaystyle=$
$\displaystyle\gamma\,I_{u}$ (7) $\displaystyle\frac{dD_{u}}{dt}$
$\displaystyle=$ $\displaystyle\delta\,I_{u}$ (8)

where $S_{u},I_{u},R_{u}$ and $D_{u}$ denote the fractions of Susceptible,
Infected, Recovered and Dead individuals, respectively, in each region $u$,
with $u=\\{1,2,3\\}$. The terms $S_{u}\,C_{u}\,C_{v}\,I_{v}$ in Eqs. (5) and
(6) introduce the above-mentioned coupling among the regions. In other words,
the coupling terms correspond to the infection of individuals in a given
location $u$ by infected individuals [15, 25] of all the three locations
($v=\\{1,2,3\\}$) that move along the city. Notice that we do not restrict the
$C_{u}$ values to the range $0<C_{u}<1$ since, strictly speaking, they are not
probabilities – they are related to the degree of exposure of individuals and
they depend on the specific area of the city such individuals live. For
simplicity, we considered homogeneous probabilities $\beta$, $\gamma$ and
$\delta$ for the three regions.

##  3 Results

We considered the initial conditions $R_{u}(0)=D_{u}(0)=0$ and
$I_{u}(0)=0.001$ for $u=\\{1,2,3\\}$. The normalization condition is obeyed at
each time step, i.e., $S(t)=1-I(t)-R(t)-D(t)$, where
$S(t)=S_{1}(t)+S_{2}(t)+S_{3}(t),I(t)=I_{1}(t)+I_{2}(t)+I_{3}(t),R(t)=R_{1}(t)+R_{2}(t)+R_{3}(t)$
and $D(t)=D_{1}(t)+D_{2}(t)+D_{3}(t)$. In this case, the initial condition for
the susceptible populations is $S_{1}(0)=S_{2}(0)=S_{3}(0)=S(0)/3$.

Since we have distinct behaviors regarding the social isolation in the three
regions, we have $C_{2}<C_{1}<C_{3}$, i.e., there is less exposure in drug
trafficking regions (area 2) and more in milícia controlled ones (area 3). As
we are interested in the qualitative description of the results obtained in
[21], and the parameters $C_{1},C_{2}$ and $C_{3}$ are the novelty of the
model (as well as the division in three regions), we fixed the transition
probabilities $\beta=0.32$, $\gamma=0.045$ and $\delta=0.054$, based on
estimates for COVID-19 parameters [18, 26].

Figure 1: Fractions of Susceptible, Infected, Recovered and Dead individuals
for the areas 1 (control - black curve), 2 (drug gangs - red curve) and 3
(milícias - blue curve) as functions of time for $C_{1}=1.0$, $C_{2}=0.4$ and
$C_{3}=1.3$. We can see distinct time evolutions for the 3 regions due to
distinct individuals’ degrees of exposure.

In Fig. 1 we exhibit the fractions of Susceptible, Infected, Recovered and
Dead individuals for the three areas 1, 2 and 3 as functions of time for
$C_{1}=1.0$, $C_{2}=0.4$ and $C_{3}=1.3$, obtained from the numerical
integration of Eqs. (5) - (8). Since the exposure coefficients are arbitrary
parameters, we chose as unity the one for government controlled areas and
modified the other two accordingly. We can see that if we consider the
coupling among the 3 regions, together with the distinct exposure coefficients
$C_{u}$ ($u=\\{1,2,3\\}$), the model leads to distinct evolutions in each
region $u$. In particular, one can see that the outbreak for the infections is
more pronounced in the region 3 dominated by the milícias. On the other hand,
one observe the flattening of the curve of infections in region 2, dominated
by the drug gangs. Such difference among the three regions leads to distinct
evolutions of the number of the other three compartments $S,R$ and $D$.
Specifically, the number of deaths is higher in region 3 and lower in region
2, in comparison with the region 1 controlled by the State. These results
suggest that the scenario observed in [21], i.e. distinct evolutions of the
infections and deaths in the three regions of Rio de Janeiro, is indeed due to
the distinct exposure patterns imposed by drug trafficking gangs and milícias,
in such two regions (2 and 3) that are not controlled by the State.

Figure 2: Comparative fractions of Infected individuals for the areas 1
(control), 2 (drug gangs) and 3 (milícias) as functions of time for fixed
$C_{1}=1.0$. Left panel: fixed $C_{3}=1.3$ with $C_{2}=0.4$ (full lines) and
$C_{2}=0.2$ (dashed lines). Right panel: fixed $C_{2}=0.4$ with $C_{3}=1.3$
(full lines) and $C_{3}=1.8$ (dashed lines).

In Fig. 2 we exhibit comparative results for the fraction of Infections in the
three regions. We fixed the standard exposure coefficient for region 1 as
$C_{1}=1.0$, and present two distinct comparisons. First, we consider the
curves for fixed $C_{3}=1.3$ (as in Fig. 1) and distinct values of $C_{2}$
(Fig. 2, left-hand side), namely $C_{2}=0.4$ (full lines) and $C_{2}=0.2$
(dashed lines). One can see that lowering the exposure coefficient for the
more social isolated group (the one controlled by drug trafficking gangs in
our case study) considerably lowers the infection rates for individuals
belonging to this group, while maintaining the other two groups fairly
unchanged. On the other hand, in the right-hand side panel we consider the
curves for fixed $C_{2}=0.4$ (as in Fig. 1) and distinct values of $C_{3}$,
namely $C_{3}=1.3$ (full lines) and $C_{3}=1.8$ (dashed lines). It is very
interesting to see that in this case, when we change the degree of exposure of
the less social isolated group (the milícia controlled one, in our case
study), the infected population of all three groups increase considerably, and
roughly by the same proportion. This clearly indicates that, in a city divided
in three different groups, with different degrees of exposure, increasing the
social distancing among individuals that already have a lower degree of
exposure, albeit lowering the infection rates of this particular group, does
not change the overall scenario: the two other groups, which were already
responsible for the larger share of infected individuals are only slightly
affected by it. Meanwhile, if we increase the degree of exposure (i.e.
decrease the social distancing) of individuals which were not respecting
social isolation policies in the first place, the consequences for the entire
population are indeed drastic. It is worthwhile to mention that this does not
refer to our case study alone - in any populated area where one is able to
subdivide the inhabitants into groups with three different degrees of
exposure, this conclusion holds. Moreover, since we assumed a mean-field-like
approximation modulated by the exposure coefficients, we did not take into
account the geographical properties of Rio de Janeiro City and its controlled
regions; therefore, this aspect takes a more general form.

From Fig. 2 one can also see some of the more general infection behaviors that
always appear when looking at epidemic models: the flattening (height decrease
and enlargement) of the infection curves over time for higher social isolation
patterns and a shorter pandemic time-frame with increased infection (and,
consequently, death) rates for lower social distancing. It is meaningful that
these behaviors are reproduced as well in a population subdivided as we
presented here, with these differences showing simultaneously between the
three different groups; i.e. group 2 has flatter infection rate curve over
time than groups 1 and 3. From the changes shown in Fig. 2 for the three
curves when different values for parameters $C_{2}$ or $C_{3}$ were chosen,
one can infer that the existence of a group which does not respect social
distancing measures has a very strong effect on individuals that are not under
a strict social distancing policy.

##  4 Final Remarks

In this work we considered a Susceptible-Infected-Recovered-Dead (SIRD) model
in order to study the evolution of COVID-19 in the Rio de Janeiro city. Taking
into account that the organized crime dominates several neighborhoods of the
city, we divided the city in three areas, namely the areas controlled by the
government (area 1), by the drug trafficking gangs (area 2) and by the
paramilitary groups called milícias (area 3). Taking into account that the
behavior of individuals in each region is distinct, determined by the
government in area 1 and by the organized crime in areas 2 and 3, we
considered distinct exposure coefficients for the three regions. In addition,
we coupled the 3 regions through these exposure coefficients, considering a
not-independent SIRD model in each region.

At this point, some limitations can be discussed. First of all, we considered
a mean-field-like approach, where each individual can interact with all
others. In this case, spatial features were not considered in the model. Our
work does not consider explicitly an upper bound for the capacity of the
healthcare system. Under-reporting is another feature that is not modeled
here. We have also assumed permanent immunity, i.e., after recovering from
COVID-19 individuals cannot be infected anymore. In addition, as we are
interested in explaining a phenomenon observed in the beginning of the
epidemic evolution, we did not implement other containment measures besides
the exposure restrictions imposed by each controlling agent.

It is important to say that our target was not to describe exactly the
evolution of COVID-19 in Rio de Janeiro city, but we are interested in the
qualitative description of an observed phenomenon, namely the distinct
evolution of the number of cases and deaths in the three distinct regions, as
observed by economists in recent works [21, 22]. In such a case, despite the
above-mentioned limitations, the model presented here is capable to explain in
a simple way how the organized crime shaped the COVID-19 evolution in Rio de
Janeiro.

Furthermore, since we employ a mean-field-like approach, the geographical
features and subdivisions of our case study were not taken into account, which
in a certain way can be actually considered a feature instead of a liability:
it shows, in a first-order approximation, the consequences of the epidemic
progression in individuals belonging to different exposure patterns. It can be
inferred from our results that a group of individuals that respects social
distancing policies has less impact on the overall pandemic picture – besides,
of course, lowering their own infection and death risks – than a group of
individuals that presents higher degrees of exposure. The increase of exposure
in the latter group considerably enhances the infection rate, while the
decrease of exposure in the former, social isolated, group has a much lesser
impact on the progression of the disease.

## Acknowledgments

The authors acknowledge financial support from the Brazilian scientific
funding agencies CNPq (Grants 303025/2017-4 and 311019/2017-0) and FAPERJ
(Grant 203.217/2017).

## References

  * [1] W. O. Kermack, A. G. McKendrick, A contribution to the mathematical theory of epidemics, Proceedings of the Royal Society of London, Series A, 115 (772):700–721 (1927).
  * [2] R. M. Anderson, R. M. May, Infectious Diseases of Humans: Dynamics and Control (Oxford University Press, Oxford, 1991).
  * [3] N. T. J. Bailey, The Mathematical Theory of Infectious Diseases and its Application (Hafner Press, New York, 1975).
  * [4] The WHO official COVID-19 website (https://www.who.int/emergencies/diseases/novel-coronavirus-2019).
  * [5] The worldometer data bank for COVID-19 statistics (https://www.worldometers.info/coronavirus/)
  * [6] R.M. Anderson and R.M. May, Directly Transmitted Infectious Diseases: Control by Vaccination, Science 215, 1053 (1982).
  * [7] T. Odagaki, Exact properties of SIQR model for COVID-19, Physica A 564, 125564 (2021).
  * [8] T. Odagaki, Self‑organized wavy infection curve of COVID‑19, Scientific Reports 11, 1936 (2021).
  * [9] N. Crokidakis, Modeling the early evolution of the COVID-19 in Brazil: Results from a Susceptible–Infectious–Quarantined–Recovered (SIQR) model, International Journal of Modern Physics C 31, 2050135 (2020).
  * [10] M. El Fatini, R. Pettersson, I. Sekkak and R. Taki, A stochastic analysis for a triple delayed SIQR epidemic model with vaccination and elimination strategies, Journal of Applied Mathematics and Computing 64, 781 (2020).
  * [11] G. Hackbart, Heterogeneous SIQR Models with Mass Testing and Targeted Quarantine and the Spread of Infectious Diseases, Working paper (July, 2020), Available at SSRN: $https://ssrn.com/abstract=3643816$ or $http://dx.doi.org/10.2139/ssrn.3643816$.
  * [12] N. Crokidakis, S. M. D. Queirós, Probing into the effectiveness of self-isolation policies in epidemic control, Journal of Statistical Mechanics P06003 (2012).
  * [13] H. Hethcote, M. Zhien and L. Shengbing, Effects of quarantine in six endemic models for infectious diseases, Mathematical Biosciences 180, 141 (2002).
  * [14] H. Yah and J. Li, SIQR dynamics in a random network with heterogeneous connections with infection age, Journal of Nonlinear Sciences & Applications 14, 196 (2021).
  * [15] G. Weisbuch, Urban exodus and the dynamics of COVID-19 pandemics, Physica A 569, 125780 (2021).
  * [16] X. Zhang, C.Z. Shan Jin and H. Zhu, Complex dynamics of epidemic models on adaptive networks, Journal of Differential Equations 266, 803 (2019).
  * [17] M. U. G. Kraemer et. al., The effect of human mobility and control measures on the COVID-19 epidemic in China, Science 10.1126/science.abb4218 (2020).
  * [18] N. Crokidakis, COVID-19 spreading in Rio de Janeiro, Brazil: Do the policies of social isolation really work?, Chaos, Solitons & Fractals 136, 109930 (2020).
  * [19] E. D. Arias and C. D. Rodrigues, The myth of personal security: Criminal gangs, dispute resolution, and identity in Rio de Janeiro’s favelas, Latin American Politics And Society 48, 53 (2006).
  * [20] B. Magaloni, E. Franco-Vivanco and V. Melo, Killing in the slums: Social order, criminal governance, and police violence in Rio de Janeiro, American Political Science Review 114, 552 (2020).
  * [21] R. Bruce, A. Cavgias, L. Meloni, Filling the Void? Organized Crime and COVID-19 in Rio De Janeiro, Working paper (August 21, 2020), Available at SSRN: $https://ssrn.com/abstract=3678840$ or $http://dx.doi.org/10.2139/ssrn.3678840$.
  * [22] L. Brotherhood, T. Cavalcanti, D. Da Mata, C. Santos, Slums and Pandemics (August 2, 2020), Available at SSRN: $https://ssrn.com/abstract=3665695$ or $http://dx.doi.org/10.2139/ssrn.3665695$
  * [23] R. M. da Silva, Using the SIRD model to characterize the COVID-19 spreading in the states of Paraná, Rio Grande do Sul, and Santa Catarina, SciELO Preprints, $DOI:https://doi.org/10.1590/SciELOPreprints.764$
  * [24] S. Chatterjee, A. Sarkar, S. Chatterjee, M. Karmakar, R. Paul, Studying the progress of COVID-19 outbreak in India using SIRD model, Indianian Journal of Physics (2020), $https://doi.org/10.1007/s12648-020-01766-8$
  * [25] M. A. Pires et. al., What is the potential for a second peak in the evolution of SARS-CoV-2 in emerging and developing economies? Insights from a SIRASD model considering the informal economy, arXiv:2005.09019 (2020).
  * [26] Fei Zu et. al., Clinical course and risk factors for mortality of adult inpatients with COVID-19 in Wuhan, China: a retrospective cohort study, Lancet 395, Issue 10229, 1054-1062 (2020).


# The influence of the secular perturbation of an intermediate-mass companion:
II. Ejection of hypervelocity stars from the Galactic Center

Xiaochen Zheng （郑晓晨） Department of Astronomy, Tsinghua University, Beijing,
China  Douglas N. C. Lin (林潮) Department of Astronomy and Astrophysics,
University of California, Santa Cruz, USA  Institute for Advanced Studies,
Tsinghua University, Beijing, China  Shude Mao (毛淑德) Department of Astronomy,
Tsinghua University, Beijing, China  National Astronomical Observatories of
China, Chinese Academy of Sciences, 20A Datun Road, Beijing, China

(Received Sep 1, 2020; Revised Apr 2, 2021; Accepted Apr 6, 2021)

###### Abstract

There is a population of stars with velocities in excess of 500 km s-1
relative to the Galactic center. Many, perhaps most, of these hyper-velocity
stars (HVSs) are B stars, similar to the disk and S stars in a nuclear cluster
around a super- massive black hole (SMBH) near $\rm Sgr~{}A^{\star}$. In the
paper I of this series, we showed that the eccentricity of the stars emerged
from a hypothetical disk around the SMBH can be rapidly excited by the secular
perturbation of its intermediate-mass companion (IMC), and we suggested IRS
13E as a potential candidate for the IMC. Here, we show that this process
leads to an influx of stars on parabolic orbits to the proximity of $\rm
Sgr~{}A^{\star}$ on a secular timescale of a few Myr. This timescale is much
shorter than the diffusion timescale into the lost cone through either the
classical or the resonant relaxation. Precession of the highly-eccentric
stars’ longitude of periastron, relative to that of the IMC, brings them to
its proximity within a few Myr. The IMC’s gravitational perturbation scatters
a fraction of the stars from nearly parabolic to hyperbolic orbits, with
respect to the SMBH. Their follow-up close encounters with the SMBH induce
them to escape with hyper-velocity. This scenario is a variant of the
hypothesis proposed by Hills based on the anticipated breakup of some
progenitor binary stars in the proximity of the SMBH, and its main objective
is to account for the limited lifespan of the known HVSs. We generalize our
previous numerical simulations of this process with a much wider range of
orbital configurations. We demonstrate the robustness and evaluate the
efficiency of this channel of HVS formation. From these numerical simulations,
we infer observable kinematic properties for the HVSs.

methods: numerical — stars: distances — stars: early-type — stars: kinematics
and dynamics — (stars:) planetary systems — stars: solar-type — Galaxy: center
— Galaxy: disk

††journal: ApJ††software: REBOUND (Rein & Liu, 2012; Rein et al., 2019)

##  1 Introduction

In recent years, a population of hyper-velocity stars (HVSs) has been found in
the Galactic halo (Brown et al., 2005; Brown, 2015; Zheng et al., 2014; Huang
et al., 2017). Many, if not most of them are B stars with limited main
sequence life span ($\lesssim 10^{1-2}$ Myr). They are moving relative to the
Galactic center with speeds in excess of 500 km s-1. Their inferred travel
timescale from the Galactic Center (GC) region to their present-day location
(over a distance of $\sim 10^{1-2}$ kpc) is a fraction of the main sequence
life span.

Several hypotheses have been proposed for the origin of the kinematic
properties of HVSs. The most-widely cited Hills (1988) hypothesis is based on
the pre-discovery prediction of HVSs as ejected members of tidal disrupted
close binary systems during their close encounters with super-massive black
holes (SMBH) located at $\rm Sgr~{}A^{\star}$. In order to escape from its
intense gravitational field, these binaries must acquire nearly parabolic
orbits prior to their close encounters with the SMBH (Bromley et al., 2006).
Such large eccentricities may be excited through classical or resonant
relaxation (Rauch & Tremaine, 1996) analogous to the diffusion of stars into
the loss cone in the context of tidal disruption events (Frank & Rees, 1976).
Since these diffusion processes generally require timescales considerably
longer than a few Myr, older, low-mass stars may provide a greater supply of
progenitors for the HVSs. Although there is a population of low-mass main-
sequence-star candidates with large radial velocities (Li et al., 2015; Du et
al., 2019), recent GAIA data indicate that not all but most of them are likely
to be bound to the Galaxy (Boubert et al., 2018).

Some alternative scenarios for the dynamical origin of the HVSs attribute the
scattering culprits to a hypothetical intermediate-mass companion (IMC) (Yu &
Tremaine, 2003) or some stellar-mass black-hole (O’Leary & Loeb, 2008)
companion(s) to $\rm Sgr~{}A^{\star}$. Through three-body interaction, some
HVSs, including surviving close binaries may be ejected (Lu et al., 2007).
Levin (2006) discussed the possibility of ejecting high-velocity stars from
the GC via an inspiring intermediate-mass black hole within a few dynamical
friction timescales. Another class of hypothesis accredits the HVSs’ formation
to the evolution of very close binary stars (Tauris & Takens, 1998). When the
more massive primary stars undergo an off-center explosion, their less massive
companions may be ejected from the proximity of SMBH due to the recoil motion
(Baruteau et al., 2011). Similar to Hill’s hypothesis, most of these scenarios
do not directly account for why the HVSs are mostly B stars other than the
fact that massive stars are mostly in binary and multiple systems (Chini et
al., 2012).

In an attempt to address the origin of the hyper-velocity B stars, we pay
attention to their age similarities with a population of very young massive
and Wolf-Rayet stars with estimated ages of $4-6$ Myr inside the central
$\lesssim 0.5$ pc (Ghez et al., 2003). At least half of these stars are
orbiting around the central SMBH in a relatively flat clockwise disk. These
stars are dubbed as clockwise stars (CWSs) (Löckmann & Baumgardt, 2009). The
possibility of other stars in an inclined counter-clockwise disk has also been
suggested (Paumard et al., 2006) albeit the existence of the second disk
remains controversial (Lu et al., 2006, 2009).

It has been suggested that the coeval CWSs probably formed in a common natal
disk with nearly circular orbits (Goodman, 2003; Levin & Beloborodov, 2003).
But, several of these stars are observed with large eccentricities ($\gtrsim
0.7$) (Lu et al., 2006; Paumard et al., 2006; Gillessen et al., 2009; Do et
al., 2013). Local dynamical relaxation processes cannot excite such large
eccentricities within the short time interval of a few Myr (Freitag et al.,
2006; Alexander et al., 2007; Madigan et al., 2011). Yu et al. (2007)
considered the dynamical perturbation of an inward-migrating intermediate-mass
companion (a black hole or a dark cluster) under the influence of dynamical
friction. They found that, along the perturber’s migrating path, many stars of
this stellar population could be captured into its mean motion resonance and
endure the excitation both in eccentricity and inclination. In order to
accommodate the young age of the eccentric stars in the Galactic center, this
scenario requires a short migration time scale. As a variant, Zheng et al.
(2020) suggested in Paper I of this series that the IMC’s secular perturbation
can significantly excite the eccentricity of the young stars on a time scale
comparable to their age, in the absence of any migration. Since this
interaction is long-range over many orbits, the IMC can be either a compact
stellar cluster or a point mass. We speculate that the culprit is IRS 13E
which has an estimated mass $M_{\rm IMC}\simeq 10^{4}M_{\odot}$ and is located
at a projected separation of $\sim 0.13$ pc from the $\rm Sgr~{}A^{\star}$
(Krabbe et al., 1995; Genzel et al., 2003; Maillard et al., 2004; Schödel et
al., 2005; Tsuboi et al., 2017), albeit its integrity and mass estimates
remain controversial(Wang et al., 2020; Zhu et al., 2020).

We show that the IMC’s perturbation on the stars accumulates in time near its
secular resonance, where the precession of its orbit around the SMBH matches
that of the stars induced by it. Persistent and protected angular-momentum
extraction from the stars by the IMC substantially increases their
eccentricities. During the depletion of the stars’ natal disk, the location of
secular resonance sweeps towards $\rm Sgr~{}A^{\star}$ over an extended
region, inducing high eccentricity for a large fraction of CWSs. This effect
is analogous to that induced by Jupiter on the asteroids in the solar system
(e.g., Zheng et al., 2017a, b).

The numerical results from Paper I also show that some disk stars become
unbound to the SMBH with $e_{\ast}>1$. In this paper, we carry out a
systematic study on the ejection probability of young stars formed in disks
with nearly circular orbits around $\rm Sgr~{}A^{\star}$. We consider various
contributions to the gravitational potential, including that from the super-
massive black hole (SMBH), IMC, the Galactic bulge, disk, and halo, as well as
that from a gaseous disk out of which the disk stars emerge (§2). We carry a
series of numerical simulations on the secular interaction between the IMC and
stars emerging from a hypothetical gaseous disk with circular orbits. We show
the rapid excitation of their eccentricity and the probability of close
periastron passage near the SMBH in §3. We present the results of several
series of numerical simulations for a range of initial geometry and potential
boundary condition in §4. We derive the ejection probability, the spatial
distribution, and trajectories of these HVSs. Finally in §5, we summarize our
results and discuss their implications.

##  2 Numerical models

A detailed description of the method and model parameters are given in the
paper I of this series. We briefly recapitulate various components of the
gravitational potential used in our calculation. We numerically integrate the
Newtonian equation of motion with the open-source N-body code REBOUND (Rein &
Liu, 2012) and its built-in MERCURIUS integrator (Rein et al., 2019).

###  2.1 Prescription for the gravitational potential

Table 1: Parameters for various numerical simulations

models | ${s^{\prime}_{\rm IMC}(\rm AU)}^{a}$ | $i_{\rm IMC,\star}(^{\circ})^{b}$ | $i_{\star}(^{\circ})^{c}$ | ${\rm Gas\ Disk\ Potential}^{d}$ | ${\rm Gas\ Disk\ Gap}^{e}$ | Halo Potential f
---|---|---|---|---|---|---
NGD | 0 | 0 | 0 | $-$ | $-$ | Spherical
fiducial | 0 | 0 | 0 | $+$ | $+$ | Spherical
SOFT | 5000 | 0 | 0 | $+$ | $+$ | Spherical
CCW | 0 | 180 | 0 | $+$ | $-$ | Spherical
CCW-NGD | 0 | 180 | 0 | $-$ | $-$ | Spherical
I30 | 0 | 30 | 0 | $+$ | $-$ | Spherical
I30-NGD | 0 | 30 | 0 | $-$ | $-$ | Spherical
I60 | 0 | 60 | 0 | $+$ | $-$ | Spherical
I60-NGD | 0 | 60 | 0 | $-$ | $-$ | Spherical
I60-IS60 | 0 | 60 | 60 | $+$ | $-$ | Spherical
I60-TRI | 0 | 60 | 0 | $+$ | $-$ | Triaxial
I120 | 0 | 120 | 0 | $+$ | $-$ | Spherical
I120-NGD | 0 | 120 | 0 | $-$ | $-$ | Spherical

Note. — a The softening parameter of IMC’s potential. b The initial mutual
inclination between the orbits of IMC and disk stars. c The initial
inclination of the orbit of disk stars (gas disk) relative to the Galactic
disk (reference plane).d The presence of the gas disk potential.e The presence
of a gas-free gap around the IMC’s orbit. f The form of the halo potential.

For different components of the Galactic potential, we adopt the conventional
prescription for the SMBH, IMC, Galactic bulge, disk, and halo such that

$\begin{split}\Phi_{\rm SMBH}&=-\frac{GM_{\rm SMBH}}{r}+{\rm post\ Newtonian\
terms},\\\ \Phi_{\rm IMC}&=-\frac{GM_{\rm IMC}}{\sqrt{\mid\vec{r}-\vec{r}_{\rm
IMC}\mid^{2}+s^{\prime 2}_{\rm IMC}}},\\\ \Phi_{\rm bulge}&=-\frac{GM_{\rm
bulge}}{r+a_{\rm bulge}},\\\ \Phi_{\rm disk}&=-\frac{GM_{\rm
disk}}{\sqrt{R^{2}+\left(a_{\rm disk}+\sqrt{z^{2}+d_{\rm
disk}^{2}}\right)^{2}}},\\\ \Phi_{\rm halo}&=-\frac{GM_{\rm halo}}{r}{\rm
ln}\left(1+\frac{r}{a_{\rm halo}}\right).\end{split}$ (1)

In the above expression, $r$, $R$, and $z$ are the spherical radius,
cylindrical radius and the distance above the disk; $r_{\rm IMC}$ is the
distance of the IMC from the GC. In Table 1 of Paper I, we specify $a_{\rm
bulge}(=0.6$ kpc), $a_{\rm disk}(=5$ kpc), $d_{\rm disk}(=0.3$ kpc), and
$a_{\rm halo}(=20$ kpc) are the scaling lengths for the bulge, disk, and halo
respectively; $M_{\rm SMBH}(=4\times 10^{6}M_{\odot})$ and $M_{\rm
IMC}(=10^{4}M_{\odot})$ are the mass of the SMBH and IMC, $M_{\rm
bulge}(=10^{10}M_{\odot})$, $M_{\rm disk}(4\times 10^{10}M_{\odot})$, and
$M_{\rm halo}(=10^{12}M_{\odot})$ are the mass scaling factor for the bulge,
disk, and halo, respectively (Miyamoto & Nagai, 1975; Hernquist, 1990; Navarro
et al., 1997). For the central SMBH, we include the post-Newtonian corrections
(1st and 2nd orders) near its Schwarzchild radius, same as the formulae in
Kidder (1995). As it makes negligible effects on most orbits, the $2^{\rm
nd}$-pN is only considered for those highly eccentric orbits, where the
closest distance to the central SMBH is smaller than $\sim 10^{-5}$ pc, to
avoid large energy changes during close encounters with the SMBH as discussed
in Rodriguez et al. (2018).

In one of the models (SOFT), we consider the possibility of IMC being a dense
cluster with a Plummer potential and a softening parameter $s^{\prime}_{\rm
IMC}(=5\times 10^{3}$ AU). In all models, IMC’s initial semi-major axis and
eccentricity are $a_{\rm IMC}=0.35$ pc and $e_{\rm IMC}=0.3$ respectively.

We neglect the departure from spherical symmetry for the Galactic bulge.
Although a tri-axial mass distribution in the Galactic bulge can lead to IMC’s
precession which modifies the secular evolution of single and binary stars’
orbits (Petrovich & Antonini, 2017; Bub & Petrovich, 2020), it is negligible
in comparison with the IMC’s secular perturbation on the relatively massive
main-sequence stars within the central parsec from the Galactic center during
their lifespan (i.e., over a few to tens Myr). However, a non-spherical
component in the mass distribution in the Galactic halo can significantly
modify the trajectories of HVSs on a dynamical timescale. We consider this
possibility in model I60-TRI with a triaxial potential (Yu & Madau, 2007) such
that

$\begin{split}\Phi_{\rm halo}&=-\frac{GM_{\rm halo}}{r^{\prime}}{\rm
ln}\left(1+\frac{r^{\prime}}{a_{\rm halo}}\right),\\\ {\rm with}\ \ \
r^{\prime}&=p_{\rm halo}^{1/3}q_{\rm
halo}^{1/3}\left(x^{2}+\frac{y^{2}}{p_{\rm halo}^{2}}+\frac{z^{2}}{q_{\rm
halo}^{2}}\right)^{1/2},\end{split}$ (2)

where $p_{\rm halo}$ and $q_{\rm halo}$ describe the asphericity of the halo
potential. If the $p_{\rm halo}=q_{\rm halo}=1$, the halo potential is
spherical. In I60-TRI model, we set $p_{\rm halo}=0.8$ and $q_{\rm halo}=0.7$,
with z axis normal to the Galactic disk, the same as the assumption in Yu &
Madau (2007).

###  2.2 Gaseous and stellar disks

For the gaseous disk, we adopt a prescription similar to our previous models
of protostellar disks. In the absence of a gap, the precession of the embedded
stars is driven by the gaseous disks’ central gravity

$f_{\rm emb,gas}=-4\pi G\Sigma Z_{k},$ (3)

where $Z_{k}=1.094$ is the same as in Ward (1981); Nagasawa et al. (2005);
Zheng et al. (2017a). This force is also applied to an IMC with an inclined
orbit or a counter orbiting IMC. For a co-orbiting IMC with a semi-major axis
$a_{\rm IMC}(=0.3-0.4$ pc) and eccentricity $e_{\rm IMC}(=0.2-0.4)$, we assume
it opens up a gap in the disk with a width

$a_{\rm in/out}=a_{\rm IMC}(1\mp e_{\rm IMC})\left[1\mp\left(\frac{M_{\rm
IMC}}{3M_{\rm SMBH}}\right)^{1/3}\right].$ (4)

Such that the central gravity from the gaseous disk on it becomes

$\begin{split}&f_{\rm gap,gas}=2\pi
G\Sigma\sum_{l=0}^{\infty}\left[\frac{(2l)!}{2^{2l}(l!)^{2}}\right]^{2}\left(\frac{2}{4l+1}\right)\\\
&\times\left[(2l)\left(\frac{r}{a_{\rm
out}}\right)^{2l+1/2}-(2l+1)\left(\frac{a_{\rm
in}}{r}\right)^{2l+1/2}\right].\end{split}$ (5)

For illustrative purpose, we assume a power-law surface density ($\Sigma$)
distribution constructed for the minimum mass solar nebula (Hayashi et al.,
1985) such that

$\Sigma=\Sigma_{0}\left(\frac{r}{R_{0}}\right)^{-3/2}e^{-t/\tau_{\rm
dep}}~{}\rm g/cm^{2},$ (6)

where $\tau_{\rm dep}(=3.16$ Myr) is the depletion timescale of the gas
nebula, $\Sigma_{0}(=600$ g cm-2) and $R_{0}(=10^{3}~{}\rm AU)$ are the
fiducial surface density and radii.

We assume the stars emerge from the gaseous disk and the gaseous disk
preserves its spin orientation. In one model (I60-IS60), we assume the gaseous
disk is inclined to the Galactic disk by $60^{\circ}$. In all other models,
the gaseous and Galactic disks are assumed to be in the same plane. For models
where the orbital planes of the stars, IMC, and/or the Galaxy are inclined,
the inclination of the stars relative to the gaseous disk ($i_{\ast}$) evolves
with time. The force in the ($z$) direction normal to the gaseous-disk plane
is approximated with

$f_{\rm z,gas}=-2\pi G\Sigma(R)\frac{{\rm sign}(z)(z/H)^{2}}{1+(z/H)^{2}}$ (7)

where the gaseous disk’s aspect ratio $H/R=(H_{0}/R_{0})(R/R_{0})^{1/4}$ with
$H_{0}/R_{0}=10^{-2}$.

Figure 1: The orbit evolution of a potential HVS. Color represents the time
evolution. The white cross refers to the center of our Galaxy, the grey dots
label the orbit of an IMC, with eccentricity 0.3.

###  2.3 Numerical models

In order to investigate the effects of sweeping secular resonance on the
evolution of GC stars and its parameter dependence, we test a series of
numerical models with various setups. Table 1 lists all those models,
including simulations with three different gas-disk depletion timescales. In
all cases, the stars are assumed to reside in the gaseous disk with circular
orbits initially.

We assess the influence of the sweeping secular resonance by turning off the
gravity due to the gaseous disk in a comparative model NGD. We also simulated
a SOFT model in which the IMC’s gravity is approximated by a Plummer
potential. In the fiducial, NGD, and SOFT models, we assume the $\rm
Sgr~{}A^{\star}$, companion, stars and gaseous disk are all in the same plane,
and their orbital angular momentum vectors are aligned. In the CCW model, we
consider the possibility a counter orbiting IMC as inferred from some
observations of IRS 13E (Genzel et al., 2003). In a series of 3-d cases, we
assume 30∘, 60∘, and 120∘ misalignment angle between the IMC’s and the star’s
initial orbital plane.

##  3 The secular evolution of a representative star

In Paper I, we show that IMC’s secular perturbation rapidly and robustly leads
to the eccentricity excitation of the disk stars under a wide range of
kinematic configuration. However, to the lowest order in $e_{\ast}$ and
$e_{\rm IMC}$, secular interaction leads to angular momentum exchanges with
negligible energy exchange. Under this condition, the stars’ $e_{\ast}$ is
excited without significant changes in $a_{\ast}$, and the perturbed stars
remain bound to the SMBH. However, as $e_{\ast}$ approaches to unity, higher-
order perturbation through mean-motion resonances (MMRs) intensify. MMRs
induce energy exchange between the stars and the IMC (Murray & Dermott, 2000).

It is also possible for stars’ initial $a_{\ast}=a_{0}\gtrsim(1-e_{\rm
IMC})a_{\rm IMC}/2$ to cross the IMC’s orbit as their $e_{\ast}$ approaches to
unity. The stars’ longitude of periastron ($\varpi_{\ast}$) precess mainly
under the influence of IMC’s secular perturbation, the gravity due to the
Galactic bulge and a hypothetical gaseous disk, a rate

$\frac{d\varpi_{\star}}{dt}=\frac{d\varpi_{\star}}{dt}^{(\rm
sec)}+\frac{d\varpi_{\star}}{dt}^{(\rm
bulge)}+\frac{d\varpi_{\star}}{dt}^{(\rm emb,gas)}.$ (8)

In the limit of small $e_{\ast}$, the magnitude and sign of each term on the
right hand side can be determined in terms of the model parameters (see Paper
I). At the instant of orbit crossing (when the IMC is located at $\theta_{\rm
IMC}$ and $r_{\ast}=r_{\rm IMC}$), the stars’ azimuthal position is
$\theta_{\ast}\simeq\pi\pm\varpi_{\ast}$. An close encounter between the star
and IMC would occur if $|\theta_{\ast}-\theta_{\rm IMC}|<<1$ or equivalently
$|\varpi_{\ast}-\theta_{\rm IMC}|\simeq\pi$. The strong interaction during
their close encounters (with minimum distance between the stars and the IMC,
$s_{\ast}<<a_{\ast}$) can also lead to an energy exchange between the stars
and the IMC (Lu et al., 2007).

The energy exchange between the IMC and the stars in its MMRs or endure close
encounters with it leads to modest changes in their $e_{\ast}$ and periastron
distance relative to the SMBH ($r_{\rm min}$). They become unbound from the
SMBH when their $e_{\ast}\geq 1$. We illustrate this possibility with a
representative star with an initial $a_{\ast}=a_{0}=0.2$ pc in the fiducial
model (Table 1), the orbital evolution is shown in Figure 1.

The color changes from yellow to dark in Figure 1. It delineates the time
evolution. For this representative star, its initial orbit is circular. Its
eccentricity is gradually excited by the IMC’s secular resonance, as we
suggest above. In the last few orbits shown here, its eccentricity approaches
to unity and its apoastron distance is outside the IMC’s periastron distance.
Differential precession brings the star to the IMC’s azimuthal proximity where
they undergo close encounter. The energy exchange between the star and the IMC
liberates the star from SMBH’s gravitational confinement. The star’s escape
speed is determined by its modified impact parameter relative to SMBH after
its close encounter with the IMC. Since the star ensures strong gravitation
force from the SMBH, we carry out a convergence test to ensure numerical
accuracy.

##  4 Ejection of high velocity stars

The possibility of launching hyper-velocity stars (HVSs) arises from the rapid
eccentricity excitation of disk stars within 1 pc via the SSR mechanism. After
a fraction of stars attain nearly parabolic orbits, their periastron distance
$r_{\rm min}$ is significantly reduced. Through the Hills (1988) mechanism,
any compact binary stellar systems may break up and toss out HVSs. Besides,
IMC’s nonlinear secular perturbation may also lead to the ejection of single
stars to become HVSs. Since the Hills mechanism has already been extensively
studied (Kenyon et al., 2008), we focus our investigation in this paper on the
IMC-induced ejection of single stars.

Figure 2: Left panels: The closest distance ($r_{\rm min}$) of 5000 disk stars
to the $\rm Sgr~{}A^{\star}$ as a function of their initial location
($a_{0}$). The gray dotted lines represent the tidal disruption limit for a
main-sequence star with mass $12~{}M_{\odot}$, and the gray dashed lines label
the Schwarzchild radius the SMBH. Right panels: The magnitude of the total
velocity of disk stars as a function of distance from the SMBH around 6 Myr.
Colors label the value of orbital eccentricity at the end of the simulation:
the yellow colors refer to the low-eccentricity excitation, while dark-blue
colors represent hyperbolic orbits, and green colors indicate the simulated
stars are excited to nearly parabolic orbit. Several critical mean motion
resonance locations to the IMC are labels at the left panels. The grey dashed
lines at the right panels define the criterion of an HVS, which must satisfy
$v>500~{}\rm km/s$ and $a>100$ pc.

###  4.1 The influence of sweeping secular resonance

To examine this possibility, we consider several models, which are variations
of the fiducial model (Figure 2). In all these models, we use $\Sigma_{0}=600$
g cm-2, $\tau_{\rm dep}=3.16$ Myr, $M_{\rm IMC}=10^{4}M_{\odot}$, $a_{\rm
IMC}=0.35$ pc and $e_{\rm IMC}=0.3$. For each model, we randomly place a
population of $5000$ disk stars with an initial semi-major axis $a_{0}$
randomly distributed between 0.05 pc and $a_{\rm in}(\sim 0.22$ pc) (Equation
4). We first analyze the results of the NGD model which neglect the effect of
the sweeping secular resonance. In this case, large eccentricity excitation
occurs only for stars with orbits which are initially around 2:1 and 5:2 mean
motion resonances with that of the IMC (at $\sim 0.22$ and $\sim 0.19$ pc
respectively). Very few stars are ejected into hyperbolic orbits with
asymptotic speed $v>500$ km s-1, after their $r_{\rm min}$ is reduced to $\sim
1$ AU (top panel in Figure 2).

Figure 3: Close encounters immediately prior to ejection. Along the x axis,
HVSs’ closest distance from the IMC $r_{\rm IMC,\ast}=|{\bf r}_{\rm IMC}-{\bf
r}|$ is measured in units of pc. The y axis is their
$\Delta\theta_{\rm\ast,IMC}=\theta_{\rm IMC}-\omega_{\ast}$ immediately after
the instant of their closest encounters (for stars on parabolic orbits, $\pi$
corresponds to colocation with the IMC). The ejected HVSs’ distances from $\rm
Sgr~{}A^{\star}$ at the end of the calculation are indicated by color.

In the fiducial model (middle panel in Figure 2), the evolving gaseous disk
introduces the effect of sweeping secular resonance. After 6 Myr ($\sim
2~{}\tau_{\rm dep}$), a population of disk stars in the region $0.15-0.22$ pc
from the SMBH attains hyperbolic orbits ($e_{\star}>1$) and reaches distance
$\gtrsim 100$ pc with asymptotic speed with a magnitude $v>500$ km s-1 mostly
in the direction away from the SMBH. We refer these stars as the candidate
HVSs. We have not included, in our models, compact binary systems which are
the progenitors of HVSs under the Hills hypothesis. This ejection channel for
HVSs through sweeping secular resonance occurs rapidly during the epoch of
disk depletion. It can account for HVSs’ young ages and high masses. Moreover,
a minute fraction of stars attain $r_{\rm min}$ comparable to the tidal
disruption radius (Frank & Rees, 1976), $r_{\rm tide}\simeq(M_{\rm
SMBH}/m_{\star})^{1/3}s_{\ast}(\approx 3.5\times 10^{-6}$ pc). For all those
very close approaching events, we re-check them with higher accuracy. These
rare close encounters raise the possibility of infrequent tidal disruption
events.

The HVSs’ asymptotic speed in the fiducial model is $\sim 1\%$ of the escape
speed for hypothetical stars on circular orbits at $r_{\rm min}$. However, it
is slightly larger than IMC’s orbital speed at IMC’s periastron. This
“slingshot” is accomplished in three steps (see discussions in §3 and
illustrative Figure 1). 1) Most of the “heavy lifting” is accomplished through
IMC’s secular perturbation which induces some stars to attain nearly parabolic
orbits (with very little orbital energy). 2) These marginally bound stars
receive some impulsive “gravitational assist” from the IMC during their close
encounters. Although the fractional change on the stellar orbits ($\sim(M_{\rm
IMC}a_{\rm IMC}/M_{\rm SMBH}r_{\rm IMC,\ast})^{1/2}$ where $r_{\rm
IMC,\ast}=|{\bf r}_{\rm IMC}-{\bf r}|$ is the stars’ distance from the IMC) is
modest, the IMC and SMBH’s gravity together may significantly modify some
stars’ energy and impact parameter relative to the SMBH. 3) Some of the
scattered stars undergo an additional close encounter with the SMBH and escape
from its gravitational confinement with $e_{\star}>>1$.

We select a group of HVSs and plot their distance $r_{\rm IMC,\ast}$ and
angular separation $\Delta\theta_{\rm\ast,IMC}=|\theta_{\rm
IMC}-\omega_{\ast}|$ from the IMC during or immediately after the instant of
their closest encounters (Figure 3). Dispersion of
$\Delta\theta_{\rm\ast,IMC}$ from $\pi$ is due to the scattering during the
closest encounters and it can attain non-negligible values for stars with
orbits marginally overlapping that of the IMC. The ejected HVSs’ distances
from $\rm Sgr~{}A^{\star}$ at the end of the calculation are indicated by
color. The lack of correlation between the minimum star-IMC distance and HVSs’
distances from $\rm Sgr~{}A^{\star}$ is an indication that the IMC is not the
primary culprit of hyper-velocity ejection. In contrast to the correlation
between $r_{\rm min}$ and the HVSs’ final distances at the end of the
simulation (Figures 2 and 4) provides supporting evidences to suggest that
their scattering by the SMBH after their impulsive close encounters with the
IMC was the event that led to their escape with hyper velocities.

Figure 4: Same as Figure 2, with a comparison between 5000 disk stars in the
CCW and the CCW-NGD models.

In model SOFT, we adopt a Plummer model for the companion’s gravitational
potential. The softening parameter $s_{\rm IMC}^{\prime}\sim 0.07~{}a_{\rm
IMC}$ is marginally smaller than its Hills radius $r_{\rm H}=(M_{\rm
IMC}/3M_{\rm SMBH})^{1/3}a_{\rm IMC}\sim 0.09~{}a_{\rm IMC}$. This modified
potential weakens the impact of star-IMC close encounters without sacrificing
the IMC’s distant secular perturbation on the stars. In addition, we turn off
the merging possibility between the IMC and disk stars only in this SOFT
model. The bottom panel of Figure (2) indicates that $e_{\ast}$ for a large
fraction of stars is excited to $\sim 1$. With their nearly zero total orbital
energy, close encounters with IMC’s relatively shallow potential are adequate
to perturb these stars from parabolic to hyperbolic orbits. The result of the
SOFT model suggests that the IMC needs not to be an IMBH, and the distributed
potential of a compact stellar cluster is a sufficient agitator for the HVSs.

###  4.2 Counter orbiting IMC

We also simulate a CCW model with the identical gas and stellar disks as the
fiducial model. The IMC’s $M_{\rm IMC}$, $a_{\rm IMC}$, and $e_{\rm IMC}$ are
also identical except its orbit is in the counter-clockwise direction (i.e.,
its orbital angular momentum is anti-parallel to that of the stars and gas
disk). A wealth of disk stars in more extended regions are significantly
excited (with $e_{\star}>>1)$ by this counter-orbiting IMC, albeit in the 2-D
limit. However, those disk stars initially distributed within $\sim 0.13$ pc
are less perturbed compared to the fiducial model and SOFT model in Figure 2).
Moreover, the mean motion resonances are less pronounced in this model, as
shown in Figure 4. The CCW model can produce many more ejections and HVSs than
the fiducial/SOFT models. The frequency of HVSs is close to $1\%$ among all
disk stars (shown in Figure 7). This result is consistent with the analytic
calculation in the paper I of this series, where we showed that a counter-
orbiting IMC could indeed induce secular resonances over a wide region.
Although the presence of gas disk introduces additional precession which
initially reduces the effectiveness of IMC’s secular perturbation, the
influence of a counter orbiting IMC is reinstalled after the disk depletion.
In another series of CCW-NGD model, a counter orbiting IMC is imposed without
the gas disk potential (i.e., with $\Sigma_{0}=0$). The results in Figure 4
show that it is just as capable to excite $e_{\star}$ and eject HVS’s.

Figure 5: Same as Figure 2 with a comparison between the I30, I30-NGD, I60,
and I60-NGD models.

###  4.3 IMC on an inclined orbit

Figure 6: The $(e,i)$ distribution of simulated stars around 6 Myr as a
function of distance to the GC in the I30, I30-NGD, I60, I60-NGD, I120, and
I120-NGD models. Colors label the eccentricity of disk stars, and $i$ refers
to the inclination of the orbits of disk stars and IMC relative to the
Galactic disk plane (reference plane), range from $0^{\circ}$ to
$180^{\circ}$, where $90^{\circ}$ means that the orbits of disk stars are
perpendicular to the Galactic disk. The black dots represent the IMC’s final
semi-major axis and inclination (relative to the reference plane).

In most of the above models, we assume the simulated disk stars formed in the
mid-plane of a gas disk, and they are coplanar to a co-orbiting or a counter
orbiting IMC. Due to the uncertainties in the observed orbit of IRS 13E, a
more general analysis is warranted. Although we can generalize the analytic
calculation in the paper I, many degrees of freedom expand its complexity.
Instead, we adopt additional illustrative series of numerical simulations. For
models with evolving gas disks, we use the same value of $\Sigma$ distribution
(but without a gap), $\Sigma_{0}$, and $\tau_{\rm dep}$ as the fiducial model.

In the I30, the orbits of all the disk stars are initially set to be in the
same plane with $i_{\ast}=0^{\circ}$. We also adopt a mutual initial
inclination of $i_{\rm IMC}-i_{\star}=30^{\circ}$ for IMC’s orbit . The IMC’s
mass, semi-major axes, and eccentricity ($M_{\rm IMC}$, $a_{\rm IMC}$, and
$e_{\rm IMC}$) in this and other models with inclined stellar orbits are
identical to those in the fiducial model. Similarly, the model I60 and I120
are set with the initial mutual inclination between disk stars and the IMC, to
be $60^{\circ}$ and $120^{\circ}$, respectively. For all these inclined
models, disk stars are initially embedding in the gaseous disk and coplanar to
the Galactic disk (refers to as reference plane), but the IMC is in an
inclined orbit. Accordingly, the related NGD model, I30-NGD, I60-NGD, and
I120-NGD models, describe the same misalignment between the stars and IMC but
without considering the existence of gas disk.

Compare to the 2-d cases, the sweeping secular resonance of an inclined IMC
can still significantly excite disk stars albeit less effective, especially
for the highly inclined orbit, e.g., I60 model. Even though there is no gas
disk under the coplanar circumstance, the IMC’s secular perturbation and mean
motion resonances can stimulate some hyper-velocity escapees in the NGD model.
However, in the slightly inclined case without gas disk, model I30-NGD, only a
few disk stars can be incredibly excited with a periapse distance smaller than
$10^{-4}$ pc. As a result, no HVS was ejected. The situation is improved in
the I60-NGD model. In the absence of a gaseous disk, the sweeping secular
resonance effect does not contribute to any HVS population in this case. The
mutual inclination between the IMC and the disk stars is larger than the
threshold value of Lidov-Kozai secular resonance around a point mass
potential, which is generally set as $39^{\circ}$. Even though the Galactic
bulge potential introduces additional precession, more disk stars are excited
and ejected in the I60-NGD model than the I30-NGD model.

In Figure 6, we focus on the inclination distribution of disk stars under each
inclined model after $\sim 6$ Myr evolution. For a clearer contrast, the left
panels of Figure 6 show the results of inclined cases with the gas disk, while
the right panels are the corresponding models without the gas disk. As
expected, the presence of a gas disk not only introduces sweeping secular
resonance which leads to efficient $e_{\ast}$ excitation and ejection of HVSs
but also damps the stars’ excited inclination relative to the disk and
preserve their initial $i_{\ast}$, especially for those residuals at
$r_{\ast}<0.1$ pc (see also in the paper I).

Figure 7: The HVSs frequency in each model. The blue dots label the models
without gas disk, while the yellow dots represent the models with gas disk.
The x-axis, $i_{\rm IMC,\star}$, refers to the initial mutual inclination
between the orbits of IMC and disk stars.

In the paper I, we show that the characteristic timescales for secular
evolution for the Galactic center stars depend weakly on $a_{\ast}$. Since
they are comparable to the simulated duration ($\sim 6$ Myr) for all the
models, typical stars, including those with $a_{\ast}\lesssim 0.1$ pc, went
through approximately one to a few secular cycles in both $e_{\ast}$ and
$(i_{\ast}-i_{\rm IMC})$. This pattern is evident in the relatively low
$e_{\ast}$ and the gradient of $(i_{\ast}-i_{\rm IMC})$ after $\sim 6$ Myr.
The range of $e_{\ast}$ modulation is small for $a_{\ast}<0.1$ pc. However,
the eccentricity of some more distant stars (with $a_{\ast}>0.1$ pc) are
significantly excited. These stars also endure large changes in $i_{\ast}$
from their initial values. Due to their contribution to the potential, changes
in $i_{\ast}$ in models with a gas disk are generally smaller than those
without gas disk (i.e., the NGD models).

Figure 8: The top panels are the same as Figure 6, but with a comparison
between the I60-IS60 and the I60-TRI models. $i_{,\rm gas}$ refers to the
inclination of the orbits of disk stars (and IMC) relative to the plane of the
gas disk, range from $0^{\circ}$ to $180^{\circ}$. The bottom panels keep
tracing these escapees’ final trajectories in the I60-IS60 and I60-TRI models
around 60 Myr.  Figure 9: Left panels: It is the instantaneous trajectories of
ejected stars with distance $>1$ kpc in the I60-TRI model around 60 Myr. The
z-axis is the azimuth angle, $\phi$, which is the angle between ${\bf z}$ and
${\bf v}$. The y-axis is the angle $\theta$ shows between the direction of
${\bf r}$ and ${\bf v}$. Right panels: The $\theta({\bf r,v})$ of these
ejected stars as a function of time (from $\sim 6$ Myr to $\sim 60$ Myr). The
color indicates the distance from the Galactic center.  Figure 10: Spatial
distributions of 10000 simulated stars in the I60, I60-IS60, and I60-TRI
models at 6 Myr (left panels: a, c, and e) and 60 Myr (right panels: b, d, and
f, only stars with the velocity $>500$ km/s are plotted). All disk stars are
projected in the standard Galactic coordinates. The initial orbital plane of
disk stars is presented with thick grey line. Their distances to the GC are
labeled with colors. For better comparison, those faraway, ejected HVSs are
marked with larger sizes, proportional to the logarithm of their distances to
the GC.

The potential for the gas disk modifies the frequency of both longitudinal and
nodal precessions. The disk potential diminishes on a timescale $\tau_{\rm
dep}$ which is also comparable to the duration of the simulation and the
eccentricity excitation timescale in models I30, I60, and I120. Direct
comparisons between these models with their diskless counter parts (models
I30-NGD, I60-NGD, and I120-NGD) indicate that changes in both $e_{\ast}$ and
$i_{\ast}$ are reduced for the close-in stars (with $a_{\ast}<0.1$ pc) by the
disk potential (Figures 5 and 6). Outside $\sim 0.1$ pc, changes in the disk
potential during its depletion significantly enlarges the radial extent of
IMC’s secular resonance (cf. equivalent models without the evolving gas
disks). Significant eccentricity excitation (leading to $r_{\rm min}<10^{-3}$
pc) and inclination changes occur for stars with $a_{0}\sim 0.15-0.22$ pc,
similar to the fiducial and SOFT models. This region is populated by stars
with a wide range of inclination and relatively high eccentricity. This result
may be linked to the observed inclination diversity for the young stars around
$\rm Sgr~{}A^{\star}$ (Paumard et al., 2006; Lu et al., 2006, 2009). While the
eccentricity of many stars is excited to values comparable to that observed
among the CCWs, few stars acquire $r_{\rm min}$ sufficiently small $(\lesssim
10^{-4}~{}\rm pc)$ to undergo close encounter with the SMBH and then be
ejected as HVSs. Although the frequency of HVSs is increased by the presence
of the disk, it is also reduced by IMC’s inclined orbit (Figure 7).

###  4.4 HVSs’ trajectory in the Galactic halo

At a distance greater than $\sim 1$ kpc, the HVSs’ orbital trajectories are
under the influence of the Galactic disk and halo. In all the above inclined
models, we assume the initial orbits of all the stars (including those
eventually ejected as HVS’s) are in the same plane as the Galactic disk. We
carried out another model in which we consider the possibility that the
stellar disk is inclined $60^{\circ}$ to the Galactic disk (Paumard et al.,
2006) in Figures 8 and 10. It shows that HVSs’ asymptotic trajectories are
indeed randomly distributed rather than cluster around the Galactic disk’s
pole.

In all but one of the simulations, we adopt a spherically symmetric potential
for the Galactic halo. Tri-axial mass distribution in the Galactic halo (Wegg
et al., 2019) can lead to the randomization of HVSs’ orbits relative to
Galactic plane and the reorientation of their trajectories from radially
outward direction relative to the Galactic center (Yu & Madau, 2007; Brown,
2015; Du et al., 2019). In Figures 8, we plot the stars’ orbital inclination
relative to their natal disk and their Galactico-centric distance. Their
position in the Galactic longitude and latitude, viewed from the Galactic
Center is show in Figure 10, in a similar way as Lu et al. (2010).

In model I60-TRI, we use Equation 2 to approximate a tri-axial mass
distribution and potential. We plot the distribution of the angles

$\theta({\bf r,v})={\rm cos}^{-1}\left({{\bf r\cdot v}\over|r||v|}\right)\ \ \
\ \ \ {\rm and}\ \ \ \ \ \ \phi({\bf z,v})={\rm cos}^{-1}\left({{\bf z\cdot
v}\over|z||v|}\right)$ (9)

distribution for the HVSs after 60 Myr (left panel of Figure 9). This sample
of 14 HVSs is chosen with $r>10^{2}$ pc and $v>500$ km s-1 as in Figure 2. Six
HVSs which are escaping the Galactic halo reach $30-60$ kpc with relatively
small $\theta({\bf r,v})$, i.e. their trajectories have not deviated
significantly from the radial direction. Seven HVSs’s distance remain within
$\sim 10$ kpc. Their $\theta({\bf r,v})$ are widely distributed as shown by Yu
& Madau (2007). Although the tri-axial potential has randomized their
trajectory, there is no obvious indication of clustering along the minor axis
of the halo potential.

We also plot the evolution of $\theta({\bf r,v})$ (right panel of Figure 9)
for the same HVS sample. The most distant stars are also the fastest moving
stars. Their trajectories have not deviated significantly from radial motion
away from the Galactic center. However, the orientation of the closer HVS
appears to evolve under the influence of non-spherical halo potential. Thus,
we do not expect the observed proper motion of all the HVSs to be directed
away from the Galactic center.

We also find that in the CCW model (not plotted here), most of the escaping
stars retain their initial retrograde orbits (with $i_{\star}=180^{\circ}$)
though the magnitude of their $e_{\star}$ is extensively excited. However,
some escaping disk stars flip their orbital angular momentum vector to follow
the orbital direction of the IMC. Since both the disk stars and the IMC are in
the same orbital plane in the CCW model, no star attains an inclined orbit.
Nevertheless, misalignment between the orbit of the disk stars with respect to
the Galactic plane and a tri-axial mass distribution in the halo can also
change the spatial distribution of HVS’s in this model.

##  5 Summary & Discussion

The main motivation of this investigation is to explore the origin of a
population of hyper-velocity B stars (HVSs). These stars appear to be escaping
from the Galaxy with velocity $>500$ km/s. Their special distribution in the
Galaxy is widely spread out. The similarities of their spectral class and life
expectancy ($\sim 10^{7}$ yr) with those of the S stars and disk stars near
$\rm Sgr~{}A^{\star}$ naturally raise the possibility that they may have been
ejected from the proximity of the SMBH at the center of the Milky Way.

The central conundrum is the mechanisms that enabled the HVSs to gain such
enormous kinetic energy and to venture out of the deep gravitational potential
imposed by the SMBH. A widely accepted scenario proposed by Hills (1988)
suggests that an infusion of close binary stars to the proximity of an SMBH
can lead to their tidal disruption and the ejection of one component with
hyper-velocity. The main drawback of this scenario in the present context of
the Galactic HVSs is the long timescales required for the diffusion of binary
stars through stellar relaxation processes (Lightman & Shapiro, 1977; Rauch &
Tremaine, 1996; Merritt, 2013).

In the paper I of this series, we show a promising mechanism for rapid
eccentricity excitation near $\rm Sgr~{}A^{\star}$. We attribute the main
cause to be secular perturbation from an intermediate-mass companion to the
super-massive black hole. We speculate the IMC culprit to be IRS 13E. We show
the robustness of this effect for a wide range of IMC’s orbital configuration
especially if we consider the evolution of the gravitational potential during
the depletion of their natal disk. We provide analytic calculation and
numerical simulation to show that recently formed (several Myr ago) stars can
indeed attain the large observed eccentricities. We showed that the IMC’s
secular perturbation is over a considerable distance and extensive time such
that it can be either a compact stellar cluster or an IMBH with a mass $\sim
10^{4}M_{\odot}$.

In this paper, we extend the IMC scenario for the HVSs’ past dynamical
evolution. With numerical simulations, we show that the IMC’s secular
perturbation is effectively exciting the eccentricity of some young stars to
near unity. Through this action at a distance, the orbital energy of these
nearly parabolic stars essentially reduces the minute fraction of their
gravitational potential energy at their periastron. Some stars’ periastron
approaches to (but do not cross) their tidal disruption radius while their
apoastron crosses the IMC’s orbit.

Due to the combined perturbation by the IMC and the Galactic bulge, the stars’
eccentric orbits precess. Some stars undergo close encounters with the IMC and
endure impulsive 3-body interaction which modifies their trajectory,
especially the impact parameter, $1-e_{\ast}$, and energy of their orbital
energy relative to the SMBH. In the next close encounter with the SMBH, some
stars are slingshot to escape IMC’s gravitational confinement as HVSs.

We examine the robustness of this process for a range of IMC’s orbital
configurations. The ejection probability ranges from a few to a fraction of a
percent depending on whether the orbits of the IMC and the disk stars are in
the same plane (in either a counter or co-orbiting configuration) or inclined
planes. Although the presence and depletion of a natal disk promote the
ejection of the HVSs, they are not essential. The asymptotic trajectories of
the HVSs are widely distributed rather than concentrated in a plane over the
sky. If the halo potential is tri-axial, the instantaneous HVSs’ direction of
motion need not be pointed away from the $\rm Sgr~{}A^{\star}$, albeit the
fastest moving HVSs may follow such pathways.

This IMC-secular-perturbation mechanism differs from the binary-break-up
mechanism suggested by Hills since all the stars in our models are single.
Nevertheless, three body interaction (between the star, IMC, and SMBH) plays a
pivotal role in the liberation of the HVSs (Yu & Tremaine, 2003). In this
work, we include an additional ingredient of IMC’s secular perturbation and
sweeping secular resonance as the main mechanism for rapidly supplying
recently formed stars to the neighborhood of the SMBH. The main advantage of
this scenario for the origin of HVSs is that it can naturally account for the
youthfulness of the HVSs.

Since the periastron distance of the escaping stars are generally well outside
the tidal disruption radius (Figures 2 and 4), some young close binary stars
may survive their close encounter with the SMBH as they are ejected to become
HVSs. The detection of any such binary HVSs can provide a discriminative test
since they are excluded by the Hill’s mechanism (Lu et al., 2007).

Finally, the results presented here are robust for a wide range of IMC’s
kinematic properties. Similar companions may be common in other galaxies,
especially if galaxies evolved through mergers (a basic concept of
$\Lambda$CDM scenario) and central massive black holes in galaxies are
ubiquitous. These massive black holes undergo orbital decay toward the center
of their host galaxies under the action of dynamical friction and gas drag.
During their approach toward coalescence, similar dynamical effects between
migrating massive black holes may excite the eccentricity of their surrounding
stars as well as stars along their paths. The subsequent rapid injection of
stars to the proximity of the massive black holes can significantly modify the
rate of tidal disruption events, the formation and coalescence of a binary
star and stellar-mass black holes (as sources of gravitational wave events),
and the ejection of hyper-velocity stars and stellar-mass black holes. We end
with speculation that such processes may even lead to the possibility of
freely floating intermediate-mass black holes outside any host galaxies. The
follow-up investigation of these related issues will be presented elsewhere.

Simulations in this paper made use of the REBOUND code which is publicly
available at REBOUND website. We thank Tuan Do, Stefan Gillessen, Mark Morris,
and Hangci Du for useful conversation. This work is partly supported by the
National Key Research and Development Program of China (No. 2018YFA0404501 to
SM), by the National Science Foundation of China (Grant No. 11821303,
11761131004 and 11761141012 to SM) and by the China Postdoctoral Science
Foundation (Grant No. 2017M610865 to XCZ).

## References

  * Alexander et al. (2007) Alexander, R. D., Begelman, M. C., & Armitage, P. J. 2007, ApJ, 654, 907, doi: 10.1086/509709
  * Baruteau et al. (2011) Baruteau, C., Cuadra, J., & Lin, D. N. C. 2011, ApJ, 726, 28, doi: 10.1088/0004-637X/726/1/28
  * Boubert et al. (2018) Boubert, D., Guillochon, J., Hawkins, K., et al. 2018, MNRAS, 479, 2789, doi: 10.1093/mnras/sty1601
  * Bromley et al. (2006) Bromley, B. C., Kenyon, S. J., Geller, M. J., et al. 2006, ApJ, 653, 1194, doi: 10.1086/508419
  * Brown (2015) Brown, W. R. 2015, ARA&A, 53, 15, doi: 10.1146/annurev-astro-082214-122230
  * Brown et al. (2005) Brown, W. R., Geller, M. J., Kenyon, S. J., & Kurtz, M. J. 2005, ApJ, 622, L33, doi: 10.1086/429378
  * Bub & Petrovich (2020) Bub, M. W., & Petrovich, C. 2020, ApJ, 894, 15, doi: 10.3847/1538-4357/ab8461
  * Chini et al. (2012) Chini, R., Hoffmeister, V. H., Nasseri, A., Stahl, O., & Zinnecker, H. 2012, MNRAS, 424, 1925, doi: 10.1111/j.1365-2966.2012.21317.x
  * Do et al. (2013) Do, T., Martinez, G. D., Yelda, S., et al. 2013, ApJ, 779, L6, doi: 10.1088/2041-8205/779/1/L6
  * Du et al. (2019) Du, C., Li, H., Yan, Y., et al. 2019, ApJS, 244, 4, doi: 10.3847/1538-4365/ab328c
  * Frank & Rees (1976) Frank, J., & Rees, M. J. 1976, MNRAS, 176, 633, doi: 10.1093/mnras/176.3.633
  * Freitag et al. (2006) Freitag, M., Amaro-Seoane, P., & Kalogera, V. 2006, ApJ, 649, 91, doi: 10.1086/506193
  * Genzel et al. (2003) Genzel, R., Schödel, R., Ott, T., et al. 2003, ApJ, 594, 812, doi: 10.1086/377127
  * Ghez et al. (2003) Ghez, A. M., Duchêne, G., Matthews, K., et al. 2003, ApJ, 586, L127, doi: 10.1086/374804
  * Gillessen et al. (2009) Gillessen, S., Eisenhauer, F., Fritz, T. K., et al. 2009, ApJ, 707, L114, doi: 10.1088/0004-637X/707/2/L114
  * Goodman (2003) Goodman, J. 2003, MNRAS, 339, 937, doi: 10.1046/j.1365-8711.2003.06241.x
  * Hayashi et al. (1985) Hayashi, C., Nakazawa, K., & Nakagawa, Y. 1985, in Protostars and Planets II, ed. D. C. Black & M. S. Matthews, 1100–1153
  * Hernquist (1990) Hernquist, L. 1990, ApJ, 356, 359, doi: 10.1086/168845
  * Hills (1988) Hills, J. G. 1988, Nature, 331, 687, doi: 10.1038/331687a0
  * Huang et al. (2017) Huang, Y., Liu, X. W., Zhang, H. W., et al. 2017, ApJ, 847, L9, doi: 10.3847/2041-8213/aa894b
  * Kenyon et al. (2008) Kenyon, S. J., Bromley, B. C., Geller, M. J., & Brown, W. R. 2008, ApJ, 680, 312, doi: 10.1086/587738
  * Kidder (1995) Kidder, L. E. 1995, Phys. Rev. D, 52, 821, doi: 10.1103/PhysRevD.52.821
  * Krabbe et al. (1995) Krabbe, A., Genzel, R., Eckart, A., et al. 1995, ApJ, 447, L95, doi: 10.1086/309579
  * Levin (2006) Levin, Y. 2006, ApJ, 653, 1203, doi: 10.1086/507830
  * Levin & Beloborodov (2003) Levin, Y., & Beloborodov, A. M. 2003, ApJ, 590, L33, doi: 10.1086/376675
  * Li et al. (2015) Li, Y.-B., Luo, A. L., Zhao, G., et al. 2015, Research in Astronomy and Astrophysics, 15, 1364, doi: 10.1088/1674-4527/15/8/018
  * Lightman & Shapiro (1977) Lightman, A. P., & Shapiro, S. L. 1977, ApJ, 211, 244, doi: 10.1086/154925
  * Löckmann & Baumgardt (2009) Löckmann, U., & Baumgardt, H. 2009, MNRAS, 394, 1841, doi: 10.1111/j.1365-2966.2009.14466.x
  * Lu et al. (2006) Lu, J. R., Ghez, A. M., Hornstein, S. D., et al. 2006, in Journal of Physics Conference Series, Vol. 54, Journal of Physics Conference Series, 279–287, doi: 10.1088/1742-6596/54/1/044
  * Lu et al. (2009) Lu, J. R., Ghez, A. M., Hornstein, S. D., et al. 2009, ApJ, 690, 1463, doi: 10.1088/0004-637X/690/2/1463
  * Lu et al. (2007) Lu, Y., Yu, Q., & Lin, D. N. C. 2007, ApJ, 666, L89, doi: 10.1086/521708
  * Lu et al. (2010) Lu, Y., Zhang, F., & Yu, Q. 2010, ApJ, 709, 1356, doi: 10.1088/0004-637X/709/2/1356
  * Madigan et al. (2011) Madigan, A.-M., Hopman, C., & Levin, Y. 2011, ApJ, 738, 99, doi: 10.1088/0004-637X/738/1/99
  * Maillard et al. (2004) Maillard, J. P., Paumard, T., Stolovy, S. R., & Rigaut, F. 2004, A&A, 423, 155, doi: 10.1051/0004-6361:20034147
  * Merritt (2013) Merritt, D. 2013, Dynamics and Evolution of Galactic Nuclei
  * Miyamoto & Nagai (1975) Miyamoto, M., & Nagai, R. 1975, PASJ, 27, 533
  * Murray & Dermott (2000) Murray, C. D., & Dermott, S. F. 2000, Solar System Dynamics (Cambridge University Press), doi: 10.1017/CBO9781139174817
  * Nagasawa et al. (2005) Nagasawa, M., Lin, D. N. C., & Thommes, E. 2005, ApJ, 635, 578, doi: 10.1086/497386
  * Navarro et al. (1997) Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490, 493, doi: 10.1086/304888
  * O’Leary & Loeb (2008) O’Leary, R. M., & Loeb, A. 2008, MNRAS, 383, 86, doi: 10.1111/j.1365-2966.2007.12531.x
  * Paumard et al. (2006) Paumard, T., Genzel, R., Martins, F., et al. 2006, ApJ, 643, 1011, doi: 10.1086/503273
  * Petrovich & Antonini (2017) Petrovich, C., & Antonini, F. 2017, ApJ, 846, 146, doi: 10.3847/1538-4357/aa8628
  * Rauch & Tremaine (1996) Rauch, K. P., & Tremaine, S. 1996, New A, 1, 149, doi: 10.1016/S1384-1076(96)00012-7
  * Rein & Liu (2012) Rein, H., & Liu, S. F. 2012, A&A, 537, A128, doi: 10.1051/0004-6361/201118085
  * Rein et al. (2019) Rein, H., Tamayo, D., & Brown, G. 2019, MNRAS, 489, 4632, doi: 10.1093/mnras/stz2503
  * Rodriguez et al. (2018) Rodriguez, C. L., Amaro-Seoane, P., Chatterjee, S., et al. 2018, Phys. Rev. D, 98, 123005, doi: 10.1103/PhysRevD.98.123005
  * Schödel et al. (2005) Schödel, R., Eckart, A., Iserlohe, C., Genzel, R., & Ott, T. 2005, ApJ, 625, L111, doi: 10.1086/431307
  * Tauris & Takens (1998) Tauris, T. M., & Takens, R. J. 1998, A&A, 330, 1047
  * Tsuboi et al. (2017) Tsuboi, M., Kitamura, Y., Tsutsumi, T., et al. 2017, ApJ, 850, L5, doi: 10.3847/2041-8213/aa97d3
  * Wang et al. (2020) Wang, Q. D., Li, J., Russell, C. M. P., & Cuadra, J. 2020, MNRAS, 492, 2481, doi: 10.1093/mnras/stz3624
  * Ward (1981) Ward, W. R. 1981, Icarus, 47, 234, doi: 10.1016/0019-1035(81)90169-X
  * Wegg et al. (2019) Wegg, C., Gerhard, O., & Bieth, M. 2019, MNRAS, 485, 3296, doi: 10.1093/mnras/stz572
  * Yu et al. (2007) Yu, Q., Lu, Y., & Lin, D. N. C. 2007, ApJ, 666, 919, doi: 10.1086/520622
  * Yu & Madau (2007) Yu, Q., & Madau, P. 2007, MNRAS, 379, 1293, doi: 10.1111/j.1365-2966.2007.12034.x
  * Yu & Tremaine (2003) Yu, Q., & Tremaine, S. 2003, ApJ, 599, 1129, doi: 10.1086/379546
  * Zheng et al. (2017a) Zheng, X., Lin, D. N. C., & Kouwenhoven, M. B. N. 2017a, ApJ, 836, 207, doi: 10.3847/1538-4357/836/2/207
  * Zheng et al. (2017b) Zheng, X., Lin, D. N. C., Kouwenhoven, M. B. N., Mao, S., & Zhang, X. 2017b, ApJ, 849, 98, doi: 10.3847/1538-4357/aa8ef3
  * Zheng et al. (2020) Zheng, X., Lin, D. N. C., & Mao, S. 2020, The Astrophysical Journal, 905, 169, doi: 10.3847/1538-4357/abc8e5
  * Zheng et al. (2014) Zheng, Z., Carlin, J. L., Beers, T. C., et al. 2014, ApJ, 785, L23, doi: 10.1088/2041-8205/785/2/L23
  * Zhu et al. (2020) Zhu, Z., Li, Z., Ciurlo, A., et al. 2020, ApJ, 897, 135, doi: 10.3847/1538-4357/ab980d