Observational signatures of the black hole Mass Distribution in the Galactic Center

TL;DR

This study simulates the star cluster around SgrA* to analyze the black hole mass distribution, its impact on galactic center dynamics, and the potential gravitational wave signals detectable by LISA.

Contribution

It introduces a detailed simulation framework for the black hole mass function in the Galactic Center and links it to observable signatures and gravitational wave event rates.

Findings

Black hole mass function influences SgrA* growth and dynamics.

Estimated LISA detection rate of inspiraling stellar-mass black holes is 10^-5 per year.

Massive initial black hole populations dominate inspiral event rates.

Abstract

We simulate the star cluster, made of stars in the main sequence and different black hole (BH) remnants, around SgrA* at the center of the Milky Way galaxy. Tracking stellar evolution, we find the BH remnant masses and construct the BH mass function. We sample 4 BH species and consider the impact of the mass-function in the dynamical evolution of system. Starting from an initial 6 dimensional family of parameters and using an MCMC approach, we find the best fits to various parameters of model by directly comparing the results of the simulations after $t = 10.5$ Gyrs with current observations of the stellar surface density, stellar mass profile and the mass of SgrA*. Using these parameters, we study the dynamical evolution of system in detail. We also explore the mass-growth of SgrA* due to tidally disrupted stars and swallowed BHs. We show that the consumed mass is dominated for the BH component with larger initial normalization as given by the BH mass-function. Assuming that about 10% of the tidally disrupted stars contribute in the growth of SgrA* mass, stars make up the second dominant effect in enhancing the mass of SgrA*. We consider the detectability of the GW signal from inspiralling stellar mass BHs around SgrA* with LISA. Computing the fraction of the lifetime of every BH species in the LISA band, with signal to noise ratio $\gtrsim 8$, to their entire lifetime, and rescaling this number with the total number of BHs in the system, we find that the total expected rate of inspirals per Milky-Way sized galaxy per year is $10^{-5}$. Quite interestingly, the rate is dominated for the BH component with larger initial normalization as dictated by the BH mass-function. We interpret it as the second signature of the BH mass-function.

Tables

Table 1. Table 1 : Map between the BH mass and its fractional normalization in density/Mass compared with the ZAMS normalization

Mass($M_{\odot }$)[i]	F[i]
8.15	6.95 $\times {10}^{-3}$
16.24	3.11 $\times {10}^{-2}$
23.92	5.01 $\times {10}^{-2}$
34.50	7.62 $\times {10}^{-3}$

Table 2. Table 2 : Critical value of semi-major axes above which we do not have any inspiral phase.

$m_{BH}\left(M_{\odot }\right)$	$a_{crit}\left(pc\right)$
$8.15$	$0.01$
$16.24$	$0.029$
$23.92$	$0.04$
$34.50$	$0.05$

Table 3. Table 3 : Expected rate of the inspiral as well as the number of LISA sources per one MW like galaxy.

$m_{BH}\left(M_{\odot }\right)$	$R_{obs}\left(y{r}^{-1}\right)$	$N_{lisa}\left(T_{obs}=10yr\right)$
$8.15$	$2.0\times {10}^{-6}$	$2.0\times {10}^{-5}$
$16.24$	$1.54\times {10}^{-5}$	$1.54\times {10}^{-4}$
$23.92$	$2.9\times {10}^{-5}$	$2.9\times {10}^{-4}$
$34.50$	$5.6\times {10}^{-6}$	$5.6\times {10}^{-5}$

Equations

$$\frac{\partial f_c(h,t)}{\partial t}=-\frac{\partial {\text{F}}_c(h,t)}{\partial h}-{\nu }_c(h,t)f_c(h,t)+S(t),$$

$$h(E)=16{\pi }^2{\int }_{\Phi (0)}^{E}d{E}^{\prime }{\int }_0^{r_{max}(E^{\prime })}{r}^2\sqrt{2\left(E^{\prime }-\Phi (r)\right)}dr,$$

$$\Phi (r)$$

$${\text{F}}_c(h,t)=A_c{f}_c+D\frac{\partial f_c}{\partial h},$$

$$A_c(h)$$

$$D(h)$$

$$\nu =\frac{\mu (E)}{\alpha +\ln \left(1/{\mathcal{R}}_{LC}\right)},\alpha \simeq {\left(q^2+q^4\right)}^{1/4},$$

$$q(E)\equiv \frac{\mu (E)P(E)}{{\mathcal{R}}_{LC}(E)},{\mathcal{R}}_{LC}(E)\equiv \frac{L_{LC}^2}{L_{cir}^2}=\frac{2G{M}_{\bullet }{r}_{LC}}{L_{cir}^2(E)},$$

$$\mu (E)\equiv \frac{8{\pi }^2}{g(h(E))}{\int }_0^{r_{max}}\frac{⟨\Delta v_{\perp }^2⟩r^{2}dr}{\sqrt{2(E-\Phi (r))}},$$

$$⟨\Delta v_{\perp }^2⟩$$

$$\rho (r)={\rho }_0{\left(\frac{r}{R_{scale}}\right)}^{-\gamma }{\left[1+{\left(\frac{r}{R_{scale}}\right)}^{\alpha }\right]}^{(\gamma -\beta )/\alpha }.$$

$${\rho }_{BH,0}$$

