The Astrophysical Journal, 794:106 (23pp), 2014 October 20 doi:10.1088/0004-637X/794/2/106
C
2014. The American Astronomical Society. All rights reserved. Printed in the U.S.A.
ON THE DISTRIBUTION OF STELLAR REMNANTS AROUND MASSIVE BLACK HOLES: SLOW MASS
SEGREGATION, STAR CLUSTER INSPIRALS, AND CORRELATED ORBITS
Fabio Antonini
Canadian Institute for Theoretical Astrophysics, University of Toronto, 60 St. George Street, Toronto, Ontario M5S 3H8, Canada; antonini@cita.utoronto.ca
Received 2014 February 19; accepted 2014 August 13; published 2014 September 29
ABSTRACT
We use N-body simulations as well as analytical techniques to study the long-term dynamical evolution of stellar
black holes (BHs) at the Galactic center (GC) and to put constraints on their number and mass distribution. Starting
from models that have not yet achieved a state of collisional equilibrium, we find that timescales associated with
cusp regrowth can be longer than the Hubble time. Our results cast doubts on standard models that postulate high
densities of BHs near the GC and motivate studies that start from initial conditions that correspond to well-defined
physical models. For t he first time, we consider the distribution of BHs in a dissipationless model for the formation
of the Milky Way nuclear cluster (NC), in which massive stellar clusters merge to form a compact nucleus. We
simulate the consecutive merger of ∼10 clusters containing an inner dense sub-cluster of BHs. After the formed NC
is evolved for ∼5 Gyr, the BHs do form a steep central cusp, while the stellar distribution maintains properties that
resemble those of the GC NC. Finally, we investigate the effect of BH perturbations on the motion of the GC S-stars
as a means of constraining the number of the perturbers. We find that reproducing the quasi-thermal character of
the S-star orbital eccentricities requires 1000 BHs within 0.1 pc of Sgr A*. A dissipationless formation scenario
for the GC NC is consistent with this lower limit and therefore could reconcile the need for high central densities
of BHs (to explain the S-stars orbits) with the “missing-cusp” problem of the GC giant star population.
Key words: galaxies: nuclei – Galaxy: center – Galaxy: formation – stars: black holes –
stars: kinematics and dynamics
Online-only material: color figures
1. INTRODUCTION
Over the whole Hubble sequence, massive nuclear clus-
ters (NCs) are observed at the center of many galaxies. The
frequency of nucleation among galaxies less luminous than
∼10
10.5
L
is close to 90% as determined by ACS HST ob-
servations of galaxies in the Virgo and Fornax galaxy clusters
(Carollo et al. 1998;B
¨
oker et al. 2002;C
ˆ
ot
´
eetal.2006; Turner
et al. 2012). The study of NCs is of great interest for our un-
derstanding of galaxy formation and evolution as indicated by
the fact that a number of fairly tight correlations are observed
between their masses and global properties of their host galax-
ies such as velocity dispersion and bulge mass (Ferrarese et al.
2006; Wehner & Harris 2006; Graham & Spitler 2009; Scott &
Graham 2013; Leigh et al. 2012). Intriguingly, similar scaling
relations are obeyed by massive black holes (MBHs), which
are predominantly found in massive galaxies that, however,
show little evidence of nucleation (e.g., Graham & Spitler 2009;
Neumayer & Walcher 2012). The existence of such correlations
might indicate a direct link among large galactic spatial scales
and the much smaller scale of the nuclear environment and sug-
gests that NCs contain information about the processes that have
shaped the central regions of their host galaxies.
How NC formation takes place at the center of galaxies is
still largely debated (e.g., Hartmann et al. 2011; Gnedin et al.
2014; Carlberg & Hartwick 2014; Mastrobuono-Battisti et al.
2014). Relatively recent work has shown that “dissipationless”
models can reproduce without obvious difficulties the observed
properties (Turner et al. 2012) and scaling relations (Antonini
2013) of NCs. In these models an NC forms through the
inspiral of massive stellar clusters into the center because of
dynamical friction where they merge to form a compact nucleus
(e.g., Tremaine et al. 1975; Capuzzo-Dolcetta & Miocchi 2008;
Capuzzo-Dolcetta 1993). Alternatively, NCs could have formed
locally as a result of radial gas inflow into the galactic center
accompanied by efficient dissipative processes (Schinnerer
et al. 2008; Milosavljevi
´
c 2004). Naturally, dissipative and
dissipationless processes are not exclusive and both could be
important for the formation and evolution of NCs (Hartmann
et al. 2011; Antonini et al. 2012; De Lorenzi et al. 2013).
The Milky Way NC, being only 8 kpc away, is currently
the only NC that can be resolved in individual stars and
for which a kinematical structure and density profile can be
reliably determined (Genzel et al. 2010). This offers the unique
possibility to resolve the stellar population in order to study the
composition and dynamics close to an MBH and put constraints
on different NC formation scenarios. The Milky Way NC has
an estimated mass of ∼10
7
M
(Launhardt et al. 2002; Sch
¨
odel
et al. 2009), and it hosts a massive black hole of ∼4 × 10
6
M
(Genzel et al. 2003; Ghez et al. 2008; Gillessen 2009), whose
gravitational potential dominates over the stellar cusp potential
out to a radius of roughly 3 pc—the MBH radius of influence. A
handful of other galaxies are also known to contain both an NC
and an MBH, which typically have comparable masses (Seth
et al. 2008). Population synthesis models suggest that roughly
80% of the stellar mass i n the inner parsec of the Milky Way
is in (>5 Gyr) old stars (Pfuhl et al. 2011) although the light is
dominated by the young stars. This appears to also typically be
the case in most NCs observed in external galaxies (Rossa et al.
2006).
Over the last decade, observations of the Galactic NC have
led to a number of puzzling discoveries. These discoveries
include the presence of a young population of stars (the
S stars) near Sgr A* in an environment extremely hostile to star
formation (paradox of youth; Morris 1993; Sch
¨
odel et al. 2002)
and a significant paucity of red giant stars in the inner half of a
1
The Astrophysical Journal, 794:106 (23pp), 2014 October 20 Antonini
parsec (conundrum of old age; Merritt 2010). Number counts of
the giant stars at the Galactic center (GC) show that their visible
distribution is in fact quite inconsistent with the distribution
of stars expected for a dynamically relaxed population near a
dominating Keplerian potential (Buchholz et al. 2009;Doetal.
2009; Bartko et al. 2010): instead of a steeply rising Bahcall
&Wolf(1976) cusp, there is a ∼0.5 pc core. The lack of
a Bahcall–Wolf cusp in the giant distribution casts doubts on
dynamically relaxed, quasi-steady-state models of the GC that
postulate a high central density of stars and stellar black holes
(BHs). In these models the central distribution of stars and BHs
is determined by just a handful of parameters: the MBH mass,
the total density outside the relaxed region, and the slope of the
initial mass function (IMF; Merritt 2013). Given the unrelaxed
form of the density profile of stars, making predictions about
the distribution of the stellar remnants becomes a much more
challenging, time-dependent problem susceptible to the initial
conditions and to the (yet largely unconstrained) formation
process of the NC (Antonini & Merritt 2012).
Understanding the distribution of the “stellar remnants” in
systems similar to the Milky Way’s NC is crucial in many
respects. Examples include randomization of the S-star orbits
via gravitational encounters (Perets et al. 2009), warping of the
young stellar disk (Kocsis & Tremaine 2011), and formation
of X-ray binaries (Muno et al. 2005). Stellar nuclei similar to
that of the Milky Way are also the location of astrophysical
processes that are potential gravitational wave (GW) sources
both for ground- and space-based laser interferometers. These
include the merger of compact object binaries near MBHs
(Antonini & Perets 2012) and the capture of BHs by MBHs,
called “extreme mass ratio inspirals” (EMRIs; Amaro-Seoane
et al. 2012). The efficiency of these dynamical processes and
rate estimates for GW sources are very sensitive to the number
of BHs near the center. Therefore, a fundamental question is
whether, when given a prediction for the initial distribution of
stars and BHs, the system is old enough that the heavy remnants
had time to relax and segregate to the center of the Galaxy.
Motivated by the above arguments, we consider the long-term
evolution of BH populations at the center of galaxies, starting
from different assumptions regarding their initial distribution.
Since the stellar BHs at the GC are not directly detected, time-
dependent numerical calculations, like the ones presented below,
are crucial for understanding and making predictions about the
distribution of stellar remnants at the center of galaxies.
In Section 2 we explore the evolution of models in which stars
and BHs follow initially the same spatial distribution, which
is far from being in collisional equilibrium. Contrary to some
previous claims (Preto & Amaro-Seoane 2010), we find that in
these models the time to regrow a cusp in both the BH and the
star distribution is longer than the age of the Galaxy. For realistic
number fractions of BHs, our simulations demonstrate that over
the age of the Galaxy, the presence of a heavy component has
little effect on the evolution of the stellar component.
In Sections 3, 4, and 5, we discuss the evolution of BHs
in a globular cluster merger model for NCs. We present the
results of direct N-body simulations of the merger of globular
clusters containing two mass populations: stars and BHs. These
systems were in an initial state of mass segregation with the BH
population concentrated toward the cluster core. Each cluster
was placed on a circular orbit with a galactocentric radius
of 20 pc in a N-body system containing a central MBH. We
find that the inspiral of massive globular clusters in the center
of the Galaxy constitutes an efficient source term of BHs in
these regions. After about 10 inspiral events, the BHs are highly
segregated to the center. After a small fraction of the nucleus
relaxation time (as defined by the main stellar population), the
BHs attain a nearly steady-state distribution; at the same time
the stellar density profile exhibits a ∼0.2 pc core, similar to
the size of the core in the distribution of s tars at the GC. Our
results indicate that standard models, which assume the same
initial phase space distribution for BHs and stars, can lead to
misleading results regarding the current dynamical state of the
Galactic center.
We discuss the implications of our results in Section 6.In
particular, we show that in order to reproduce the quasi-thermal
form of the observed eccentricity distribution of the S-star orbits,
about 1000 BHs should be present inside ∼0.1 pc of Sgr A*.
This number appears to be consistent with the number of BHs
expected in a model in which the Milky Way NC formed through
the orbital decay and merger of about 10 massive clusters.
Our main results are summarized in Section 7.
2. SLOW MASS SEGREGATION AT
THE GALACTIC CENTER
In this section we study the long-term dynamical evolution
of multi-mass models for the Milky Way NC. The primary goal
of this study is to understand the evolution of the distribution of
stars and BHs over a time of order the central relaxation time
of the nucleus, starting from initial conditions that are far from
being in collisional equilibrium.
2.1. Evolution toward the Steady State
We consider four mass groups representing main-sequence
stars (MSs), white dwarfs (WDs), neutron stars (NSs), and BHs.
After the quasi steady state is attained, the stars are expected to
follow a central r
−3/2
cusp, while the heavier particles will have
a steeper r
−2
density profile (e.g., Alexander 2005). We assume
that all species have the same phase space distribution initially
as would be expected for a violently relaxed system. This is
the assumption that was made in most previous papers (e.g.,
Freitag et al. 2006; Hopman & Alexander 2006; Merritt 2010).
We specify the mass ratio, m
wd
/m
= 0.6, m
ns
/m
= 1.4,
and m
bh
/m
= 10, between the mass group particles and
respective number fractions, f
wd
= N
wd
/N
, f
ns
= N
ns
/N
,
and f
bh
= N
bh
/N
.
Number counts of the old stellar population at the GC are
consistent with a density profile of stars that is flat or slowly
rising toward the MBH inside its sphere of influence and within
a radius of roughly ∼0.5 pc (Buchholz et al. 2009;Doetal.2009;
Bartko et al. 2010). Outside this radius the density falls off as
r
−2
. Merritt (2010) showed that a core of size ∼0.5 pc is a natural
consequence of two-body relaxation acting over 10 Gyr, starting
from a core of radius ∼1 pc. It is therefore of interest to study the
evolution of the BH distribution for a time of order the age of the
Galaxy and starting from a density distribution with a parsec-
scale core. We adopt the truncated broken-power-law model:
ρ(r) = ρ
0
r
r
0
−γ
i
1+
r
r
0
α
(γ
i
−γ
e
)/α
ζ (r/r
cut
), (1)
where ζ (x) = (2/sech(x) + cosh(x)), α is a parameter that de-
fines the transition strength between inner and outer power
laws, r
0
is the scale radius, and r
cut
is the truncation radius
of the model. The values adopted for these parameters were
2
The Astrophysical Journal, 794:106 (23pp), 2014 October 20 Antonini
r
0
= 1.5 pc, α = 4, γ
e
= 1.8, and r
cut
= 6 pc. We included
a central MBH of mass M
•
= 4 × 10
6
M
and generated the
models N-body representations via numerically calculated dis-
tribution functions. The central slope was set to γ
i
= 0.6, the
smallest density slope index consistent with an isotropic distri-
bution for the adopted density model and potential.
The normalizing factor ρ
0
was chosen in such a way that the
corresponding density profile reproduces the coreless density
model
ρ(r) = 1.5 × 10
5
r
1pc
−1.8
M
pc
−3
(2)
outside the core. This choice of normalizing constant gives a
mass density at 1 pc similar to what is inferred from observations
(e.g., Oh et al. 2009) and gives a total mass in stars within this
radius of ∼1.6 × 10
6
M
. The fact that our models are directly
scalable to the observed stellar density distribution of stars at
the GC is important if we want to draw conclusions about the
current dynamical state of stars and BHs at the GC. We note,
for example, that the merger models of Gualandris & Merritt
(2012) had core radii that were substantially larger than the
MBH influence radius. As also noted by these authors, this
simple fact precluded a unique scaling of their models to the
Milky Way—at least in the Galaxy’s current state in which the
stellar core size (∼0.5 pc) is much smaller than the Sgr A*
influence radius (∼3 pc).
We run three simulations with N = 132 k particles. These
simulations differ with each other by the adopted number
fractions of the four mass groups: (1) f
wd
= f
ns
= f
bh
= 0; (2)
f
wd
= 10
−1
, f
ns
= 10
−2
, f
bh
= 10
−3
; and (3) f
wd
= 2 × 10
−1
,
f
ns
= 2 × 10
−2
, f
bh
= 5 × 10
−3
. The latter two set of values
correspond roughly to the number fractions expected from a
standard and from a top-heavy IMF, respectively. A fraction of
f
bh
= 10
−3
is what is expected for a standard (Kroupa-like)
IMF, and it is the value typically adopted in previous studies
(e.g., Hopman & Alexander 2005, 2006). Although a larger
fraction of stellar remnants might be possible, for instance,
if the Galactic center always obeyed a top-heavy initial mass
function, the observationally constrained mass-to-light ratio of
the inner parsec limits the BH fraction to only a few percent and
is more consistent with a ratio and a total mass of BHs predicted
by a standard IMF (L
¨
ockmann et al. 2010). We evolved these
systems for a time equal to the relaxation time, T
r
infl
, computed
at the sphere of influence of the MBH. The relaxation time was
evaluated using the expression (Spitzer 1987):
T
r
=
0.34σ (r)
3
G
2
mln Λρ(r)
, (3)
where ρ is the total local mass density, and m is the average
particle mass. For the Coulomb logarithm we used ln Λ =
ln(r
infl
σ
2
/2Gm
) ≈ 10, with σ the 1d velocity dispersion
outside r
infl
= GM
•
/σ
2
.
To scale the N-body time length t o the Milky Way, we consider
that the relaxation time at the influence radius of Sgr A*,
r
infl
≈ 3 pc, is T
r
infl
≈ 25 Gyr, assuming a stellar mass of 1 M
(Merritt 2010; Antonini & Merritt 2012). Thus, when scaling to
the GC, a time of 0.4 T
r
infl
corresponds to roughly 10 Gyr.
We evolved the initial conditions with the direct N-body
integrator φGRAPEch (Harfst et al. 2008). The code im-
plements a fourth-order Hermite integrator with a predictor-
corrector scheme and hierarchical time-stepping. The code
combines hardware-accelerated computation of pairwise inter-
particle forces (using the Sapporo library, which emulates the
GRAPE interface utilizing GPU boards; Gaburov et al. 2009)
with a high-accuracy chain regularization algorithm to f ollow
the dynamical interactions of field particles with the central
MBH particle. The chain radius was set to 10
−2
pc, and we used
a softening = 10
−6
pc. The relative error in total energy was
typically ∼10
−4
for the accuracy parameter η = 0.01.
Figure 1 shows the evolution of the N-body models over
one relaxation time. The heavy particles segregate to the center
because of dynamical friction. After the central mass density of
BHs becomes comparable to the density in the other species,
the evolution of the BH population starts being dominated by
BH–BH self-interactions; at the same time the lighter species
evolve in response to dynamical heating from the BHs, which
causes the local stellar densities to decrease and Lagrangian
radii t o expand. As shown below, the same heating rapidly
converts the initial density profile into a steeply rising density
cusp with slope, γ ≡−d log ρ/d log r ≈ 3/2. The inclusion of
a BH population has therefore two effects on the main-sequence
population: it lowers the stellar densities and at the same time
it accelerates the evolution of the density of stars toward the
γ = 3/2 steady-state form.
The lower panels of Figure 1 display the density profile of
stars and BHs over 10 Gyr of evolution. These plots show that,
starting with a fraction of BHs that corresponds to a standard
IMF, (1) after ∼10 Gyr the density of BHs can remain well
below the density of stars at all radii, and (2) even after 10 Gyr
of evolution, the density distribution of stars looks very dif-
ferent from what expected for a dynamically relaxed popula-
tion around an MBH. These findings are in agreement with the
Fokker–Plank simulations of Merritt (2010) but in contrast with
more recent claims that mass segregation can rebuild a stel-
lar cusp in a relatively small fraction of the Hubble time (e.g.,
Preto & Amaro-Seoane 2010, and the Introduction of Chen &
Amaro-Seoane 2014). Figure 2 displays the evolution of the ra-
dial profile of the density profile slope. Comparing the evolution
observed in models with and without BHs, we see that a cusp in
the main-sequence population develops earlier in models with
BHs. Figure 2 shows that for f
bh
= 10
−3
and f
bh
= 5 × 10
−3
,
a stellar cusp only develops after ∼0.6 T
r
infl
and ∼0.4 T
r
infl
,re-
spectively. Therefore, over the timescales (10 Gyr) and radii
(r 0.01 pc) of relevance, the inclusion of a BH population has
little or even no influence on the evolution of the lighter popula-
tions. This latter point is more clearly demonstrated in Figure 3,
which directly compares the Lagrangian radii evolution of our
f
bh
= 0 model with models with BHs. The stellar populations
evolve similarly in these models independently of f
bh
until ap-
proximately 0.6 and 0.4 × T
r
infl
for f
bh
= 10
−3
and 5 × 10
−3
,
respectively. After this time, heating of the lighter species by the
heavy particles starts becoming important, causing the density of
the former to decrease and deviate from the evolution observed
in the single-mass component model. However, the transition
to this phase clearly occurs after the models have been already
evolved for a time comparable to (for f
bh
= 5 ×10
−3
) or longer
than (for f
bh
= 10
−3
) the age of the Galactic NC.
1
Given the results of the simulations presented in this section,
we can schematically divide mass segregation in two phases: in
phase 1 the density of BHs is smaller than the density of stars and
the models evolve mainly because of scattering off the stars—the
BHs inspiral to the center because of dynamical friction and the
stellar distribution relaxes as in a single-mass component model.
1
The mean stellar age in the Galactic NC is estimated to be ∼5 Gyr (Figer
et al. 2004).
3
The Astrophysical Journal, 794:106 (23pp), 2014 October 20 Antonini
Figure 1. Top panels show the Lagrangian radii for the four stellar species during the N-body simulations with BH number fractions: f
bh
= 10
−3
(left panels) and
f
bh
= 5 × 10
−3
(right panels). Top tick marks give times after scaling to the Milky Way; we adopted a relaxation time at the Sgr A* influence radius of 25 Gyr (e.g.,
Merritt 2010). Bottom panels show the density profile of stars and BHs at t = (0, 2.5, 5, 7.5, 10) Gyr; the central density increases with time. Clearly, even after a
time of the order of 10 Gyr, the distribution of stars and BHs in our models can be very different from the relaxed multi-mass models that are often used to describe
the center of galaxies. The vertical line marks the MBH influence radius.
(A color version of this figure is available in the online journal.)
Figure 2. Evolution of the density slope, γ ≡−d log ρ/d log r, of the main-sequence density profile in multi-mass N-body models (center and right panels) compared
with a model with only one mass component (left panel). The continued curves show profiles at (0.2, 0.4, 0.6, 0.8, 1) × T
r
infl
; increasing line width corresponds to
increasing time. The dashed curve corresponds to the initial model. Adding a BH component accelerates the growth of a density cusp in the stellar component.
However, the time to regrow a cusp in these models is always longer than 0.2 T
r
infl
, i.e., 5 Gyr when scaled to the Milky Way, a time longer than the mean stellar age
of the Galactic NC.
4
The Astrophysical Journal, 794:106 (23pp), 2014 October 20 Antonini
Figure 3. Evolution of the f
bh
= 0 model Lagrangian radii and density profile compared with the evolution of models with f
bh
= 10
−3
(left panel) and f
bh
= 5 ×10
−3
(right panel). In the bottom panels the density profile is plotted at t = (0, 5, 10, 15, 20, 25) Gyr of evolution. The blue dashed curves give the evolution of the stellar
distribution Lagrangian radii in the models with BHs, and the bottom panels show the respective stellar density profiles at t = (0, 5, 10) Gyr. Over a time of order
10 Gyr, the evolution of the density profile of stars in the f
bh
= 10
−3
model is not much affected by the presence of the BHs. In the model with f
bh
= 5 × 10
−3
the
evolution toward the steady state is faster, and after 10 Gyr the stars have formed a cusp. Vertical lines give the MBH influence radius.
(A color version of this figure is available in the online journal.)
In phase 2, when the density of BHs becomes comparable and
larger than the density of stars, BHs and stars evolve because
of scattering off the BHs, which causes the models to rapidly
evolve toward the steady state.
Perhaps the most interesting aspect of our simulations is the
long timescale required by the BH population to segregate to
the center through dynamical friction (phase 1 above) and reach
a (nearly) steady-state distribution—a time comparable to the
relaxation time as defined by the dominant stellar population.
In what follows we show that these predictions agree well with
the evolution expected on the basis of theoretical arguments.
2.2. Analytical Estimates
In order to understand the evolution of the distribution of
BHs observed in the N-body simulations, we evolved the
population of massive remnants using an analytical estimate
for the dynamical friction coefficient. The stellar background
was represented as an analytic potential that was also allowed to
evolve in time according to the evolution observed in the stellar
distribution during the N-body simulations.
We began by generating random samples of positions and ve-
locities from the isotropic distribution function corresponding to
the density model of Equation (1). The orbital equations of mo-
tion were then integrated forward in time in the evolving smooth
stellar potential and included a term that describes the orbital
energy dissipation due to dynamical friction. The dynamical
friction acceleration was computed using the expression
a
fr
≈−4πG
2
m
bh
ρ(r)
v
v
3
ln Λ
v
0
dv
4πf (v
)v
2
+
∞
v
dv
4πf (v
)v
2
ln
v
+ v
v
− v
− 2
v
v
, (4)
where v is the velocity of the inspiraling BH, and f (v
)theve-
locity distribution of field stars. The second term in parenthesis
5
