ON THE DISTRIBUTION OF STELLAR REMNANTS AROUND MASSIVE BLACK HOLES: SLOW MASS SEGREGATION, STAR CLUSTER INSPIRALS, AND CORRELATED ORBITS

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⁻² 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; Do et al. 2009; Bartko et al. 2010). Outside this radius the density falls off as r⁻². 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:

$$\begin{equation} \rho (r)=\rho _0 \left(\frac{r}{r_0} \right)^{-\gamma _i} \left[ 1+\left(\frac{r}{r_0} \right)^{\alpha }\right]^{(\gamma _i-\gamma _{e})/{\alpha }}\zeta (r/r_{\rm cut}), \end{equation} \tag{ 1 }$$

where ζ(x) = (2/sech(x) + cosh(x)), α is a parameter that defines the transition strength between inner and outer power laws, r₀ is the scale radius, and r_cut is the truncation radius of the model. The values adopted for these parameters were r₀ = 1.5 pc, α = 4, γ_e = 1.8, and r_cut = 6 pc. We included a central MBH of mass M_• = 4 × 10⁶ M_☉ and generated the models N-body representations via numerically calculated distribution functions. The central slope was set to γ_i = 0.6, the smallest density slope index consistent with an isotropic distribution for the adopted density model and potential.

The normalizing factor ρ₀ was chosen in such a way that the corresponding density profile reproduces the coreless density model

$$\begin{equation} \rho (r)=1.5\times 10^5 \left(\frac{r}{1\ {\rm pc}} \right)^{-1.8} \, M_{\odot } \ {\rm pc}^{-3} \end{equation} \tag{ 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⁶ 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⁻¹, f_ns = 10⁻², f_bh = 10⁻³; and (3) f_wd = 2 × 10⁻¹, f_ns = 2 × 10⁻², f_bh = 5 × 10⁻³. 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⁻³ 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öckmann et al. 2010). We evolved these systems for a time equal to the relaxation time, $T_{r_{\rm infl}}$, computed at the sphere of influence of the MBH. The relaxation time was evaluated using the expression (Spitzer 1987):

$$\begin{equation} T_{r}=\frac{0.34\sigma (r)^3}{G^2 \langle m\rangle {\rm ln}\, \Lambda \rho (r)}, \end{equation} \tag{ 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σ²/2Gm_⋆) ≈ 10, with σ the 1d velocity dispersion outside r_infl = GM_•/σ².

To scale the N-body time length to the Milky Way, we consider that the relaxation time at the influence radius of Sgr A*, r_infl ≈ 3 pc, is $T_{r_{\rm infl}}\approx 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_{\rm infl}}$ corresponds to roughly 10 Gyr.

We evolved the initial conditions with the direct N-body integrator ϕGRAPEch (Harfst et al. 2008). The code implements a fourth-order Hermite integrator with a predictor-corrector scheme and hierarchical time-stepping. The code combines hardware-accelerated computation of pairwise interparticle 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 follow the dynamical interactions of field particles with the central MBH particle. The chain radius was set to 10⁻² pc, and we used a softening = 10⁻⁶ pc. The relative error in total energy was typically ∼10⁻⁴ 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 to expand. As shown below, the same heating rapidly converts the initial density profile into a steeply rising density cusp with slope, γ ≡ −dlog ρ/dlog 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.

Zoom In Zoom Out Reset image size

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 different from what expected for a dynamically relaxed population 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 stellar 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 radial 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⁻³ and f_bh = 5 × 10⁻³, a stellar cusp only develops after ${\sim }0.6\,T_{r_{\rm infl}}$ and ${\sim }0.4\,T_{r_{\rm infl}}$, respectively. 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 populations. 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 approximately 0.6 and $0.4 \times T_{r_{\rm infl}}$ for f_bh = 10⁻³ and 5 × 10⁻³, 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⁻³) or longer than (for f_bh = 10⁻³) the age of the Galactic NC.¹

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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. 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.

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 velocities from the isotropic distribution function corresponding to the density model of Equation (1). The orbital equations of motion 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

$$\begin{eqnarray} \boldsymbol{a_{\rm fr} } &\approx & -4\pi G^2 m_{\rm bh} \rho (r) \frac{\boldsymbol{v}}{v^3} \left ( {\rm ln}\, \Lambda \int _{0}^v dv_\star 4\pi f(v_\star) v_\star ^2\right. \nonumber\\ &&\left. +\int ^{{\infty }}_v dv_\star 4\pi f(v_\star) v_\star ^2 \left[ {\rm ln} \left(\frac{v_\star +v}{v_\star -v} \right)-2\frac{v}{v_\star } \right] \right), \end{eqnarray} \tag{ 4 }$$

where v is the velocity of the inspiraling BH, and f(v_⋆) the velocity distribution of field stars. The second term in parenthesis of Equation (4) represents the frictional force due to stars moving faster than the test mass (Chandrasekhar 1943). Such "non-dominant" terms are neglected in the standard Fokker–Plank treatment in which the dynamical friction coefficient is obtained by integrating only over field stars with velocity smaller than that of the test particle (e.g., Hopman & Alexander 2006; Alexander & Hopman 2009; Merritt 2010). Antonini & Merritt (2012) showed that this approximation breaks down in a shallow density profile of stars around an MBH where such terms can become dominant, as there are a few or even no particles moving more slowly than the local circular velocity.

The N-body integrations show that the stellar distribution changes with time in a quasi-self-similar way—the stellar density profile break radius shrinks progressively with time while the outer profile slope is maintained roughly unchanged. In order to account for such evolution, we computed at each time the best-fitting density model of Equation (1) to the density profile of stars in the N-body system at that time. We used this density model to compute gravitational potential, distribution function, and corresponding dynamical friction coefficient. This procedure allowed us to include the evolution of the stellar background when evolving the BH population. Our integrations are unique in the sense that they are the first ones to simultaneously include (1) a correct estimate of the dynamical friction coefficient, which takes into account the contribution of stars moving faster than the inspiraling BH, and (2) a realistic treatment of the evolution of the stellar background under the influence of gravitational encounters. However, since our analysis does not take into account BH–BH self-interactions, our integrations are only valid until the density of BHs remains well below the density of stars. In this respect, our approach is limited to the early evolution of the system, when the BHs only represent a negligible perturbation on the evolution of the light component.

The upper panel of Figure 4 shows the density profile of BHs obtained from the semi-analytical modeling described above and compares it with the results from the direct N-body simulation with the BH fraction f_bh = 10⁻³. The central density of BHs increases with time at a rate that is comparable in the two models. The plot shows that the spatial distribution of stellar-mass BHs near the GC might not have reached a steady-state form—at least if their initial distribution was similar to what used in our models. In fact, even after a time of ≈8 Gyr, the central density of BHs is still substantially lower than the density of stars. As our analysis demonstrates, the persistence of such low densities of BHs is a direct consequence of the long timescale of inspiral in a density core near an MBH. This latter point is illustrated in the lower panel of Figure 4, which shows the trajectory of a 10 M_☉ BH at the GC. The rate of orbital decay slows down as the BH reaches ≈r₀/2 because of the lack of low-velocity stars in the core. After the BH reaches this radius, dynamical friction becomes very inefficient and the decay of the BH orbit proceeds at a rate that is comparable to the rate at which the core radius in the stellar distribution shrinks because of gravitational encounters—a time of order the relaxation time of the nucleus.

Zoom In Zoom Out Reset image size

In this section we used direct N-body integrations as well as analytic models to describe the evolution of multi-mass models of the Milky Way NC characterized by an initial parsec-scale core in the density distribution. Calculations similar to those described here were recently performed by Preto & Amaro-Seoane (2010) and Gualandris & Merritt (2012).