$${\rho }_{BH,0}^{(i)}={\rho }_{BH,0}\frac{{\int }_{m_i}^{m_{i+1}}{\Phi }_{BH}(m)mdm}{{\int }_{m_{min,BH}}^{m_{max,BH}}{\Phi }_{BH}(m)mdm}.$$

$$\gamma$$

$${\log }_{10}(R_{\mathrm{scale}}[pc])$$

$${\log }_{10}(M_{\ast ,\mathrm{init}})$$

$${\log }_{10}(M_{\bullet ,\mathrm{init}})$$

$${\mathrm{frac}}_{\mathrm{source}}$$

$${\log }_{10}(r_{\mathrm{source}})$$

$$\gamma$$

$${\log }_{10}(R_{\mathrm{scale}}[pc])$$

$${\log }_{10}(M_{\ast ,\mathrm{init}})$$

$${\log }_{10}(M_{\bullet ,\mathrm{init}})$$

$${\mathrm{frac}}_{\mathrm{source}}$$

$${\log }_{10}(r_{\mathrm{source}})$$

$$T_{rel}(r)=\frac{{\sigma }^2}{D[{\left(\Delta V_{||}\right)}^2]},$$

$$D[{\left(\Delta V_{||}\right)}^2]\simeq 2.94\frac{G^2\rho m\ln \Lambda }{\sigma },$$

$$⟨m_{BH}⟩(r)=\frac{{\sum }_i{m}_{BH,i}{\rho }_{BH}^i(r)}{{\sum }_i{\rho }_{BH}^i(r)}.$$

$${\sigma }_i^2(r)=\frac{G}{{\rho }_i(r)}\times {\int }_r^{\infty }\frac{{\rho }_i(r^{\prime })M(r^{\prime })}{r^{\prime 2}}d{r}^{\prime },i=(star,BH).$$

$$\frac{d}{dt}\left(a(t)\right)$$

$$\frac{d}{dt}\left(e^2(t)\right)$$

$${\left(SNR\right)}^2=2\sum _{n=1}^{n_{max}}\int \frac{h_{c,n}^2}{f_{n}S(f_n)}d\ln f_n,$$

$${\left(SNR\right)}^2=\left(\frac{512}{5d^2}\right)T_{obs}\frac{{\left(G{M}_c\right)}^{10/3}}{c^8}{\left(2\pi f_{orb}\right)}^{4/3}\sum _{n=1}^{n_{max}}\frac{g_n(e)}{n^{2}S(f_n)}.$$

$$W(i)=\left(\frac{\Delta t_i}{t_i}\right)\times {\mu }_i,$$

$$W_{tot}(i)=\left(\frac{N_{BH}(r\le r_{sam})}{N_{tot}}\right)\times W(i),$$

$$R_{obs}$$

Full text

1 Introduction

The Supermassive black hole (SMBH) at the Galactic center, Sgr A*, is surrounded by a cluster of stars and their remnant black holes (BHs), which migrate inwards by dynamical friction. An exhaustive consideration of this process was performed in Refs. [1, 2, 3] which considered the formation of a steep density cusp with high velocity dispersion around the galactic center. Mass segregation of heavier objects concentrates them around the SMBH. This generates gravitational waves (GWs) which are observable at cosmological distances [4, 3, 5, 6, 7, 8]. In addition, the close encounters between the stellar mass BHs can also create GWs sources [10, 11, 12, 13, 14, 15, 9]. Another observational signature of dense stellar cluster around galactic centers is X-ray binaries due to the tidal capture of stars by BHs and neutron stars (NSs). This was recently investigated in Ref. [16] and applied to the $\gamma$ ray excess at the center of the Milky Way (MW) [2, 17, 18, 19, 20, 16]. A further interesting observable may originate from mergers of BHs and neutron stars at the Galactic center [21, 22, 23].

There are various ways to model the behavior of such a system. N-body and Monte Carlo codes are computationally very expensive. Using the distribution function is favoured computationally for describing different elements and finding the time evolution of the entire system at different time snapshots. The evolution of the system is given by the Fokker Planck equation, coupled to the Poisson equation. In Ref. [24], the Fokker Planck approach was carefully studied confirming explicitly that the results from this approach are well matched with the results from numerical simulations.

Here we simulate the star cluster at the center of the Milky Way galaxy using the Fokker Planck approach presented in the Phase Flow, hereafter (PF), code which is a library in the publicly available code AGAMA 111https://github.com/GalacticDynamics-Oxford/Agama [24]. Our system is made of a SMBH, Sgr A*, main sequence stars and four species of the stellar-mass BHs as the remnant to the main sequence stars. For the first time, in this context, we use the publicly available tool COSMIC 222https://github.com/COSMIC-PopSynth/COSMIC (Compact Object Synthesis and Monte Carlo Investigation Code) to find the BH remnant mass as a function of the zero age main sequence (ZAMS) mass focusing only on $Z=0.001$. Work is in progress to consider the impact of different metallicities in the analysis as well as enlarging the populations of stellar-mass BHs to more BHs. We consider the continuous star formation in our setup as well. We study a six dimensional parameter space, as our initial condition, and evolve the system for a total period of $t=10.5$ Gyrs. We find the best fit to these parameters by directly comparing the results after $t=10.5$ Gyrs with the most recent observations. Using the best fit values, we study the evolution of the stellar cluster in details. For the first time, we directly consider the impact of the BH mass-function in the simulation of PF. One of the most important places that the impact of BH mass-function directly appears is in the consumed mass to SgrA* from different BH species where we find that the contribution from the second heaviest BH is dominated over the rest of the BHs as well as the stars in the main sequence. Assuming that 10% of the disrupted stars contribute in the growth of the central BH mass, we find that stars are the second most dominant contributor in the mass growth of SgrA*.

