Methods

Merger tree

Where we got the merger data from and how we extracted it. - Using Andrea Kulier’s cosmological simulation, we plotted galaxy masses as a function of time. - We took the data from the three most massive galaxies (ie, most massive at redshift z=0). - We curve-fit each galaxy’s mass to polynomial functions of time. - We also curve-fit each galaxy’s (accreted+seeded) central black hole mass to exponential functions of time. - We added those functions to the Fortran hermite code - We implemented the SO potential-density pair (using the m200, r200 from the fitting functions above, getting half-mass radius from NFW profile, and using it as half-mass radius for the SO profile).

Simulation of Cen et al. (have to look up the year)

AR-Chain code

Summary of the code and the modifications we made.

For the numerical simulations presented here, we used a modified version of the algorithmic chain integrator AR-Chain developed by \citet{Mikkola06}. It uses algorithmic chain regularization for high-precision integration of few-body dynamics, and is capable of handling velocity-dependent forces efficiently. It includes relativistic post-Newtonian terms up to order PN2.5 \citep{Mikkola08}.

Galaxy background potential

Phase-space diffusion

Weak encounters with background stars will let the SMBHs diffuse through phase space while they are orbiting within the gravitational potential of the galaxy. The diffusion can be expressed as change in velocity of an SMBH by \(\Delta \vec{v}\) per unit time. We can split this change into a component along the direction of motion of the SMBH, and one perpendicular to that. Following \citet{Binney08}, the diffusion coefficients can be expressed as \[\begin{aligned} D[\Delta v_\parallel] & = & -\frac{4\pi G^2\rho(r)M_\bullet\ln\Lambda}{\sigma^2}f(\chi),\label{eq:df}\\ D[(\Delta v_\parallel)^2] & = & \frac{4\sqrt{2}\pi G^2\rho(r)M_\bullet\ln\Lambda}{\sigma}\frac{f(\chi)}{\chi},\\ D[(\Delta \vec{v}_\bot)^2] & = & \frac{4\sqrt{2}\pi G^2\rho(r)M_\bullet\ln\Lambda}{\sigma}\left[\frac{\mbox{erf}(\chi)-f(\chi)}{\chi}\right],\end{aligned}\] where \(\Delta v_\parallel \equiv \Delta \vec{v}\cdot\vec{v}/v\) is the velocity change in direction of motion, and \(\Delta \vec{v}_\bot \equiv \Delta \vec{v} - \Delta v_\parallel \cdot\vec{v}/v\) is the velocity change perpendicular to the direction of motion. Here, \(M_\bullet\) is the mass of the black hole, and \(\chi = \frac{v}{\sqrt{2}\sigma(r)}\). The function \(f(\chi)\) is given by \[f(\chi) = \frac{1}{2\chi^2}\left(\mbox{erf}(\chi)-\frac{2\chi}{\sqrt{\pi}}\exp\left(-\chi^2\right)\right).\] We approximate the factor \(\Lambda\) in the Coulomb logarithm as \[\Lambda = \left(\frac{M_{NSC}}{M_\bullet}\right)\left(\frac{r}{r_h}\right).\] We can identify Eq. \ref{eq:df} as the dynamical friction term, that is, if we assumed \(D[(\Delta v_\parallel)^2] = D[(\Delta \vec{v}_\bot)^2] = 0\), we would get Chandrasekhar’s dynamical friction formula. The second term introduces a variance of the friction term, and even allows the SBHs to be accelerated when the velocity of a SBH gets sufficiently small. The third term introduces a change in velocity perpendicular to the direction of motion of the SBH. It is a randomly oriented vector, and hence causes the SBHs to execute a random walk in phase space. The last two terms will establish that the SBHs are ultimately in energy equipartition with the background stars. The velocity changes \(\Delta v_\parallel\) and \(\Delta\vec{v}_\bot\) per unit time \(\Delta t\) can be computed with the above equations. Both changes are normal distributed, where the mean, \(\mu\), and the variance, \(\Sigma\), of the distributions are given by \[\begin{aligned} \mu_\parallel &=& D[\Delta v_\parallel]\Delta t,\\ \Sigma_\parallel &=& D[(\Delta v_\parallel)^2]\Delta t,\\ \mu_\bot &=& 0,\\ \Sigma_\bot &=& D[(\Delta \vec{v}_\bot)^2]\Delta t.\end{aligned}\] We compute the diffusion coefficients for each black hole at each time step, and modify its velocity on a Monte Carlo basis. For each time step we draw a random orientation before adding the perpendicular velocity change to the respective SBH. Hence, the SBH’s modified velocity, \(v_f\), is computed using \[\begin{aligned} \vec{v}_f &=& \vec{v}_0 + \Delta v_\parallel \hat{v}_\parallel + \Delta v_\bot \hat{v}_\bot,\\ \Delta v_\parallel &=& \mathcal{N}(\mu_\parallel, \Sigma_\parallel),\\ \Delta v_\bot &=& \mathcal{N}(\mu_\bot, \Sigma_\bot).\end{aligned}\] The change of energy, \(\mbox{d}E_{BH}\), of the orbiting black hole due to phase-space diffusion is given back to the stellar background potential, with \(\mbox{d}E = -\mbox{d}E_{BH}\). As a consequence of this energy transfer, inspiralling black holes will cause an expansion of the NSC. For this purpose we calculate the change in potential energy, \(\mbox{d}W\), of the stellar system using \[\begin{aligned} E &=& T + W = \frac{1}{2}W,\\ \mbox{d}W &=& -2\,\mbox{d}E_{BH},\end{aligned}\] where we made use of the virial theorem \(2T+W =0\). With this change in potential energy we can calculate a new radius for the stellar background potential at each integration step. For the Plummer sphere the new scale radius can be calculated as \[a_{new} = a\left(1+\frac{32a\,\mbox{d}W}{3\pi GM_{NSC}^2}\right)^{-1}.\]

Gravitational wave recoils

The code AR-Chain includes PN terms up to order 2.5. The SMBHs can therefore merge via gravitational wave emission. We include gravitational wave recoils following the prescription outlined in \citet{Kulier_2015}, which is based on the fitting formula by \citet{Lousto12}. To save computational time, we assume that a merger will be inevitable when the separation between two SBHs gets smaller than 10000 Schwarzschild radii. At the moment of the merger, we assume that the spin vectors of the two SBHs are randomly aligned. This results in kick velocities of up to several thousand kms\(^{-1}\), with a median kick of \(\approx 290\)kms\(^{-1}\). Since our simulations focus on NSC with relatively low escape velocities, this implies that a majority of the merging SBHs escape from the NSCs.

Black holes can also eject each other via strong three-body interactions. We remove SBHs from the simulations once they move beyond 1kpc from the NSC, assuming that it will take them more than a Hubble time to find their way back into the center of the host galaxy.

Simulation setup

Injection, merging, escape, effective radius, set of galaxies