Preto & Amaro-Seoane (2010) studied the evolution of models with two mass species: stars and BHs. These authors concluded that mass segregation of the heavy component speeds up cusp growth in the lighter component by factors up to 10 in comparison with the single-mass case. This conclusion is somewhat in agreement with the results of our simulations, which also show that a stellar cusp, extending out to roughly r_infl, regrows faster in models with BHs (e.g., Figure 2). However, for realistic numbers of BHs we find that the timescale of cusp regrowth is only a factor of two shorter than in the single-mass component models. Preto & Amaro-Seoane (2010) argued that the timescales associated with cusp regrowth are clearly shorter than the Hubble time for nuclei similar to that of the Milky Way—even though the relaxation time, as estimated for a single-mass stellar distribution, exceeds one Hubble time. On the basis of our study, we conclude instead that over one Hubble time, and if a standard IMF is adopted, adding a heavy component has relatively little effect on the evolution of the main-sequence component (e.g., Figure 3). Even for a top-heavy IMF, which results in initial larger densities of BHs, the time for cusp regrowth is longer than the mean stellar age in the Galactic center (∼5 Gyr). The reason for this is that because of the inefficient dynamical friction force in a density core around an MBH, the central density of BHs remains well below the density of stars for a time of order the relaxation time of the nucleus. The time required to regrow a cusp in the stellar distribution appears to be longer than the Hubble time for galaxies similar to the Milky Way.

Gualandris & Merritt (2012) simulated the merger between galaxies with MBHs containing four mass groups, representative of old stellar populations. They followed the evolution of the merger products for about three relaxation times and found that the density cores formed during the galaxy mergers persisted and that the distribution of the stellar-mass black holes evolved "against an essentially fixed stellar background." Gualandris & Merritt (2012) also integrated the exact same Fokker–Planck models as in Preto & Amaro-Seoane (2010) and argued that the accelerated cusp growth described by these latter authors is seen to be present only at small radii, r ≲ 0.05r_infl. At radii larger than these, adding the BHs has the effect of lowering the density of the stellar component at all times. Gualandris & Merritt (2012) argued that Preto & Amaro-Seoane (2010) were misled by looking at the very-small-radius regime in their Fokker–Plank solutions, where the cusp in the main-sequence component stands out. Our study shows that the BH population has indeed two effects on the main-sequence population: it lowers the "mean" density of stars (the point stressed in Gualandris & Merritt 2012), and it accelerates the redistribution of the stars in the phase space toward the γ = 3/2 steady state (as found in Preto & Amaro-Seoane 2010). So in a sense, the BHs both "create" and "destroy" a cusp: although the presence of a BH population can significantly accelerate the timescale of cusp regrowth in the stellar distribution, the scattering off the heavy (BH) component causes the density of stars to decrease at radii larger than ∼0.05r_infl.

Sufficiently massive and compact clusters can decay toward their parent galaxy central region in a time much shorter than the Hubble time. An approximation of the time for clusters (within the half mass radius of the stellar bulge) to spiral to the center is given by (Antonini et al. 2012):

$$\begin{equation} \Delta t_{h} \approx 3 \times 10^{10} \ {\rm yr} \tilde{R}_{{\rm eff}}^{1.75} \tilde{m}^{-1} , \end{equation} \tag{ 5 }$$

where $\tilde{R}_{{\rm eff}}$ is the galactic effective radius in kiloparsecs, and $\tilde{m}$ is the globular cluster mass in units of 10⁶ M_☉. Within Δt_h the forming nucleus has a luminosity comparable to that of the surviving clusters. Equation (5) predicts that a significant fraction of the globular cluster population in faint and intermediate luminosity stellar spheroids would have spiraled to the center by now, while in giant ellipticals, because of their larger characteristic radii, the time required to grow an NC might be longer than 10¹⁰ yr. We stress that the inspiral time obtained using Equation (5) gives only a crude approximation (likely an overestimate) of the real dynamical friction timescale. Nevertheless, it is reasonable to draw the conclusion that the observed lack of nuclei in galaxies more massive than about 10¹⁰ M_☉ could be due to the longer infall times in these galaxies, due to their larger values of R_eff. Figure 5 presents a test of this idea. This figure gives effective radii versus masses, M_gal, for galaxies belonging to the Virgo cluster that either have (filled black circles) or do not have (star symbols) a central NC. Dashed curves give the value of R_eff obtained by setting Δt_h = 10¹⁰ yr in Equation (5) and adopting various masses of the sinking object. The figure shows that only in galaxies with M_gal ≲ 10¹⁰ M_☉, massive globular clusters would have enough time to spiral into the center, merge, and form a compact nucleus. The observed absence of compact nuclei in giant ellipticals could be therefore interpreted as a consequence of the long dynamical friction timescale of globular clusters in these galaxies.

Zoom In Zoom Out Reset image size

We add that the density profile of stars in giant ellipticals is often observed to be flat or slowly rising inside the influence radius of the MBH. As shown in Section 2 this implies a very long dynamical friction timescale inside the MBH influence radius due to the absence of stars moving more slowly than the local circular velocity (Antonini & Merritt 2012). Massive clusters orbiting within the core of a giant elliptical galaxy do not reach the center even after 10¹⁰ yr. In addition, because of the strong tidal field produced by the MBH, globular clusters can only transport little mass to the very central region of the galaxy. Both these effects, i.e., long inspiral times and little mass transported to the center, have been suggested to suppress the formation of NCs in bright galaxies, in agreement with observations (Antonini 2013).

A useful reference time for our study is the relaxation time computed at the radius containing half of the mass of the system, R. Setting ln Λ = 12, and ignoring the possible presence of an MBH, the half mass relaxation time is

$$\begin{equation} T_h= 2.1 \times 10^5\frac{ [ r_h(\rm pc)]^{3/2} N^{1/2}}{(m/M_{\odot })^{1/2}}\frac{10}{\ln (0.4 N)}\ {\rm yr}, \end{equation} \tag{ 6 }$$

where N is the total number of stars.

In the absence of an MBH, collisional relaxation leads to mass segregation and core collapse. In a preexisting NC, the presence of an MBH inhibits core collapse, causing instead the formation of a Bahcall–Wolf cusp, n ∼ r^−7/4, on the two-body relaxation timescale (Preto et al. 2004; Merritt 2009). Nuclear clusters belonging to the Virgo galaxy cluster have a half-mass relaxation time that scales with the total absolute magnitude of the host galaxy, M_B, as (Merritt 2009):

$$\begin{equation} {\rm log}\left(T_h/{\rm yr}\right)=9.38\hbox{--}0.43(M_B+16). \end{equation} \tag{ 7 }$$

Galaxies with luminosities less than ∼4 × 10⁸ L_☉ have NCs with relaxation times that fall below 10 Gyr. These galaxies have NCs with masses ≲  10⁷ M_☉ and half mass radii r_h ≲ 10 pc. These limiting values appear close to those characterizing the Milky Way NC, suggesting that only spheroids fainter than the Milky Way have collisionally relaxed nuclei. Relaxation times for nuclei with masses ≳  10⁷ M_☉ are therefore too long for assuming that they have reached a collisionally relaxed state, but they are still short enough that gravitational encounters would substantially affect their structure over the Hubble time. This is in agreement with the results of Section 2 and also appears to be consistent with absence of a Bahcall–Wolf cusp in the distribution of stars at the GC.

Globular clusters with N ∼ 10^{6 − 7} have relaxation times T_h ∼ 10^{9 − 10} yr. Most Galactic globular clusters are therefore relaxed systems. The timescale required for the BH population to segregate to the cluster center and form a subcluster dynamically decoupled from the host stellar cluster is approximately the core-collapse time for the initial BH cluster (e.g., Banerjee et al. 2010):

$$\begin{equation} T_{ms}\approx \frac{m}{m_{bh}}T_{cc}, \end{equation} \tag{ 8 }$$

where m_bh is the mass of a stellar black hole, and T_cc is the core-collapse time of the host stellar cluster, which is about T_cc ≈ 15 × T_h for a Plummer model. After ∼T_ms the central density of BHs becomes large enough that BH–BH binary formation takes place through three or four body interactions (Heggie & Hut 2003). The formed BH binaries then "harden" through repeated superelastic encounters that lead to the ejection of BHs from the cluster core until eventually only a few BHs are left (e.g., Banerjee et al. 2010).