As another application of this approach, we study the GW from the inspiral phase of stellar mass BHs around SgrA*. Starting from a sample of 10000 systems, we consider dynamical evolution of the semi-major axes as well as the eccentricity including the impact of the angular momentum diffusion. We compute the fraction of the lifetime of every systems in the LISA band with the signal to noise ratio of 8 or above to their actual lifetime and use this as a weighting factor in our estimation. Furthermore, we also rescale this with the actual number of the BHs interior to $r=10pc$ as the upper limit of our sampling. We compute the expected rate of inspiral for the Milky Way (MW) like galaxies. Our results are in good agreement with the inferred inspiral rate computed from totally other methods.

Our simulations are based on various assumptions and simplifications as listed below:

$(1)$ We assume spherical symmetry for the MW core [29].

$(2)$ We neglect gas inflow in our simulation, the modification of stellar orbits through hydrodynamical drag [31] and through non-spherical contributions (e.g. disk like) to the global gravitational potential [32].

$(3)$ We neglect resonant relaxation. This includes both of scalar and vector resonances as considered in some details at [37, 38, 39, 40]. The importance of scalar resonant relaxation (RR) is still under debate, with controversy in the literature about its importance. While Ref. [33] pointed out that the scalar RR affects the tightly bound stars and enhances the tidal disruption rates, the detail consideration in Ref. [34] showed that relativistic precession suppresses RR and so it is not very efficient. This conclusion is consistent with Ref. [35] where the authors showed that the effect of scalar RR is small (see Fig. 17 of this paper). Vector RR, on the other hand, is harder to predict and it is not obvious whether it might have important effects on the evolution of the system. If important, it would affect the spherical symmetry of the system [36].

$(4)$ Another neglected effect is close binaries near SMBH. As discussed briefly in [55] the impact of binaries on the cusp dynamics is small. As the vast majority of these binaries are soft the total transferred energy by them to the cusp is small. Dynamical evolution of binaries interacting with the stars around SMBH was further performed in Ref. [41] where the author showed that most of the binaries are disrupted either by the loss-cone effects, happening far away from the SMBH, or through the three body encounters, very close to the center. The main conclusion was that binaries are disrupted before they could experience an exchange interaction. Although very close to the GC secular Kozai evolution affects the dynamics of binaries and periodically change their inclinations and eccentricities, these populations are continuously destroyed by different processes such as the evaporation, mergers and disruption [42, 43, 44, 45, 46]. An extensive study of the binary population around SMBH including all of these effects is beyond the scope of this paper and is left for a future work.

The outline of our paper is as follows. In section 2 we review the Fokker Planck approach. In section 3 we describe the initial conditions of the simulation. Using an MCMC approach, we find the best fit for a 6 dimensional family of parameters and use them as our initial conditions throughout our simulations. In section 4 we study the two body relaxation in detail and consider the time evolution of Sgr A*, as well as the mass profile of the stars and BHs. We also compute the mass-radius relation for the BHs at a few different times and initial slopes. In section 5, we consider the detectability of this system with LISA . We conclude in section 6.

2 Fokker-Planck approach for the relaxation process

Here we describe our simulation setup. We use the PF code which is based on a Fokker Planck approach. We compute the time evolution of the distribution function of a muti-body system made of the stars on the main sequence and different BH species being influenced by a central super-massive BH.

[TABLE]

here $f_{c}$ refers to the phase space distribution function of every species, with sub-index $c$ referring to stars and BHs. Eq. (2.1) contains a few additional parameters, $h(E),\textit{F}_{c}(h,t)$, $\nu_{c}(h,t)$ and $S(t)$.

In the following, we provide a brief introduction to each of these components,

$\bullet$ Phase space volume, $h(E)$: refers to the volume enclosed by an energy hypersurface $E$ [24],

[TABLE]

where $\Phi(r)$ describes the total gravitational potential,

[TABLE]

and $\rho_{i}(r)=4\pi\int_{h(\Phi(r))}^{\infty}dh^{\prime}\frac{f_{i}(h^{\prime})}{g(h^{\prime})}\sqrt{2[E(h^{\prime})-\Phi(r)]}$ refers to the mass density for every species. Finally, $g(h)\equiv\frac{dh(E)}{dE}$ denotes the density of states [24].

$\bullet$ ** Mass Flux, $\textit{F}_{c}(h,t)$:** describes the mass flux through $h$ and is given by,

[TABLE]

where $A_{c}(h)$, $D(h)$ refer to the advection and diffusion coefficient, respectively,

[TABLE]

The summation in Eqs. (2.5) and (2.6) is over different species. While the advection term is proportional to the mass of individual species, the diffusion coefficient is exactly the same for all of them. In our analysis, we adopt a Coulomb logarithm [25] $\ln{\Lambda}\simeq\ln{\left(\frac{M_{\bullet}}{M_{*}}\right)}\simeq 10$.

$\bullet$ ** Sink term, $\nu_{c}(h,t)$:** SMBH acts like a sink as it consumes stars and compact objects which arrive sufficiently close to it. This happens in one of the two ways. Stars get disrupted by the tidal force of the SMBH at the distance $R_{tid}\equiv R_{\star}\left(M_{\bullet}/M_{\star}\right)^{1/3}$ where $R_{\star}$ and $M_{\star}$ refer to the radius and the mass of a typical star, respectively. The mass and the radius of the star are correlated, and we use the fitting formula in Ref. [26]. Not the entire mass of the disrupted star is added to the mass of SMBH. Only a fraction of its mass contributes in enhancing the SMBH’s mass. A new parameter, $f_{dis}$, is introduced to refer to the fraction of mass of a typical star that would be added to the central BH mass, $f_{dis}M_{\star}$. In the following, we consider two different examples with $f_{dis}=(1\%,10\%)$. Compact objects in the form of stellar BHs, on the other hand, are being eaten by the SMBH within the capture radius $R_{cap}=8GM_{\bullet}/c^{2}$. Here we assume that the entire mass of the swallowed BHs are added to the SMBH mass.

In the following, we use $R_{tid}$ and $R_{cap}$ for the so called loss-core radius, hereafter $r_{LC}$. Consuming the neighboring objects by the SMBH lead to an increase in its mass. Here we follow the approach of Ref. [47] to model this effect using the one dimensional orbit averaged Fokker-Planck equation. Using this approach, the loss term is estimated as,

[TABLE]

where we have removed the sub-index $c$ for the sake of simplicity. $q$ refers to the loss-cone filling factor and is given by [24],

[TABLE]

where $\mathcal{R}_{LC}(E)$ refers to the width of the loss-cone and is given by Eq. (2.8). $L^{2}_{cir}(E)$ and $L^{2}_{LC}$ denote the angular momentum of the circular orbits and the loss-cone, respectively. In addition, $P(E)$ refers to the period of nearly radial orbits with energy $E$.

Finally $\mu(E)$ denotes the orbit averaged diffusion coefficient and is given by,

[TABLE]

where $r_{max}$ refers to the maximum radius accessible to star with total energy $E$ and is given by $\Phi(r_{max})=E$. In addition, $\langle\Delta v^{2}_{\perp}\rangle$ is given by,

[TABLE]

In Sec. 4.2 we consider a signature of the above sinking process which is shown as an enhancement of the SMBH mass with time.

$\bullet$ ** Source term, $S(t)$:** We have also added a source term in the Fokker Planck equation, $S(t)$, which refers to the continuous star formation. Here for simplicity we take this function to be constant in time and we consider the source to happen interior to a given radius with a density profile proportional to $\sqrt{1/r-1/r_{source}}$. So we add two free parameters in our fittings in the next section as the mass of the source as well as the enclosed radius $r_{source}$. Work is in progress to generalize this to the case with an extended source with varies both in time and radii [67].

3 Initial density profile

3.1 Parameterization of the initial density profile

One of the key ingredients in the Fokker Planck approach is the initial condition for the density profile of different species in the system. Initial conditions are very important since the subsequent evolution of the system strongly and non-linearly depends on their selection. However, since we are dealing with very long period of evolution we have to be very careful in identifying the initial conditions. The best way is to use some parametric choices for the initial density profile and figure out the values for different parameters by comparing the “evolved” profiles with the current constraints in the galactic center. Having this said, in what follows, we first make some selections for the parametric initial conditions. Then, we determine different parameters using the direct comparisons with observations.

Hereafter, we take the following functional form for the initial density profile of different species,

[TABLE]

There are five different parameters in Eq. (3.1). Not all of the parameters above are important, though. Here we fix two of these parameters, $\alpha=4$ and $\beta=5$ and find the rest of them by a direct comparison with observations. We have checked that changing the above values for $\alpha$ and $\beta$ do not affect our results. We are therefore left with three parameters in the density profile. As we will mention, due to computational reasons, in our fitting we replace the density normalization with the mass normalization. This together with $R_{scale}$ and $\gamma$ would be determined later.

While the scaling radius, $a$, and the density slope, $\gamma$, are taken to be the same for different species, the overall normalization in the density, or equivalently in the mass, depends on the mass and abundance of them. Therefore we need to fix the mass and the number of different species in the simulation box. We choose 5 different species in the entire of our analysis. One in the form of star in the main sequence with a characteristic mass $m_{\star}=1M_{\odot}$. This is due to the fact that heavier stars have shorter lifetimes. Therefore the dynamically dominant component of stars at the present time would be the lower mass stars. We have also considered four BH species with different masses and abundances.

Using the COSMIC code, we find the BH remnant mass as a function of the zero age main sequence (ZAMS) mass. This strongly depends on the value of the initial metallicity, Z, that we take to be rather small $Z=0.001$. Using the remnant BH-ZAMS map and by assuming a Kroupa initial mass function (IMF) for the stars, we find the IMF for the BHs. Figure 1 presents the IMF of BHs compared with the initial Kroupa mass function of the stars. Using the BH IMF, we make some bins in the mass of BHs and find the mean mass as well as the normalization of them with respect to the stars. This is done in a few steps. First we compute the overall normalization in the entire of BHs using the Kroupa mass function for the stars.

[TABLE]

where $\Phi^{*}\equiv dN_{\star}/dM_{\star}$ denotes the Kroupa mass function for the stars and $m_{rem}$ refers to the remnant mass as a function of ZAMS inferred from the COSMIC code. In the numerator, we take $m_{min}=18.59M_{\odot}$ as a threshold to get BH as a remnant (using the COSMIC code). We also limit ourselves to masses less than $m_{max}=100M_{\odot}$. The lower limit on the denominator is chosen to be $m_{min}=0.05M_{\odot}$ as well.