In galaxies similar to the Milky Way, stellar clusters with masses ≳  10⁶ M_☉ and starting from a galactocentric radius of 1 kpc have orbital decay times due to dynamical friction less than ≲ 3 Gyr (Equation (5)). The cluster's dynamical friction time is therefore typically long compared with the timescale over which the BHs would segregate to the center of the cluster.

It is possible that the BH population will evaporate through superelastic encounters before the cluster reaches the center of the galaxy. This could lead to the formation of an NC with a much smaller abundance of BHs relative to stars than what is predicted by standard initial mass functions. On the other hand, for massive clusters after the BHs are already segregated to the center, the encounter driven evaporation timescale of the BH sub-cluster typically requires an additional few Gyr of evolution to complete (Dowing et al. 2010, 2011). Moreover, recent theoretical studies (Morscher et al. 2013; Sippel & Hurley 2013), together with several observational evidences (Maccarone et al. 2007, 2011; Brassington et al. 2010; Strader et al. 2012), show that old globular clusters may still contain hundreds of stellar BHs at present, which suggests that BH depletion might not be as efficient as previously thought. This indicates that for many large clusters (the ones most relevant to NC formation), most of the BHs will not be ejected before inspiral has occurred. The above arguments convinced us that the inspiral of massive clusters in the central region of the Galaxy could serve as a continuous source term of BHs in these regions.

In Antonini et al. (2012) we used N-body simulations to study how the presence of an MBH at the center of the Milky Way impacts the globular cluster merger hypothesis for the formation of its NC. We determined the properties of the stellar distribution in a galactic nucleus forming through the infall and merging of globular clusters. We showed that a model in which a large fraction of the mass of the Milky Way NC arose from infalling globular clusters is consistent with existing observational constraints. Here we replaced the single-mass globular cluster models of Antonini et al. (2012) with systems containing (in addition to the stellar component) a remnant population of BHs.

These simulations were performed by using ϕGRAPE (Harfst et al. 2006), a direct-summation code optimized for running on GRAPE accelerators (Makino & Taiji 1998). This integrator is equivalent to ϕGRAPEch, which we used in Section 2, but without the regularized chain. The accuracy and performance of the code are set by the time-step parameter η and the softening length . We used a Plummer softening for the gravitational force between particles, and we did not model binary formation in the calculation reported below. We set η = 0.01 and = 10⁻⁴ pc. With these choices, energy conservation was typically ≲ 0.01% during each merging event. The simulations were carried out using the 32 node GRAPE cluster at the Rochester Institute of Technology and also on Tesla C2050 graphics processing units on the Sunnyvale cluster at the Canadian Institute for Theoretical Astrophysics. In the latter integrations, ϕGRAPE was used in serial mode with sapporo (Gaburov et al. 2009). Each simulation required between three and four months total of computational time.

Table 1 summarizes the parameters of the N-body models. We performed four simulations. Runs A1 and A2 are high-resolution simulations (N ∼ 10⁶) that explore the dynamics of one single globular cluster inspiral. In simulations B and C, the total number of N-body particles was greatly reduced in order to more efficiently follow the consecutive inspiral and merger of 12 dynamically evolved clusters. In these latter simulations, each inspiral simulation was started after the stars from the previously disrupted cluster were set to a state of collisionless dynamical equilibrium and the number of BHs in the cluster remnant dropped to <10. This corresponds to a time of 1–3 × 10⁷ yr for each inspiral event to complete, with the longer times corresponding to the earlier infalling clusters. The clusters were initially placed on circular orbits at a distance of 20 pc from the center. The choice of circular orbits was motivated by the well-known effect of orbital circularization due to dynamical friction (e.g., Casertano et al. 1987; Hashimoto et al. 2003). In the consecutive merger simulations (runs B and C), in order not to favor any particular direction for the inspiral, the orbital angular momenta were selected in the following way: the surface of a sphere can be tessellated by means of 12 regular pentagons, the centers of which form a regular dodecahedron inscribed in the sphere. The coordinates of the centers of these pentagons were identified with the tips of the 12 orbital angular momentum vectors. In this way, the inclination and longitude of ascending node of each initial orbit were determined.

We adopted two different N-body models to represent the central region of the Galaxy. Model 1 is obtained by an inner extrapolation of the observed density profile of stars in the Galactic nuclear bulge outside 10 pc. In these regions the Galaxy is dominated by the presence of the nuclear stellar disk, which is characterized by a flat density profile. Accordingly, we adopted the truncated shallow power-law density model:

$$\begin{equation} \rho _{\rm gx,1}(r)=\tilde{\rho }\left(\frac{r}{\tilde{r}}\right)^{-\gamma }{\rm sech}\left(\frac{r}{r_{\rm cut}}\right), \end{equation} \tag{ 9 }$$

where $\tilde{\rho }=400$ M_☉ pc⁻³ is the density at $\tilde{r}=10$ pc, and the truncation radius is r_cut = 20 pc. Hence, the initial conditions for model 1 do not include a preexisting NC, and they correspond to a shallow density cusp around a central MBH.

Model 2 was obtained by the superposition of the density model of Equation (9) and a broken power-law model representing an NC:

$$\begin{eqnarray} \rho _{\rm gx,2}(r)&=&\rho _{\rm gx,1}(r) + \rho _b\left(\frac{r}{r_{b}}\right)^{-\gamma _i} \nonumber\\ && \times \left[1+\left(\frac{r}{r_{b}}\right)^\alpha \right]^{(\gamma _i-\beta)/\alpha }{\rm sech}\left(\frac{r}{r_{\rm cut}}\right), \end{eqnarray} \tag{ 10 }$$

with ρ_b = 4.1 × 10⁴ M_☉ pc⁻³, r_b = 1.5 pc, γ_i = 0.5, β = 1.9, α = 3.73, and r_cut = 20 pc. This model corresponds approximately to the best-fitting density profile of the simulation end-products of Antonini et al. (2012).

In both model 1 and model 2 we included a central MBH of mass M_• = 4 × 10⁶ M_☉ and we generated their N-body representations via numerically calculated distribution functions.

We generated our globular cluster initial conditions following the same procedure described in Antonini et al. (2012), where a detailed description of the initial conditions of the clusters can be found. In brief, the clusters are started on circular orbits of radius r = 20 pc, and their initial masses and radii are set up in such a way as to be consistent with the galactic tidal field at that radius. The clusters are King models (King 1962) with central (King) potential W₀ = 5.8, core radius r_k = 0.5 pc, and central velocity dispersion σ_K = 35 km s⁻¹. With this set of parameters the truncated mass of the clusters was m_cl ≈ 1.1 × 10⁶ M_☉.

To these models we added a heavier mass group representing a population of stellar BHs. The relative values of the particle masses was 1:10. These represent one solar mass stars and 10 M_☉ BHs, respectively. The two mass groups had the same initial phase-space distribution. Standard population synthesis models predict that about 1% of the total mass in a stellar system will be in BHs, top-heavy mass functions result in about five times more BHs. Accordingly, in our initial models the total mass in BHs was 10⁻² and 5 × 10⁻² times 4 × 10⁶ M_☉ (the mass of the non-truncated King model) for runs A and B and run C (Table 1), respectively. The choice of scaling the total number of BHs to the initial non-truncated cluster mass is based on the fact that when the cluster reaches 100 pc (roughly the radius at which our models start to be tidally truncated) mass segregation is likely to have already occurred. under these circumstances, tidal stripping will preferentially remove stars from the outer part of the system, leading to an overabundance of BHs with respect to standard mass functions (Banerjee & Kroupa 2011). In addition, this choice resulted in a sufficiently good statistics for the remnants population. The effect of varying the initial mass in stellar BHs and their initial dynamical state inside the clusters will be investigated in a future paper.