Next, we use the BH IMF and make some regular bins in the entire range of BH mass to compute the density normalization in a given bin. Eq. (3.3) presents BH density normalization in bin (i). Here $\Phi_{BH}$ refers to the BH IMF, as shown in Figure 1, $m_{i}$ and $m_{i+1}$ denote the lower and upper limits of the BH mass in i-th bin and $m_{min,BH}=2.86M_{\odot}$ and $m_{max,BH}=40.35M_{\odot}$ refer to the minimum and maximum values of the BH mass.

[TABLE]

We also compute the mean value of the BH mass in every mass bins. This gives us a consistent link between the BHs mass/normalization and the ZAMS normalization. Table 1 presents the mean mass as well as the normalization of the BH densities for our four BH species. We emphasis here that the second column gives us the “fractional” normalization in the BH mass as refereed by F[i]. Therefore the actual normalization in the density/Mass of the BHs is given by this number multiplied with the actual normalization of the ZAMS, $M_{\star}$, which is a parameter in our fits. Since the BHs are the remnant of the stars in the ZAMS, we take the normalization of the stars to be $(F[0]=1-\sum_{i}F[i])\times M_{\star}$.

In summary, we are left with three free parameters in the density profiles. On the other hand, due to the consumption of the disrupted stars and swallowed BHs, central BH mass increases with the time. Therefore, its initial mass is another free parameter that must be determined as well. Finally, there are two more parameters associated with the source term, i.e. fractional source mass and its radius, that should be found as well. Here we split the total mass in every species to two different parts. One is the continuously formed element, $fraC_{source}F[\mu]M_{\star}$, where $fraC_{source}$ refers to the fraction of the continuously born source and $\mu=(0,1,2,3,4)$ refers to both of the stars and different BH species. And, the second part is the initial normalization $\left(1-fraC_{source}\right)F[\mu]M_{\star}$. We name the source radius as $r_{source}$. This gives us a family of 6 parameters that should be fixed using the observations.

3.2 Inferring the initial density profile using the observations

Having introduced a model with 6 free parameters, here we attempt to find each of these parameters by using a direct comparison with observations. We particularly use the data from [49] including the flux density, see the left panel in Figure 9 of [49], at different radii as well as the enclosed mass at two different radii, namely $M(r\leq 1pc)=10^{6}M_{\odot}$ and $M(r\leq 4pc)=10^{7}M_{\odot}$. Finally, we assume that after 10.5 Gyr the central BH mass is $M_{\bullet}=4\times 10^{6}M_{\odot}$. Using the above observations and by performing an MCMC analysis, we find the proper values for the above 6 parameter family. The only left over parameter is the fraction of the tidally disrupted stars added to the central BH mass. Since only a small fraction of the tidally disrupted stars are added to the central BH mass, here consider $f_{dis}=(1\%,10\%)$ and find the actual fit for both of these choices. We aim to find some ranges for the density slope, $\gamma$, logarithm of the scaling radius, $\log_{10}(R_{\mathrm{scale}})$, logarithm of the initial total mass normalization, $\log_{10}(M_{\mathrm{*,init}})$, logarithm of the initial central BH mass, $\log_{10}(M_{\mathrm{\bullet,init}})$, fraction of the source mass, $\mathrm{frac}_{\mathrm{source}}$, and the logarithm of the source radius, $\log_{10}(r_{\mathrm{source}})$.

3.2.1 Parameter fits for $f_{dis}=1\%$

Here we find the actual fit to observations when only one percent of the disrupted stars are added to the central BH mass. Figure 2 presents the posterior plot for the range of the parameters. One sigma range of the above 6 dimensional family are given as,

[TABLE]

Having determined the median and 16, 84 percentile for the above parameters, in the following, we use the median values and we run the Phase-Flow code to find the actual evolution of the system.

3.2.2 Parameter fits for $f_{dis}=10\%$

Here we find the fits for the case with 10% of the tidally disrupted stars added to the mass of the central BH mass. Figure 3 presents the posterior plot for the range of the parameters. One sigma range of the above 6 dimensional family are given as,

[TABLE]

Comparing the case with $f_{dis}=1\%$ to $f_{dis}=10\%$, we clearly see that the fitted parameters are very close. Therefore, in what follows we only consider the case with $f_{dis}=10\%$.

4 Two body relaxation using the one dimensional Fokker-Planck approach

So far we have fixed the initial conditions of simulation using the actual fits to the observations. Hereafter, we take the median values for the above 6 parameter family and simulate the system for an extended period of time, t = 10.5 Gyrs. We aim to study the dynamics of system thoroughly. For this purpose, we consider several quantities including the relaxation timescale, mass growth of the central BH, consumption of the stars and BH to the central BH and the evolution of the density profiles in details. As we will see due to the complexity of having different mass components in the same simulation the evolution is highly non-trivial.

4.1 Relaxation Timescale

We start by considering the behavior of the relaxation timescale. For a single mass component, it is defined as,

[TABLE]

where $\sigma$ refers to the velocity dispersion while $D[\left(\Delta V_{||}\right)^{2}]$ is given by,

[TABLE]

In our case, however, we are dealing with multiple masses in the same box and so we should generalize Eq. (4.1) to case with multiple mass species. There is not a unique way of doing this. Here propose an effective way to compute the relaxation timescale. We compute an averaged velocity dispersion in the numerator of Eq. (4.1) while we add different contributions in Eq. (4.2). We should however notice that this is only an effective description and in reality heavier BHs may have shorter relaxation timescales. Figure 4 presents the radial behavior of the effective relaxation timescale for few different times. While the typical relaxation timescale is rather long in the galactic nuclei, it decreases by few orders of magnitude toward smaller radii.

4.2 Evolution of central BH mass

Here we study the growth of the mass of central BH with time. Figure 5 presents the dynamical evolution of the central BH mass due to the tidal disruption of stars and BHs with time. Left panel shows a comparison between the full numeric (blue line) with the contribution from different species (both of stars and BHs with different masses) (dotted dashed Orange line). The results are completely matched. According to the plot, the fractional changes in the central BH mass is about 20% during the entire evolution of the system. Right panel shows different components individually. Our results are highly non-trivial. Although we assume a 10% contribution from the stars, they (almost) dominate the mass growth in the central BH. Although different BHs have the same overall shape, the contribution from the BH with mass $m_{BH}=23.92M_{\odot}$, the second heaviest BH, is dominated over the rest of the BHs. A quick look at Tab. 1 shows that it is because the initial normalization of this BH, inferred from the BH IMF, is relatively higher than the others. This can be interpreted as the direct impact of the IMF of BH in the evolution of system. We emphasize that this is the first time that the detailed impact of IMF of BH is shown in the tidal disruption event.

4.3 Averaged mass with radius

Next, we consider the radial profile of average masses around the central BH,

[TABLE]

Figure 6 presents the average mass computed using Eq. (4.3) at different times. From the plot it is clear that averaged mass is quickly redistributed with time early on while shows slower changes later on. Furthermore, while it increases toward the center, it get suppressed at larger radii.

4.4 3D/2D density profile of stars and comparison with simulations/observations

Another interesting quantities are the evolution of the density profile of stars in 2D and 3D. The aim is to compare our results with the recent observational results [49], in the case of 2D, as well as N-body simulations [56, 57, 58]. Figure 7 presents the 3D density profile (left) and 2D surface density (right) with time. In 3D density profile, dashed red line is a fit to Bahcall-Wolf (BW) cusp [1, 55] with the slope $\gamma=-5/4$. In 2D surface density profile, we compare the final surface density with the recent observational results in [49]. Our median parameters inferred from MCMC analysis, gives us an incredible fit to the most recent observations in the galactic center.

4.5 Dynamical evolution of density profile of BHs

Next we consider the dynamical evolution of BHs density profile with time. Likewise the case of stellar profile, we aim to compare our results with those of N-Body simulations. The main advantage of our approach is that using the Fokker Planck, we may extend the density profiles toward smaller radii which are completely inaccessible with the costly N-Body simulations. Figure 8 presents BHs density profiles with time for different BH masses. The growth of the BW cusp is clearly seen from the plot. Since different BH masses are connected through their IMF, we may want to check how does this affect their final density profiles. In figure 9, we present density profile of BHs after $t=10.5$ Gyrs. From the plot, we clearly see the impact of IMF of BHs in enhancing mostly both of the density and mass of second heaviest BH, $m=23.92M_{\odot}$ at every radius. Heaviest BH, with $m_{BH}=34.5M_{\odot}$ experiences faster movement to the center which leads to further suppression of its density profile at larger radii.

4.6 Velocity dispersion of BHs and impact of the central BH

Finally we study the behavior of velocity dispersion defined as,

[TABLE]

here $\rho_{i}(r)$ refers to the density profile of stars as well as BHs and $M(r)$ denotes the total mass interior to the radius $r$. As expected there are different contributions in the interior mass including the stars, different BH species and central BH. Figure 10 presents the radial profile of the velocity dispersion. On the left panel, we present the stellar contribution in the velocity dispersion. From the plot it is clear that at smaller radii, $r/pc<0.1$, the impact of the central BH is dominated in the velocity dispersion. However at larger radii the impact of the stars get dominated thanks to their shallower cusp profile. Right panel compares the velocity dispersion of stars and BHs. It is clearly seen that their velocity profile are quite similar to each other except a little difference for $0.1<r/pc<10$ due to the transition from being dominated by the central BH or stars.

5 GW Implications

So far we presented various aspects of having BH mass-function in the galactic center. Here we investigate another observational signature of such a mass function in the connection with the GW. We start with the detailed description of the steps involved. Then we present the forecast for the expected number of the events after a maximum operation of LISA of about $T_{obs}=10$yrs.

5.1 Set up

Here we describe the setting of our simulation for counting the number of expected events as seen with LISA. Our method is based on a post-processing of the AGAMA code where we add in the impact of GW as well as the angular momentum diffusion in the evolution of the semi-major axes, $a$, and eccentricity, $e$. We make a sample of $N_{tot}=10000$ systems made of pairs of initial $(a,e)$, 100 values in each of them. Initial values of semi-major axes are taken from a distribution of $f(a)/a$, where $f(a)$ refers to the phase-space distribution function. We allow the system to get relaxed for an amount of $t=3.7$ Gyrs where we see that it does not change dramatically afterward. We use this distribution function in our samples and find a grid of semi-major axes in the range of $(7\times 10^{-6}-10)$ pc. In addition, initial values of the eccentricity are chosen from a thermal distribution. Taking these initial conditions, we evolve the system for a maximum period of $t=6.6$Gyrs or when any given systems cross the loss-cone surface, as is determined in the following.