The mass-segregated cluster models were then created via N-body integrations, starting from the cluster equilibrium models. Figure 6 gives the evolution of the two mass components during these integrations. We let the system evolve for a few relaxation times as defined by Equation (7). The stellar BHs accumulate toward the center and by approximately one half-mass relaxation times their distribution appears to have reached an approximately steady state. At the same time the density profile defined by the stellar component undergoes a slow expansion due to heating by the BHs.

Zoom In Zoom Out Reset image size

Figure 7 shows surface density contours of the single inspiral simulation A1. After ∼10⁷ yr the stellar cluster is at about 5 pc from the center; at these distances the disruption process due to tidal stress from the central MBH begins. Rapid removal of stars from the outer part of the cluster by the galaxy and MBH tidal fields unveils its mass-segregated BH cluster. Figure 8 shows the time evolution of radius and bound mass for the globular clusters in runs A1 and A2. Our cluster models rapidly evolve to a state of dark stellar cluster, i.e., a dense cluster dominated by dark stellar remnants (Banerjee & Kroupa 2011). After this state is reached, because of the drop in the total cluster mass (see lower panels of Figure 8), the dynamical friction drag on the remaining BH cluster is largely suppressed, slowing down its orbital decay toward the center of the galaxy. Noticeable, because of the common motion around the system MBH-cluster center of mass, the MBH is significantly displaced from the galaxy center. More precisely, we found a maximum displacement of ∼5 pc in run A1 and a somewhat smaller (maximum) displacement of ∼2 pc in run A2.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 8 shows that after the stellar clusters are disrupted, the remaining dark clusters have a bound mass of ∼2.7 × 10³ M_☉ in simulation A1 and 1.7 × 10⁴ M_☉ in simulation A2, corresponding to 9 and 67 BH particles, respectively. The enhanced removal of stars and BHs decelerates the orbital evolution of the cluster because of its lower mass. During the inspiral after about 2 × 10⁷ yr the BH cluster core has collapsed to ≲ 0.05 pc. This makes the central density much higher, which prevents the complete disruption of the BH cluster.

Assume that the cluster has reached a state of "thermal equilibrium" at the center, i.e., a state in which the stars and BHs are represented by lowered Maxwellians: $m_{\rm st}\sigma _{\rm st}^2=m_{\rm bh}\sigma _{\rm bh}^2$, where σ_st (σ_bh) is the central one-dimensional velocity dispersion of the stars (BHs). If an MBH of mass M_• is present at the center of the galaxy, disruption occurs at distance

$$\begin{equation} r_\mathrm{disr} \approx 2 \left(\frac{\sigma _\mathrm{NC}}{5\sigma _K}\right)^{2/3} \left(\frac{r_\mathrm{infl}}{r_K}\right)^{1/3} r_K \end{equation} \tag{ 11 }$$

from the MBH. Then the BH cluster tidal disruption radius will be smaller than that of the stellar cluster by a factor

$$\begin{equation} \approx 0.3 \left({20} \frac{r_{K; \rm bh}}{r_{K; \rm st}} \right)^{2/3} \left(\frac{1}{10} \frac{m_{\rm bh}}{m_{\rm st}} \right)^{1/3}, \end{equation} \tag{ 12 }$$

which suggests that the BHs can end up being more centrally concentrated than the stars. A condition for this to happen is that the BH cluster must not evaporate before it has lost significant orbital energy by dynamical friction. In fact, the internal evolution of a compact BH cluster embedded in the extreme tidal field of the GC can proceed very rapidly and lead to the cluster complete dissolution on a short timescale, of order a few Myr (e.g., Banerjee & Kroupa 2011; Gürkan & Rasio 2005). In our simulations the internal dynamical evolution of the clusters has been suppressed by giving a non-negligible softening radius to the cluster particles. In fact, the adopted integrator cannot treat the postcollapse evolution of the cluster, since we used a softened potential. Thus, we terminated the simulation at ≈3 × 10⁷ yr of evolution after the clusters have been fully depleted of stars and the remaining BH clusters underwent core collapse.

The end-product spatial density profile and cumulative mass distribution of stars and BHs are given in Figure 9. In order to obtain the BH density profile, we forced the unbinding of the remnant clusters after 10⁷ yr of evolution. The unbinding of the clusters was induced by "turning off" the gravitational interaction terms between the BH cluster members and by letting the system evolve for about one crossing time. Although quite artificial, this procedure allowed us to account for the fact that the dissolution time of the cluster remnants is expected to be short relative to their dynamical friction timescale—in our models a 10⁴ M_☉ system starting from a galactocentric radius of 1 pc reaches a radius of 0.2 pc after 10⁸ yr. The BH clusters will dissolve on a timescale proportional to the half-mass relaxation time, T_ev ≈ 300 × T_h (Spitzer 1987) (if we ignore the effect of the external tidal field); from Equation (6), using r_h = 0.1 pc (Figure 6), N = 1000, and m = 10 M_☉ we find T_ev ≈ 10⁷ yr. We stress that since the current state of the art computational capability does not allow us to simulate the actual number and mass of stars and to calculate the internal evolution of the cluster, the densities of BHs obtained here should be only considered as approximate (likely an underestimate of the real density). Note also that it is unlikely that the core-collapse phase of a BH cluster can lead to the formation of an intermediate mass black hole since any BH–BH merger will eject the remnant from the cluster via asymmetric emission of gravitational wave radiation before it can accrete other BHs or surrounding stars.

Zoom In Zoom Out Reset image size

The number of BHs in our simulations, N_{BH, nb}(< r), was converted to a predicted number for a Milky Way model, using the approximate scaling

$$\begin{equation} N_{\rm bh}(<r)=N_{\rm BH, nb}(<r) \times \frac{4\times 10^6 \, M_{\odot }}{10 \, M_{\odot }} \frac{m_{\rm BH, nb}}{M_\bullet }, \end{equation} \tag{ 13 }$$

with the last factor of the second term containing the masses in the units of the N-body code. The number of BHs transported in the inner ∼1 pc is of order 100 in both A1 and A2. Our simulations result in a mass distribution characterized by a flat density core inside ∼2 pc in run A1 and ∼3 pc in run A2 and an envelope that falls off rapidly at large radii. The BH density distributions flatten within a radius comparable to the size of the core observed in the stellar density profile. This is because the BH cluster does not experience significant orbital decay after the star cluster is fully dispersed. The difference in the mass distribution in the two models A1 and A2 is caused by the difference in the enclosed mass in the background galaxy. The tidal field near the galactic center is much stronger in simulation A2 because of the presence of a preexisting compact stellar nucleus. The stronger tidal field in the galaxy model of simulation A2 results in a larger core in the density distribution due to the larger tidal disruption radius of the star cluster.

The dashed curves in Figure 9 display the density profile of the galaxy background at the end of the simulations. A comparison with the functional forms of Equations (9) and (10) used to generate the initial equilibrium models shows that the background galaxy in A1 did not evolve appreciably; in A2, instead, the final density of the galaxy appears to have slightly changed, showing higher central densities and a smaller core radius (∼0.5 pc) than the initial model. Thus, the influence of the inspirals on the preexisiting distribution of field stars in our models is negligible at large radii ≳ 3 pc, while it leads to slightly higher central densities within this radius and a smaller core radius relative to the initial model distribution (we discuss this point in more detail below in Section 5.3).

The key question to be answered by these simulations is the degree to which the density of stellar BHs near the center of the galaxy is enhanced, after the inspiral, with respect to the relative density expected in the absence of dynamical evolution. Figure 10 shows the radial profile of the ratio, ρ_bh/ρ_st, of BHs to stellar densities. The dotted curve in the figure gives the density ratio, ρ_bh/ρ_st = 0.052, of the initial model before the BHs segregate to the cluster center, which is also the relative BH density expected in the absence of the cluster dynamical evolution. Inside 5 pc, our simulations result in larger BH densities relative to what is expected if the remnant population had the same density distribution of stars at the moment the clusters reach their tidal disruption radius.

Zoom In Zoom Out Reset image size

Figure 10 suggests that the evolution of an NC formed through the merger of dynamically evolved stellar clusters will be dominated by the BHs. The condition that the evolution of the light component is dominated by scattering off the heavy component is (e.g., Gualandris & Merritt 2012)

$$\begin{equation} \rho _{\rm bh}\gtrsim (m_{\rm st}/m_{\rm bh}) \rho _{\rm st}=0.1\rho _{\rm st}. \end{equation} \tag{ 14 }$$

From Figure 10 we see that the evolution of the stars will be dominated by gravitational interactions with the BHs—as opposed to self-interactions—within a radius of size ∼3 pc, roughly the MBH influence radius. We conclude that the post-infall long-term evolution of the systems presented here will be very different from that of the models discussed in Section 2 for which about one Hubble time was required in order to first met the condition Equation (14).

In order to determine the distribution of stars and stellar remnants predicted by a dissipationless formation model for the Milky Way NC, we performed simulations that followed the repeated merger of mass-segregated massive clusters in the GC; these correspond to simulations B and C of Table 1. In these integrations the number of N-body particles was much reduced with respect to the single inspiral simulations presented in Section 5.1. The sub-sampling was a necessary compromise to keep the computational time from becoming excessively long, while still allowing the simulations to follow the successive inspiral of 12 clusters into the center of the galaxy, giving a total accumulated mass in stars of ≈10⁷ M_☉, which is roughly the observed mass of the Milky Way NC.

The process of NC formation is illustrated in Figure 11, which shows the growing central density of stars and BHs during the globular cluster inspirals. As clusters merge to the center, the peak density of the model increases and an NC forms, appearing within the model inner ≈10 pc as an excess density of stars over the background density of the galaxy. Insert panels give the radial dependence of the density profile slope of the final NC models. In the radial range 0.5 pc ≲ r ≲ 5 pc, the spatial density profile of the merger product is characterized by a density of stars that rises steeply toward the center roughly as ∼r^−1.5. At smaller radii at the end of both simulations B and C, the spatial density profile of stars flattens and the radial dependence becomes approximately dlog ρ/dlog r ≈ −1 inside 0.3–0.5 pc. The BH population exhibits a very steep density cusp dlog ρ/dlog r ≈ −2.2, outside ∼0.5 pc, and a somewhat flat profile within this radius. The merger process produces an NC that is in an state of advanced mass segregation, with the heavy component dominating the density of stars inside a radius of roughly 0.3 pc and 1 pc in simulations B and C, respectively. Since smaller systems have shorter relaxation times and undergo mass segregation more quickly, the merger process effectively reduces the mass-segregation timescale of the NC compared for instance with the models discussed in Section 2.

Zoom In Zoom Out Reset image size

The right panels of Figure 11 show the projected density of the N-body model at the end of the simulations. These profiles were fitted as a superposition of two model components, one intended to represent the galaxy and the other the NC. For both components we adopted the Sérsic law profile:

$$\begin{equation} \Sigma (R)=\Sigma _{0} {\rm exp}{\left[ -b\left(\frac{R}{R_{0}} \right)^\frac{1}{n} +b\right]}, \end{equation} \tag{ 15 }$$

with

$$\begin{equation} b=2n - \frac{1}{3} + \frac{0.009876}{n}. \end{equation} \tag{ 16 }$$

For the end-product N-body model of run B, the best-fit parameters were Σ₀ = 7.21 × 10³ M_☉ pc⁻², n = 1.03, and R₀ = 31.9 pc for the bulge and Σ₀ = 1.03 × 10⁴ M_☉ pc⁻², n = 1.88, and R₀ = 9.51 pc for the NC. The best-fit parameters of the merger product of run C were Σ₀ = 7.44 × 10³ M_☉ pc⁻², n = 1.04, and R₀ = 31.6 pc for the bulge and Σ₀ = 7.60 × 10³ M_☉ pc⁻², n = 2.08, and R₀ = 11.3 pc for the NC.

We note that the final density profile and structure of the NC depends on a variety of factors, which include the initial distribution of stars and BHs inside the clusters, the strength of the tidal field due to the galaxy and MBH, and how the distribution of previously migrated stars and BHs evolves in response to their gravitational interaction with other background stars and with infalling clusters. For instance, the initial degree of internal evolution in the cluster models together with the adopted galactic MBH mass will determine how close to the galactic center the BH and stellar clusters will get before they completely dissociate; in turn this regulates the size of the region over which the stellar distribution flattens (i.e., the core size of the NC density profile), as well as the number of stellar BHs transported in the vicinity of the MBH.

In Figure 12 the mass in BHs accumulated in the inner parsec, M_bh(< 1 pc), is shown as a function of the NC mass, M_NC. At any time the NC mass is given by the sum of the accumulated globular cluster masses. A good fit to the data for M_NC > 4 × 10⁶ M_☉ is given by M_bh(< 1 pc) = a × (M_NC/4 × 10⁶ M_☉)^b, with a = 3.33 × 10⁴ M_☉ and b = 1.44 for model B and a = 1.42 × 10⁵ M_☉ and b = 0.819 for model C; these fitting functions are plotted in Figure 12 as solid curves.

Zoom In Zoom Out Reset image size

The effect of the infalling clusters on the preexisting NC is illustrated in the left panel of Figure 13, which shows the density profile of stars in simulation B after 12 inspirals were completed and that were originally inside the 3rd, 6th, 9th, and 12th infalling cluster. The density profile of stars from the 12th merged cluster, consistently with the results of the high-resolution simulations of Figure 9, is characterized by a shallow density cusp out to roughly 1–2 pc and an outer envelope with density that falls off rapidly with radius. The dot–dashed curve in the figure illustrates the density profile of the same stars after the final NC model was run in isolation for a time equal to the time for three consecutive inspiral events to occur, ∼5 × 10⁷ yr. The similarity between the two mass distributions indicates that collisional two-body relaxation, due to random star–star and star–BH gravitational encounters, can be ignored during the inspiral simulations. The densities of stars coming from previously decayed globular clusters, however, appear very similar to each other and very different from the density of stars transported during the last inspiral, being lower and having a steeper radial dependence, dlog ρ/dlog r ≈ −1.5, in the radial range 1 pc ≲ r ≲ 10 pc.

Zoom In Zoom Out Reset image size

As illustrated in the right panel of Figure 13, this evolution was also present in the one-component inspiral simulations of Antonini et al. (2012); the initial conditions of these simulations were essentially the same as those of simulation B with the difference that the clusters only had a single-mass particle group representing stars. There are two mechanisms that drive the evolution of the density profile toward their final form: (1) stars from earlier infalling clusters are stripped at smaller radii and dominate inner regions of the NC (Perets & Mastrobuono-Battisti 2014) and, as argued in the following section, (2) gravitational scattering of the previously accumulated stars and BHs by the infalling clusters.

Consider a massive cluster of mass m_cl that moves into a system of N stars. The second-order diffusion coefficients that appear in the Fokker–Plank equation and that describe the evolution of the stellar distribution due to self-scattering ($\mathcal {D}_{EE,11}$) and scattering off of the cluster ($\mathcal {D}_{EE,1 {\rm cl}}$), scale with the mass of the perturber and the mass (m) and density of field particles as

$$\begin{eqnarray} &&\mathcal {D}_{EE,11} \simeq m^2 N;\quad \mathcal {D}_{EE,1 {\rm cl}} \simeq m_{\rm cl} ^2, \end{eqnarray} \tag{ 17 }$$

where here we have approximated the cluster as a point mass perturber. From Equation (17) we see that if $m_{\rm cl}\gg \sqrt{m^2N}$, $\mathcal {D}_{EE,1 {\rm cl}}\gg \mathcal {D}_{EE,11}$, i.e., self-scattering is negligible compared with scattering off of the massive perturber; the first-order coefficients can also be ignored since they are smaller than $\mathcal {D}_{EE,1 {\rm cl}}$ by factors of $m^2N/m_{\rm cl}^2$ and m/m_cl. Thus, if

$$\begin{equation} m_{\rm cl}> \sqrt{m^2N}, \end{equation} \tag{ 18 }$$