5.2 Dynamical Evolution

As already described above, in our simulations we make a sample of $N_{tot}=10000$ systems with the above initial conditions. Next, we consider the dynamical evolution of every systems where the dynamics are given as,

[TABLE]

where $M_{tot}\equiv m_{bh}+M_{\bullet}$. Here the second term in Eq. (5.2) refers to the angular momentum diffusion term , as is discussed in more details at [67]. We add it for taking into account the impact of the angular momentum diffusion which are averaged in the 1D Fokker Planck approach. Though the details of this term is left to [67], it is arisen from the loss-rate as described in Eqs. (2.7) and Eq. (2.8) with converting the circular orbit, $e=0$, to the elliptical orbit with non-zero eccentricity. It is fair to say that this novel term was ignored in the analysis before and as it turns out this is very important contribution in enhancing the eccentricity and so pushing the system to be potentially detectable with LISA. Though not having the same origin, it is similar to the enhancement in the rate of the GW due to the famous Kozai-Lidov oscillations in the context of the triple systems. Owing to this term, a lot of orbits experience a horizontal evolution, along the eccentricity axes, before they get to the LISA band. This term is responsible in increasing the number of the expected signals in LISA which was almost entirely ignored by the previous authors [68].

5.3 Signal to Noise Ratio

Having presented the simulation set-up as well as the dynamical evolution for the system, here we present our formalism for computing the signal to noise ratio in the LISA band [69]

[TABLE]

where $n_{max}$ refers to the maximum number of the harmonics. In our analysis we take this number to be maximum $10^{5}$. In addition, $h_{c,n}=\left(\pi d\right)^{-1}\sqrt{2\dot{E}_{n}/\dot{f}_{n}}$. $d$ also refers to the distance from the source. Finally, $S(f_{n})$ refers to the LISA noise power. Since we stop the integration before the loss-cone, we can directly check that in all of our samples we have $\dot{f}/f\ll T_{obs}\simeq 10$yrs. Therefore we can expand the integral in Eq. (5.3) as,

[TABLE]

where $M_{c}$ refers to the chirp mass and we set $d=7.8kpc$. Furthermore, in our estimation, we take the maximum period of LISA operation of about $T_{obs}=10$ yrs and we choose the criteria of having $SNR=8$ in inferring the detectable number of systems.

5.4 Estimation of rate of events

Here we combine the above individual pieces and describe our method in estimating the total rate of expected LISA events after $T=10$ yrs of LISA operations, where we have taken the maximum duration of LISA to increase the chance of detecting LISA sources.

Starting from the initial condition for every pairs of $(a,e)$ as mentioned above, we evolve every systems with time either to the time that they cross the loss-cone area, defined as $a(1-e)=8GM_{\bullet}$, or up to $t=6.6$ Gyrs. Using Eq. (5.4), we compute the signal to noise ratio for each of the above systems during the entire of their lifetime, defined from $t=0$ to either their loss-cone crossing time or up to $t=6.6$ Gyrs depending on which happens earlier. For every systems, we compute the total duration of time that the system has $SNR\geq 8$ with a frequency in the LISA band. We call this the lifetime of being in the LISA and divide it to the total lifetime of this system. For example, if the system has a total lifetime of $t=10^{-6}$ Gyrs and spends a period of $\Delta t=10^{-7}$ Gyrs in the LISA band, the weighting factor for this system would be $w=10^{-7}/10^{-6}$. Next, we should take into account the rate of the replenishment for every systems. Since the replenishment is due to the energy diffusion, we use the diffusion rate as an extra factor for every systems. Therefore we get a combined rate factor, hereafter $w(i)$, for ith system as,

[TABLE]

where $\Delta t_{i}$, $t_{i}$ and $\mu_{i}$ refer to the spent period of time in the LISA band, the total lifetime and the diffusion rate for i-th system, respectively.

At the next level, we should rescale the above number to the actual number of the BHs in the galactic center interior to the radius $r=r_{sam}=10pc$, the upper limit in our sample making, as describe in 5.1. This can be done by the factor $N_{BH}(r\leq r_{sam})/N_{tot}$ where $N_{BH}(r\leq r_{sam})\equiv M_{BH}(r\leq r_{sam})/m_{BH}$ refers to the number of the BHs with mass $m_{BH}$ interior to the radius $r=r_{sam}=10pc$ and $N_{tot}=10000$ is the total number of the samples in our system. Therefore, the final rescaled rate for i-th system would be, hereafter $W_{tot}(i)$,

[TABLE]

Furthermore, we should also eliminate the cases with semi-major axes above a critical value, hereafter $a_{crit}$, where the timescale of the GW, $\tau_{GW}\equiv|a/\dot{a}|$, is longer than the loss-cone timescale $t_{LC}\equiv\left(L_{LC}/L_{cir}\right)^{2}t_{Rel}$. Where the loss-cone timescale is defined as a timescale associated with a change in the angular momentum by the order of the loss-cone [4]. Combining with the loss-cone surface, $a(1-e)=8GM_{\bullet}$, we can estimate the critical semi-major axes above which $\tau_{GW}>t_{LC}$. As argued in [4], systems with $a>a_{crit}$ scatter either inward the loss-cone or to larger orbits without emitting the GWs. Therefore they are not source of the inspiral phase and must be removed from our sample. In Table 2 we estimate $a_{crit}$ for different BH masses.

In order to get the rate of the observable systems in our sample, $R_{obs}$, we sum over the total rescaled rate of any systems as,