an inspiraling cluster will reduce the energy relaxation timescale compared with that due to stars alone by a factor ${\sim} m_{\rm cl}^2 \ln \Lambda /m^2N \ln \Lambda ^{\prime }$. In the latter expression, Λ' represents the Coulomb logarithm estimated by taking into account the large physical size of the cluster, which implies a lower effectiveness of close gravitational scattering with respect to star–star scattering (e.g., Merritt 2013). This latter quantity can be set to

$$\begin{equation} \ln \Lambda ^{\prime } \approx \frac{1}{2} \ln \left[ 1+ \frac{p_{{\rm max}}^2}{p_0^2}\right] \end{equation} \tag{ 19 }$$

with p_max roughly one-fourth times the linear extent of the test star's orbit. Setting this size to 5 pc, and p₀ to half of the size of the cluster, ≈1 pc, the typical relaxation time becomes

$$\begin{eqnarray} && T^{{\rm eff}}_r(r)\approx \frac{0.34 \sigma ^3(r)}{G^2 \ln {\Lambda ^{\prime }} \rho } \frac{M_{\bullet }}{m_{\rm cl}^2} = 2 \times 10^7 \left(\frac{\sigma }{50 \ {\rm km\ s^{-1}}} \right)^3 \nonumber\\ && \quad \times \left(\frac{1.6}{\ln \Lambda ^{\prime }} \frac{10^3 \, M_{\odot }}{\rho } \frac{M_{\rm \bullet }}{4\times 10^6 \, M_{\odot }}\right) \left(\frac{5\times 10^5 \, M_{\odot }}{m_{\rm cl}}\right)^2 \ {\rm yr}, \quad\quad \end{eqnarray} \tag{ 20 }$$

where in deriving this expression we have used the fact that mN ≈ M_• near the sphere of influence of the MBH. For M_• ∼ 10⁷ M_☉, the infall of even a cluster of 10⁴ M_☉ would largely affect the rate of collisional relaxation, provided that the cluster spends a time of the order of $T^{{\rm eff}}_r$ in a region where the condition (18) is satisfied.

The left panel of Figure 14 shows the total mass of BHs and stars bound to the 12th infalling cluster in run B as a function of cluster orbital radius and compares these with $\sqrt{m^2N}$. The cluster spends roughly 5 × 10⁶ yr before it reaches ∼2 pc (right panel of Figure 14), after which its mass drops below $\sqrt{M_{\bullet } m}$ and self-scattering starts to be the dominant effect driving collisional relaxation. Since the time for the cluster to reach this radius is comparable to $T^{{\rm eff}}_r$, the inspiral of clusters is expected to have an important impact on the density distribution of the preexisting NC, in agreement with the results of our simulations. In the right panel of Figure 14 we show the evolution of the Lagrangian radii of the stars transported in the NC during the previous (11th) inspiral. During the first 10⁷ yr, the condition (18) is satisfied and the inspiral induces a rapid evolution of the NC mass distribution. Note that in a real galaxy, because of the smaller individual stellar mass, scattering from the perturber would dominate down to smaller radii than in our N-body model. However, in the region where the condition (18) is satisfied, the relaxation time in the simulations would be roughly the same as in the real system.

Zoom In Zoom Out Reset image size

Scattering of stars by the massive perturber will cause an initial density core to fill up and the distribution function to evolve toward a constant value. The result is a sharp increase in the density profile of stars, n ∼ r^−1.5, inside a radius approximately equal to the radius within which the condition Equation (18) is first met. After the perturber reaches a radius r_crit containing a total mass in stars smaller than its mass, a large density core is rapidly carved out as a consequence of ejection of stars from the center. This is the evolution that, for example, characterizes the mass distribution of stars during the formation and evolution of MBH binaries in galactic nuclei (see Figures 13 and 14 of Antonini & Merritt 2012). In our simulations, however, the clusters disrupt before reaching r_crit and before the second phase of cusp disruption can initiate. Thus, the net effect of infalling clusters on the preexisting NC distribution is that of inducing a short period of enhanced collisional relaxation, which causes the central density of stars and BHs to increase during the inspirals.

During and after its formation, an NC will evolve because of collisional star–star, star–BH, and BH–BH interactions that will cause its density distribution to slowly morph into the steeply rising density profile, which describes the quasi-steady-state solution of stars and BHs near an MBH. We study the evolution of the NC due to collisional relaxation by evolving the inspiral simulation end-products of simulation B for roughly half a relaxation time (or one Hubble time when scaling the N-body model to the Milky Way). In order to efficiently evolve the system for such a long timescale, the N-body model was resampled to contain a smaller number of particles. Since we are interested in the model distribution at small radii, we kept the mass of the particles the same as in the original model, so that the resolution of the simulation in the region near the center was unchanged, and we included only particles with orbital periapsis less than 10 pc and apoapsis less than 20 pc. In this way the total gravitational force acting on particles lying inside ≈10 pc was approximately unchanged with respect to the original model. This sampling procedure is equivalent to truncating the mass distribution of the model smoothly at r ≳ 10 pc. N-body integrations over a few crossing times verified the (collisionless) quasi-equilibrium state of the truncated model.

Since the dynamical friction times for globular clusters at 20 pc from the galactic center are much shorter than relaxation times, both the cluster and galaxy will not undergo a significant amount of collisional relaxation during the inspirals. In our N-body simulations B and C, the inspiral time of the globulars to reach the center and disperse around the central MBH is about 100 times shorter than the relaxation time of the NC. Thus, during the inspiral simulations, relaxation due to star–star, star–BH, and BH–BH interactions was ignored, while in the post-merger simulations presented here times were scaled to the relaxation time computed at the sphere of influence of the MBH. Note that for the sake of simplicity, the evolution was broken into two successive stages: the infall of the clusters and then evolution, due to two-body encounters, of the stellar distribution around the MBH with infalls "turned off." In reality, subsequent inspirals would be separated by times of order a Gyr and significant two-body relaxation would occur between these events.

Figure 15 shows the evolution of the mass density in each species in the post-merger integrations. Initially, the density of BHs dominates over the density in stars at r ≲ 0.1r_infl. Thus, at radii smaller than these, the evolution of the BHs is mostly driven by BH–BH self-scattering that causes their distribution to reach, in roughly $0.2\times T_{r_{\rm infl}}$, a quasi-steady-state form characterized by a steep density slope, γ ≈ 2.2, at all radii. At radii larger than ∼0.1r_infl, the stars dominate the mass density and the BH population evolves because of dynamical friction against the lighter component. The general trend is for the central mass density of BHs to steadily increase with time.

Zoom In Zoom Out Reset image size

While the BHs segregate to the center, the lighter particle mass density decreases within the radial range 0.1 ≲ r ≲ 10 pc. The evolution of the stellar distribution toward lower densities is a consequence of their gravitational interaction with the BHs: when the density in the BHs approaches locally the density in stars, heating of the light particles off the heavy particles dominates over star–star scattering, causing the density of the former to decrease. Scattering of stars off the BHs is also expected to promptly modify the star density profile at small radii, causing the formation of a "mini-cusp" at r ≲ 0.02r_infl (Gualandris & Merritt 2012). We find evidence of this phenomenon in our N-body models, which rapidly develop a n ∼ r^−1.5 cusp at r < 0.03 pc (see upper panel and corresponding insert panel of Figure 15). Outside these radii, however, even after about half a relaxation time the star distribution retains the shallow density profile at r ≲ 0.2 pc that characterized the initial model. We conclude that even after a time of order the relaxation time, model B looks different from the dynamically relaxed models that are often assumed to describe the density distribution of stars and BHs near the center of galaxies. After ${\approx }0.4\times T_{r_{\rm infl}}$, the BHs (stars) attain a central density cusp that is steeper (shallower) than in those models.

We note that the details of the final density distribution of our NC model depend on a number of factors that remain quite uncertain. For example, the final state of the BH population in the NC and the degree at which the BHs are segregated in the final model will depend on their initial number fraction and thus on the assumptions made for their initial distribution in the parent clusters.