[TABLE]

The total number of the LISA inspiral events after $T_{obs}=10$ yrs would be,

[TABLE]

5.5 Results and Discussions

Having described our novel formalism in computing the total expected numbers of the LISA events, here we present the final number for our system of 4 BHs. Table 3 presents the expected inspiral rate with LISA as well as the total number of LISA sources from the inspiralling BHs from the galactic center after $T_{obs}=10$ yrs. There are few important points that are worth highlighting here,

First, the rate for one galaxy is rather small. Therefore, it is of interest to use this novel method for different galaxies, taking into account their internal structures, and read off the final expected rate and number. If we focus on Milky Way like galaxies, with similar black holes to SgrA*, we could imagine that combining $N_{gal}\sim 10^{4}$ of them we may be able to see one inspiralling signal. Care should be taken for other type of galaxies with different structures, though. In [67], we generalize this formalism to more galaxies as well as larger samples of the BHs in hand and also different metallicities.

Second, the impact of the BH mass-function is clearly seen in the total expected number of events in LISA. Indeed, comparing Table 1 with Table 3, we clearly see that the total rate is peaked at the peak of the mass-function. This gives us a another signature of the BH mass-function this time in the GW signal.

6 Conclusions

Using an orbit averaged tool, AGAMA code, we simulated the galactic center in the presence of stars in the main sequence and 4 different BH species. We used a binary stellar evolution code, COSMIC, to create a sample of BHs from an initial low metallicity of $Z=0.001$ and constructed the BH mass-function initially. For the first time in studies of galactic nuclei, we explored the impact of the BH mass-function in the evolution of the system. We have also added the impact of the continuous star formations in the system. We used an MCMC approach for computing the initial profile of stars and BHs which are matched with the current observations after a total duration of $t=10$ Gyrs. We implemented the above results in the Phase-Flow code, a library in the AGAMA code, and studied the dynamical evolution of this system precisely. The novelty of this work involves the following results.

$\bullet$ First , the above hybrid tool of the stellar evolution with the AGAMA code enabled us to figure out the impact of the BH mass-function in the evolution of the system through its direct effect on the initial normalization of different BH species. As shown in Table 1, the normalization of the second heaviest BH is dominated over the rest of the BHs thanks to a local peak in the BH mass-function.

$\bullet$ Second, we explicitly showed that it leads to an enhancement in the contribution of the second heaviest BH in the growth of the central BH’s mass as they are more abundant in the system. Quite interestingly, as shown in the right panel of Figure 5, taking $f_{dis}=10\%$, stars make the second important contribution in the mass of the central BH.

$\bullet$ Third, detailed study of the velocity dispersion profile revealed that very close to the Central BH, stars and BHs behave similarly in the velocity dispersion. This would change in some intermediate scales where the contribution of stars begin to dominate over the contribution of the central BH which leads to a deviation from the Keplerian orbits thanks to shallower cuspy profile of stars compared with the BHs. Observationally it is very interesting and may enable us to measure the cusp profile by looking at the projected velocity dispersion close in the center. We leave further investigation of this to a future work. getting important. Finally, at larger radii stars are the dominant source of the velocity dispersion.

$\bullet$ Finally, we considered the GW in this system. Since the orbit averaged Phase Flow does not consider the angular momentum diffusion nore the GW dissipation, we made a post-processing analysis where we took into account both of these effects self consistently. For the first time in this context, we considered these effects collectively and explicitly showed an enhancement in the eccentricity from the angular momentum diffusion. We made a sample of 10000 pairs in the initial semi-major axes and eccentricity and evolved each of them for a maximum period of $t=6.6$ Gyrs. We computed the SNR for each of them during the entire of their lifetime and computed the fraction of the time that they spend in the LISA band, defined with the criteria of having $SNR\geq 8$ and in the LISA frequency band to their entire lifetime. WE weighted this number with the diffusion rate as the replenishment factor. We also rescaled that number with the actual number of the BHs interior to $r=10pc$ as the upper limit in our sample making. We presented the total rate at Eq. (5.7). Furthermore, we also eliminated the systems with $a>a_{crit}$ where $a_{crit}$ refers to the maximum semi-major axes that above it the characteristic timescale of the GW is above the loss-cone time-scale. BHs in this area either scatter off the central BH to larger orbits or get swallowed to the center without emitting any GWs. We removed these systems from our samples entirely. Our computation showed that after a maximum duration of $T_{obs}=10$ yrs, we get a rate of $10^{-6}-10^{-5}$ 1/yrs as the inspiralling rate for the BHs with different masses. The total expected rate and number of one MW like galaxy is given in Table 3. This means that in order to get one signal from the inspiral phase we need to consider of order of $N_{gal}\simeq 10^{4}$ MW like galaxies. Quite interestingly the inspiral rate is peaked at the peak of the BH mass-function. This is the second signature of the the BH mass-function proposed here.

Acknowledgment

We thank Christopher Berry, Sownak Bose, John Forbes, Dan D’Orazio and Morgan MacLeod for helpful discussions. We especially thanks Eugene Vasiliev and Fabio Antonini for very insightful discussions and for a followup collaboration. We are grateful to the anonymous referee for their very constructive comments. R.E. acknowledges the support by the Institute for Theory and Computation at Center for Astrophysics, Harvard-Smithsonian. This work was also supported in part by the Black Hole Initiative at Harvard University which is funded by a JTF grant.