Understanding the kinematical structure and shape of galactic nuclei is important for placing constraints on their formation history and evolution. We quantified the velocity anisotropy of our models by using the parameter

$$\begin{equation} \beta =1-\frac{\sigma _t^2}{2\sigma _r^2}, \end{equation} \tag{ 21 }$$

where σ_r and σ_t are the radial and tangential velocity dispersions, respectively. Figure 16 shows the radial profile of β, σ_r, and σ_t for the stellar and BH populations in simulations B and C. At the end of the inspiral simulations, model B is characterized by an approximately flat velocity anisotropy profile with β ≈ −0.5, while in model C β slightly decreases from nearly 0 to −0.2 within 1 pc and is approximately constant outside this radius. Thus, our models are tangentially anisotropic throughout 100 pc, both in the stellar and BH components. The upper panels of Figure 16 also display the radial profile of the anisotropy parameter of model B at the end of the post-merger phase, i.e., after the system was run in isolation for a time of order the relaxation time. Evidently, two-body relaxation causes β to increase and the NC to evolve toward spherical symmetry in velocity space.

Zoom In Zoom Out Reset image size

We measured the model shape in our simulations from the moment-of-inertia tensor (e.g., Katz 1991; Poon & Merritt 2004; Antonini et al. 2009). The symmetry axes are calculated as

$$\begin{equation} \tau _1=\sqrt{I_{11}/I_{{\rm max}}},\quad \tau _2=\sqrt{I_{22}/I_{{\rm max}}},\quad \tau _3=\sqrt{I_{33}/I_{{\rm max}}}, \end{equation} \tag{ 22 }$$

where I_ii is the principal moments of the inertia tensor and I_max = max{I₁₁, I₂₂, I₃₃}; particles are then enclosed within the ellipsoid x²/τ₁² + y²/τ₂² + z²/τ₃² = r². These previous steps were iterated until the values of the axial ratios had a percentage change of less than 10⁻². We also define a triaxiality parameter: T ≡ (a² − b²)/(a² − c²). The value T = 0.5 corresponds to the 'maximally triaxiality'case, while oblate and prolate shapes correspond to T = 0 and 1, respectively.

Figure 17 displays the radial profile of the axis ratios and triaxiality parameter of the model in simulation B at the end of the inspiral simulation (black solid curves) and at the end of the post-infall evolution (blue dotted curves). We also evaluated the shape of the NC component by only including the stars that were transported to the center from the infalling stellar clusters (right panels); the NC component appears strongly triaxial-like within 5 pc and mildly triaxial (quasi-oblate) outside this radius. The formed NCs in simulations C and B shared a similar morphological structure, so we only displayed results for model B. The model morphology within ∼30 pc is initially mildly triaxial and evolves because of gravitational encounters toward a more quasi-spherical shape. It is important that at the end of the post-merger phase the model still exhibits a significant degree of triaxiality, 0.1 ≲ T ≲ 0.3. In fact, such a level of asymmetry might be large enough to significantly increase the number of tidal captures of stars and stellar binaries when compared with the same rate obtained in collisionally resupplied loss cone theories where spherical geometry is often assumed both in configuration and velocity space (Merritt & Vasiliev 2011).

Zoom In Zoom Out Reset image size

Star counts using adaptive optics spectroscopy and medium-band imaging have shown that the red giants at the Galactic center, the only old stars that can be resolved in these regions, have a flat projected surface density profile close to Sgr A* (Buchholz et al. 2009; Do et al. 2009). The core in the red giant population extends out to approximately 0.3–0.5 pc from the center. However, because of the effect of projection, it is difficult to constrain the core size and three-dimensional spatial density profile, which could be slowly rising or declining toward the center.

It is possible that a cusp in the lower mass stars is present and that the observed core is the result of a luminosity function that changes within 0.5 pc. This could be due to physical collisions before or during the giant phase (Bailey & Davies 1999; Dale et al. 2009) or tidal interactions between stars and the central MBH (Davies & King 2005)—mass removal can make the luminosity that a star would otherwise reach at the tip of the red-giant phase considerably fainter and prevent these objects from evolving to become observable. While these models are possible, they seem to not fully explain the observations (Dale et al. 2009). Thus, it is important to consider the possibility that the red giants are indeed representative of the unresolved low-mass main-sequence stars and that the density core is a consequence of the NC formation history combined with a long nuclear relaxation timescale (Section 2). On this basis we can directly compare the predictions of theoretical models with the kinematics and mass distribution inferred from observations of the giant stars at the GC and draw conclusions about possible formation mechanisms for the Milky Way NC.

6.1.1. Mass Distribution

6.1.2. Kinematics

6.1.3. Kinematically Cold Sub-structures

Feldmeier et al. (2014) found indications for a sub-structure in the Galactic NC that is rotating approximately perpendicular to the Galactic rotation with ∼30 km s⁻¹ at a distance of ∼20'' or 0.8 pc from Sgr A*. In addition they found an offset of the rotation axis from the photometric minor axis and argue that this hints to infalling clusters.

We look for kinematic sub-structures in our models by using the Rayleigh (dipole; Rayleigh 1919) statistics R' defined as the length of the resultant of the unit vectors l_i, i = 1, ..., N, where l_i is perpendicular to the orbital plane of the ith particle and N is the number of particles (a total of 5715 per cluster). For each merged cluster, we computed R' over the entire course of the inspiral simulation. For a fully isotropic distribution, we expect $R^{\prime } \sim \sqrt{N}$, while R' ∼ N if the orbits are correlated.

As an example, Figure 18 gives R(= R'/N) as a function of time (scaled to the relaxation time of a Milky Way–like nucleus). Initially, after a cluster reaches the center the orbits are strongly correlated and, as expected, R ∼ 1, i.e., the stars from the infalling clusters distribute into a thin disk configuration initially. Because of two-body relaxation, and because of the perturbing effect of the later infalling clusters, R decreases with time and approaches a value more consistent with isotropy. Figure 18 shows that all clusters maintain some degree of anisotropy during the entire course of the simulation. However, the orbits of the stars from the first seven clusters are almost completely isotropic by the end of the simulation. The orbits of stars from the last four/five infalling clusters are still largely correlated after ∼3 Gyr of evolution. A linear fit to the data (R versus t) for the last decayed cluster gives R(t) = −0.11 × t/Gyr + 0.96, so that it would take about ∼10 Gyr for the stars to achieve a nearly isotropic distribution. These results are consistent with those of Mastrobuono-Battisti & Perets (2013), who found that it takes a time of order the relaxation time of the nucleus to fully randomize an initially cold disk.

Zoom In Zoom Out Reset image size

The inspiral of compact remnants into an MBH represents one of the most promising sources of gravitational wave radiation detectable by space-based laser interferometers (Amaro-Seoane et al. 2012). Event rates for such extreme mass ratio inspirals are generally estimated under the assumption that the BHs had enough time to segregate and form a steep central cusp, n ∼ r⁻², near the MBH. Such dynamical models predict inspiral rates per galaxies of ∼250 Gyr⁻¹ (Hopman & Alexander 2006). Models that include an initial parsec-scale core can result in much lower central BH densities than in the steady-state models and imply rates as low as 1–10 Gyr⁻¹ (Merritt 2010). EMRI event rates could be also severely suppressed by the Schwarzschild barrier, which limits the ability of stars to diffuse to high eccentricities onto inspiral orbits (Merritt et al. 2011).

It must be stressed that it is difficult to draw any conclusion about EMRI event rates from the models discussed in the literature because of the significant uncertainties in the underlining assumptions. Our simulations cast further doubts on results obtained from idealized time-dependent models, which rely on the assumption that BHs and stars have initially the same spatial distribution. For example, on the basis of the models discussed in Section 2 and in Merritt (2010) and Antonini & Merritt (2012), the presence of a core in the old population of stars at the GC would imply very low central densities of BHs and EMRI event rates. However, the initial conditions adopted in the simulation of Section 2 were quite artificial and not motivated by any specific physical model. We have shown that in a merger model for the formation of NCs, the resulting distribution of stellar remnants partially reflects their distribution in their parent clusters just before they reach the center of the galaxy. Thus, different, possibly more realistic, initial conditions would produce rather different central BH densities; in these models EMRI rates could be as large as (or higher than) those obtained in the steady-state models even in the presence of a core in the stellar distribution.

In order to determine the number of BHs in a Milky Way–like nucleus predicted by the inspiral simulations, we used the scaling relation given by Equation (13). We are interested in the number of compact remnants inside ≈0.01 pc, as these are the only BHs that can generate EMRIs. Since at these radii the number of BH particles in the N-body model is small, we carried out regression fits of log N_bh to log r at larger radii and extrapolate inward in order to get N_bh at the radii of interest. Following Gualandris & Merritt (2012), we performed these fits in the radial interval (r₁, r₂) such that N_bh(< r₁) = 5 and N_bh(< r₂) = 25, and assume a density profile of constant power-law index. In Figure 19 we show the (scaled) cumulative number of BHs during the post-merger simulations of Section 5.4. We find N_bh(< 0.01 pc) ≈ 200 at $t= 0.1\times T_{r_{\rm infl}}$ and N_bh(< 0.01 pc) ≈ 400 at $t= 0.3\times T_{r_{\rm infl}}$. We can directly compare the numbers so obtained with estimates from quasi-steady-state Fokker–Plank models of the GC. In their models, Hopman & Alexander (2006) assumed f = 0.01 and found N_bh(r < 0.01 pc) = 150; Freitag et al. (2006) found similar values. This is similar to the number of BHs in our model at the end of the inspiral simulations but somewhat smaller than the predicted number of BHs at the end of the post-merger simulations.

Zoom In Zoom Out Reset image size

In conclusion, a core in the density distribution of stars does not necessarily imply a low density of stellar remnants in the GC. In a merger model for NCs, BHs initially have a smaller core as a relic of their premerger mass segregation, which was not fully erased by the tidal disruption process of the cluster. After a small fraction of the relaxation time (${\sim }0.1\,T_{\rm r_{\rm infl}}$), the BHs dynamically relax because of BH–BH interactions and attain a steep central density cusp, while the core in the stellar distribution persists. After this time, the number of BHs inside 0.01 pc, the radii relevant for EMRIs, is N_bh(< 0.01 pc) ∼ 100, a value comparable to that inferred in steady-state models.

Observations of the Galactic center have revealed the existence of ∼20 young (B-type) stars orbiting within ∼0.05 pc of the central black hole, the S stars (Ghez et al. 2008; Gillessen 2009). The S stars are thought not to have formed in situ where the strong tidal field of the MBH would prevent star formation (Morris 1993). Mechanisms that invoke the formation of the S stars farther from the MBH and subsequent rapid migration to their current location are often invoked in order to explain their existence. The current paradigm is that the S stars result from the capture of individual stars by the tidal disruption of binaries following a close encounter to the central black hole (Hills 1988). This formation model can successfully account for the number of observed S stars (Perets et al. 2007), but it results in an orbital distribution that is largely biased toward high eccentricities, and therefore it is inconsistent with the quasi-thermal eccentricity distribution of the S star orbits (e.g., Antonini et al. 2010). Post-migration dynamical evolution due to gravitational perturbations from a field population of BHs has been invoked in order to bring the predicted orbital distributions more in line with observations (Perets et al. 2009; Madigan et al. 2011; Hamers et al. 2014). It has been recently shown that such a scenario could indeed lead to an orbital distribution that reproduces the random and eccentric character of the observed orbits (Antonini & Merritt 2013).

The timescale over which the S-star orbits evolve and approach the observed distribution depends on the number of objects, i.e., BHs, in the stellar cusp, in which the larger the number of BHs is, the shorter the timescale is (Antonini & Merritt 2013). Accordingly, we determine a lower limit to the number of BHs in the GC by following the dynamical evolution of the S stars for different values of N_bh and by requiring the orbits to approach the observed distribution in ∼100 Myr, the main-sequence life time of a B-type star.

We ignore two-body relaxation since the timescale of interest (∼100 Myr) is short compared with two-body (non-resonant) relaxation times near the center of the Milky Way (Merritt 2010). Following Antonini & Merritt (2013) we assume that the semi-major axis distribution, N(a), is known, and it is given by the observed values of a. The initial orbital eccentricities were assigned randomly in the range 0.93 ⩽ e ⩽ 0.99. The orbits were initialized at random times between 0 and 100 Myr and followed to a final time of 100 Myr. This setup approximately reproduces the initial conditions expected in a binary disruption scenario for the formation of the S-stars (Zhang et al. 2013). For each S star (i.e., for each value of a), 100 integrations were carried out using the same Hamiltonian model adopted in Merritt et al. (2011) and Antonini & Merritt (2013). Briefly, the Hamiltonian model accounts for the torques due to finite-N asymmetries in the field-star distribution (resonant relaxation) and 1-PN terms that result in the (Schwarzschild) precession of the argument of periastron. The direction of the torquing field was changed smoothly with time and was randomized in a time of order the precession time of the field population as described in Section VB of Merritt et al. (2011).

We consider a field population of 10 M_☉ BHs, which we distributed according to a power-law density profileρ ∼ r⁻². The number of BHs at radii less than r is

$$\begin{equation} N_{\rm bh}(<r)=N_{0.1} \left(r/0.1 \ {\rm pc}\right), \end{equation} \tag{ 23 }$$

where N_0.1 is the number of BHs within a radius of 0.1 pc; we take N_0.1 = (100, 250, 500, 700, 1000), and for each of these values we compared the results of our simulations with observations by performing a two-sample Kolmogorov–Smirnov (K-S) test on the eccentricity distributions.

In any of our cusp models, most orbits have orbital eccentricities that lie initially inside (i.e., at higher eccentricity than) the Schwarzschild barrier—the locus of points in the (E, L) plane where resonant relaxation is suppressed by fast relativistic precession of a test particle orbit (Merritt et al. 2011). As shown in Antonini & Merritt (2013), the Schwarzschild barrier behaves as a one-way permeable membrane, which is penetrable by orbits that approach it from higher eccentricities, while it is a hard barrier for orbits that approach it from above. It can be shown that the time to diffuse above the Schwarzschild barrier, toward higher angular momenta, scales approximately with the number of perturbers as $T_{d}\sim N^{-3/2}_{\rm bh}(<a)$ or as ${\sim} N^{-2}_{\rm bh}(<a)$ if the choice for the coherence time is the "vector" resonant relaxation time or the mass precession time, respectively—i.e., the test particles diffuse faster toward larger angular momenta, above the Schwarzschild barrier, the larger the number of background perturbers is (Figure 20).

Zoom In Zoom Out Reset image size

Figure 21 shows the cumulative eccentricity distributions at 5 × 10⁷ yr and at 10⁸ yr of evolution for various values of N_0.1. As the number of BHs in the cusp increases, the diffusion timescale decreases, allowing the distribution to more rapidly approach a quasi-thermal form. The insert panels give the p-value of the K-S tests as a function of N_0.1. Only for N_0.1 ≳ 1000 after ∼100 Myr does the eccentricity distribution appear to have approached a form that is consistent with observations (p-value ≳ 0.1). These results imply that in order to explain the quasi-thermal character of the S-star orbital eccentricities, the number of BHs within ∼0.1 pc of Sgr A* has to be ≳ 1000. Interestingly, this number is in agreement with the predictions of our globular cluster merger model (Figure 15), but it appears to be much larger than the number of BHs implied by the dynamical models discussed in Section 2, which for f_bh = 10⁻³ give N_0.1 ≈ 100 after 10 Gyr of evolution.

Zoom In Zoom Out Reset image size

Interestingly, a globular cluster merger model for the Milky Way NC can potentially reconcile models that require high central densities of BHs in order to explain the orbits of the S stars with the "missing cusp" problem of the giant star population.