BINARY BLACK HOLE MERGER IN GALACTIC NUCLEI: POST-NEWTONIAN SIMULATIONS

The initial conditions of the N-body models used in this work are chosen to be similar to the ones used by Berczik et al. (2005, 2006), allowing for a direct comparison to their simulations and for an accordant interpretation of our results. Here, we briefly describe the main aspects of the initial model setup and refer the reader to the latter two publications (and references therein) for further details.

The initial field particle distribution—representing the galactic nucleus—is set up following the phase-space density distribution of a King model (King 1966) with internal rotation (e.g., Longaretti & Lagoute 1996; Ernst et al. 2007, and references therein). The models have been set up with a central concentration parameter of W₀ = 6 and an initial rotation of $\omega _0\equiv \sqrt{9/(4\pi \,G\,\rho _{\mathrm c})}\,\Omega _0=1.8$, where G and ρ_c are the gravitational constant and the central concentration, respectively. The net rotation in our models is chosen to be counterclockwise, with the angular momentum vector being aligned with the z-axis. We adopt the standard N-body units (see, e.g., Aarseth 2003b, and references therein) to our numerical models by setting both the gravitational constant G and the total mass M_⋆ of the N-body system to unity, and by setting its total energy to E_⋆ = −1/4.

In this work, the galactic nucleus is represented with total numbers of N_⋆ = 25 × 10³ and 50 × 10³ field particles, respectively, each with nine different random realizations. We restrict ourselves to these relatively small particle numbers as a trade-off for the large set of simulations that we have performed in grand total. However, as Berczik et al. (2006) have shown, the hardening rate of SMBH binaries in such rotating galaxy models is in this regime essentially independent of the field particle number N_⋆. The main difference between our models and the latter ones is the $\cal {P\!N\!}$ treatment of the two SMBHs. Therefore, the results of our simulations are expected to be quasi-N-independent as well, since we use a rotation parameter of ω₀ > 1.2 (Berczik et al. 2006). The field particles have all equal mass m_⋆ = M_⋆/N_⋆. The gravitational softening length for the field particles has been set to _⋆ = 10⁻⁴ (in model units). The same softening is also used for the star–BH interaction. The two BH particles themselves are evolved without any gravitational softening, i.e., with _BH = 0.0.

We set the masses m₁ and m₂ of the two BH particles to be one per cent of the nuclei mass M_⋆, i.e., in our case m₁ = m₂ = 0.01. The two BH particles are initially placed on the y-axis at y_1,2 = ±0.3, respectively, within the z = 0 plane. The BHs are given an initial velocity v_x, corresponding to roughly 10% of the local circular velocity, as derived from the enclosed mass of the underlying field particle distribution. With this we get roughly v_1,2 ≈ ±0.07, in our model units. Note that in this configuration the two BHs are initially unbound with respect to each other.

To scale our N-body results to real galaxies, we consider the total energy of the system

$$\begin{equation} E = - \alpha \, \frac{G\,M^2}{R}, \end{equation} \tag{ 3 }$$

where G, M are the gravitational constant and the total mass of our model system (typically some fraction of the bulge and cusp components in a galactic nucleus), and R, α are model-dependent quantities, giving a scaling radius and a numerical constant of order unity, respectively. As an example, for a nonrotating King model with W₀ = 6 (which is used in our models as an initial configuration), we would have R = r_c as the core radius as defined by King (1966) and α = 0.0759.

With G = 1, we can freely choose two of the three scales (energy, radius, or mass). In standard N-body units the unit of mass is the total mass, i.e., M = 1, and the unit of energy is 4 E (E = −0.25). Therefore, the radial unit is determined by the condition R = 4 α (this value determines the value of R in N-body units, e.g., for a Plummer model, we have α = 3π/64 and R = 3π/16). The physical units of mass, length, energy, velocity, and time are then given as

$$\begin{eqnarray} \left[ M \right] & = & M,\nonumber \\ \left[ L \right] & = & \frac{R}{4\alpha },\nonumber \\ \left[ E \right] & = & 4\alpha \, \frac{G\,M^2}{R}, \nonumber \\ \left[ V \right] & = & 2\sqrt{\alpha } \left(\frac{G\,M}{R}\right)^{1/2},\nonumber \\ \left[ T \right] & = & \frac{1}{8\,\alpha ^{3/2}} \left(\frac{R^3}{G\,M}\right)^{1/2}. \end{eqnarray} \tag{ 4 }$$

The speed of light c in N-body units is then

$$\begin{eqnarray} c & = & c_0 / \left[ V \right] = c_0 \left(4\alpha \frac{G\,M}{R} \right)^{-1/2} \nonumber \\ & = & 457 \left({\frac{M}{10^{11} {M}_{\odot }}} \right)^{-1/2} \left({\frac{R}{10^3\, {\rm pc}}} \right)^{1/2}, \end{eqnarray} \tag{ 5 }$$

where c₀ is the speed of light in physical units, and for the second expression we have inserted α = −0.25. This illustrates that in the $\cal {P\!N\!}$ regime the N-body problem is not scale-free anymore—to fix the scale for c means fixing the radial scale or vice versa. In our simulations, we have chosen different values for the speed of light, i.e., c = 447, 141, 44, and 14, in order to enhance the effect of the relativistic corrections. At the same time, this variation of c scales the Schwarzschild radius R_BH = 2 G M_BH/c² of the two BHs to values of 10⁻⁷, 10⁻⁶, 10⁻⁵, and 10⁻⁴, respectively, in model units.

SMBH binaries are expected to be subject to the emission of GWs. Since the GWs drain energy and angular momentum from the binary, its orbital elements will change in the course of their dynamical evolution. Already more than four decades ago, Peters & Mathews (1963) and Peters (1964) derived the change of energy E and of the orbital elements for a Keplerian orbit, namely its semimajor axis a and eccentricity e, under the influence of GW (quadrupole) emission. The orbit-averaged expressions for two compact objects with masses m₁ and m₂ orbiting each other are given as

$$\begin{eqnarray} & & \left\langle \!\frac{da}{dt}\!\right\rangle = -\frac{64}{5} \frac{G^3\, m_1\, m_2 (m_1+m_2)}{c^5\, a^3} f(e), \\ & & \left\langle \!\frac{de}{dt}\!\right\rangle = -\frac{304}{15} \frac{G^3\, m_1\, m_2 (m_1+m_2)}{c^5\,a^4 (1-e^2)^{5/2}} \left(e+\frac{121}{304}e^3\right), \\ & & \left\langle \!\frac{dE}{dt}\!\right\rangle = -\frac{32}{5} \frac{G^4\, m_1^2\, m_2^2 (m_1+m_2)}{c^5\,a^5} f(e), \end{eqnarray} \tag{ 6 }$$

where the so-called enhancement factor f(e), given by the following function, shows a strong dependence on the orbital eccentricity e:

$$\begin{equation} f(e) = \left(1-e^{2}\right)^{-7/2}\left(1+\frac{73}{24}e^{2}+\frac{37}{96}e^{4}\right). \end{equation} \tag{ 9 }$$

We note that the rate at which energy E is lost due to the emission of GWs depends strongly on the eccentricity e of the orbit. Hence, compact binaries with high eccentricities (i.e., small a and e ≈ 1) are expected to be strong sources of GWs. The equations of Peters & Mathews (1963) describe the shrinking of the semimajor axis (corresponding to an inspiral of the two objects) and the circularization of the orbit. According to Peters & Mathews (1963) and Peters (1964), the typical timescale of coalescence due to the emission of gravitational radiation is given by

$$\begin{equation} t_{\mathrm{gr}} = \frac{5}{64} \frac{c^{5}a_{\mathrm{gr}}^{4}}{G^{3}m_{1}\,m_{2}(m_{1}+m_{2})\,f(e)}, \end{equation} \tag{ 10 }$$

where a_gr denotes the characteristic separation for GW emission.

To account for the relativistic effects in our numerical N-body simulations, we use the $\cal {P\!N\!}$ formalism. The equations of motion for a compact binary system are written in harmonic coordinates (Schutz 1985) and are defined in the inertial N-body frame of reference. As they are relativistic, they are (1) invariant under $\cal {P\!N\!}$-expanded Lorentz transformations, (2) reduce to the geodesics of the $\cal {P\!N\!}$-expanded Schwarzschild metric in the limit in which one of the masses goes to zero and neglecting the spin, and (3) are conservative if the $2.5 \,\cal {P\!N\!}$ radiation reaction term is turned off. In fact, an isolated $2 \,\cal {P\!N\!}$ binary should conserve its generalized $\cal {P\!N\!}$ integrals of motion up to ${\cal O}(1/{\mathrm c}^6)$ (Andrade et al. 2001). We implement the $\cal {P\!N\!}$ equations of motion formulated in the absolute Euclidean space and time of Newton (e.g., Blanchet 2006, Equation 168 therein). In the present work, we apply all $\cal {P\!N\!}$ corrections up to the order ${\cal O}(1/{\mathrm c}^5)$, i.e., the $2.5 \,\cal {P\!N\!}$ corrections is the highest order that we take into account. While the $2.5 \,\cal {P\!N\!}$ correction accounts for the emission of GWs, the $1 \,\cal {P\!N\!}$ and $2 \,\cal {P\!N\!}$ terms conserve energy. However, the latter two terms result in precession of the orbital pericenter. Similar to the equations of motion in the center of mass frame (Kupi et al. 2006), one can write the accelerations, e.g., for particle 1 of the binary, in the following form:

$$\begin{equation} {\bf a}_1 = - \frac{G\,m_2}{r_{12}^2} \left[ (1+{\cal A})\, {\bf n}_{12} + {\cal B}\,{\bf v}_{12} \right]\,, \end{equation} \tag{ 11 }$$

where r₁₂ is the separation of the two particles, n₁₂ is the normalized relative position vector, and v₁₂ is the relative velocity. The two functions $\cal {A}$ and $\cal {B}$ contain the different orders of the $\cal {P\!N\!}$ expansion and can be written as

$$\begin{eqnarray} {\cal A} & = & \frac{1}{c^2}\, {\cal {A}}_{1{\cal P\!N}} + \frac{1}{c^4}\, {\cal {A}}_{2{\cal P\!N}} + \frac{1}{c^5}\, {\cal {A}}_{2.5{\cal P\!N}} + {\cal O}\left(\frac{1}{c^6} \right), \\ {\cal B} & = & \frac{1}{c^2}\, {\cal {B}}_{1{\cal P\!N}} + \frac{1}{c^4}\, {\cal {B}}_{2{\cal P\!N}} + \frac{1}{c^5}\, {\cal {B}}_{2.5{\cal P\!N}} + {\cal O}\left(\frac{1}{c^6} \right). \end{eqnarray} \tag{ 12 }$$

In this notation, the first $\cal {P\!N\!}$ correction ($1 \,\cal {P\!N\!}$) is, for example, given as

$$\begin{eqnarray} {\cal A}_{1{\cal P\!N}} & = & \left[ \frac{5 G\,m_1}{r_{12}} + \frac{4 G\,m_2}{r_{12}} + \frac{3}{2}({\bf n}_{12}{\bf \cdot }{\bf v}_2)^2 - {\bf v}_1^2 \right. \nonumber \\ & & +\left. 4({\bf v}_1{\bf \cdot } {\bf v}_2) - 2{\bf v}_2^2 \vphantom{\frac{5 G\,m_1}{r_{12}}}\right] \,, \\ {\cal B}_{1{\cal P\!N}} & = & 4\,({\bf n}_{12}{\bf \cdot } {\bf v}_1) - 3\,({\bf n}_{12}{\bf \cdot }{\bf v}_2) \,. \end{eqnarray} \tag{ 14 }$$

We forgo from writing down the higher order $\cal {P\!N\!}$ corrections due to their lengthy form—especially the one of the $2 \,\cal {P\!N\!}$ correction. The complete expressions are given in, e.g., Blanchet (2006) and Kupi et al. (2006).

Our numerical implementation of the $\cal {P\!N\!}$ corrections to the accelerations—and particularly to the corresponding jerk as required for the Hermite integration scheme—has been thoroughly tested against other available codes, such as the ones used by Kupi et al. (2006) and Aarseth (2007). A detailed description of the numerical tests and the results are given in (Berentzen et al. 2008). An alternative implementation specially tailored to the treatment of single massive BH within galactic nuclei can be found in Löckmann & Baumgardt (2008). Our implementation of the $\cal {P\!N\!}$ corrections allows us to turn on and off the different $\cal {P\!N\!}$ orders separately at will. In the current set of models presented here, we apply the selected $\cal {P\!N\!}$ terms to the two BH particles at all times during the simulations.

Within the $\cal {P\!N\!}$ framework, the energy lost due to GW emission is given by (see Blanchet 2006, Equation 171 therein)

$$\begin{eqnarray} \frac{dE}{dt} &=& \frac{4}{5} \frac{G^2 m_1^2 m_2}{c^5 r_{12}^3} \left[({\bf v}_1 {\bf v}_{12}) \left(-{\bf v}_{12}^2 + 2\frac{G m_1}{r_{12}} - 8 \frac{G m_2}{r_{12}} \right) \right. \nonumber \\ &&+ \left. ({\bf n}_{12} {\bf \cdot } {\bf v}_1)({\bf n}_{12}{\bf \cdot }{\bf v}_{12}) \left(3 {\bf v}_{12}^2 - 6\frac{G m_1}{r_{12}}+\frac{52}{3}\frac{G m_2}{r_{12}}\right)\right] \nonumber \\ &&+ \ (1 \leftrightarrow 2) + {\cal {O}}(c^{-7}).\quad \end{eqnarray} \tag{ 16 }$$

Note that this latter equation gives—in contrast to the orbit-averaged expression in Equation (8)—the instantaneous energy loss due to the emission of GWs. We will apply both equations, i.e., Equations (8) and (16), for comparison in our analysis.

In this section, we describe the results of our purely Newtonian simulations, i.e., classical models without any $\cal {P\!N\!}$ corrections. For these fiducial models, the evolution of the binary BHs and the field particles is very similar to the corresponding models reported by Berczik et al. (2006). Due to the relatively high degree of rotation in the models used in this work, the initially axisymmetric nucleus models become dynamically unstable and rapidly develop a rotating triaxial (barlike) structure. As a result of this global instability and due to the dynamical friction against the stellar background, the two BHs are quickly funneled toward the density center of the nucleus, where they form a binary system in our models typically after some Δt ≈ 20, or some 30 Myr.

It actually turns out that certain binaries characteristics, such as the time when the binary forms (t_form), as well as its orbital elements, are very sensitive to the initial conditions of the stellar distribution. Since the eccentricity of the binary is a crucial quantity for its subsequent evolution toward relativistic inspiral and coalescence, it is important to obtain a statistical sample. Therefore, we decided to use nine different realizations of the same galaxy nucleus model, by sampling the underlying distribution function with different initial random seeds. This way, e.g., we find eccentricities of the binaries typically in the range between 0.4 up to 0.99. The binaries in our simulations tend to form with the higher eccentricities around 0.9. However, it cannot be decided at this point whether stochastic encounters between the BHs and the stars, or the global (nonlinear) bar-instability plays the dominant role in determining the final binary properties.

In Figure 1, we show an example for the typical time evolution of the (Keplerian) orbital elements for one of the SMBH binaries in our 50k Newtonian models. In the panels, from left to right, we show the evolution of the semimajor axis a, the eccentricity e, and the energy E of the binary, respectively. Note that the extremely large values of a and e at early times are due to the fact the two SMBHs are still unbound. Only when they have formed a gravitationally bound pair, i.e., when the binaries energy E becomes negative, the eccentricity e remains below a value of 1. As one can see in the middle panel, the eccentricity varies only mildly in the subsequent evolution, i.e., when the binary is hard and interacts with the surrounding field particles mainly by superelastic three-body scattering. Therefore, we calculate some mean eccentricity $\bar{e}$ for each model as indicated by the horizontal line in the middle panel of Figure 1. We find $\bar{e}$ to be a useful quantity to characterize the binary, but it is important to bear in mind that it provides only a first-order approximation, since the real eccentricity may in fact vary with time due to stochastic encounters with field particles.

Zoom In Zoom Out Reset image size

We find that the loss of energy (Figure 1, right panel) due to superelastic three-body encounters is almost constant after roughly t ≳ 30. The straight line shows the result of a linear least-squares fit to E(t) after the binary has formed. In Figure 2, we show the resulting slopes |dE/dt| as derived from such linear fitting as a function of (1 − e²) for our different Newtonian models.⁸ The horizontal lines in Figure 2 indicate the mean values of dE/dt for both the separate and combined sets of our Newtonian 25k and 50k models.

Zoom In Zoom Out Reset image size

Since E ∝ (1/a) for a Keplerian orbit, the rate dE/dt provides a direct measure for the hardening rate of the binary d/dt(1/a). We find that the hardening rate shows an essentially N-independence of both the eccentricity of the binaries' orbit and the number of field particles used in the models.

The latter result is in good agreement with Berczik et al. (2006) who found only small variations in the hardening rate for models with particle numbers ranging from 25k to 1M field particles, which is interpreted as a result of an efficient loss-cone refilling as found in such triaxial, rotating nuclei models.

For our two Newtonian sets of 25k and 50k models, we find mean values for |dE/dt| of about (1.01 ± 0.19) × 10⁻³ and (0.95 ± 0.17) × 10⁻³, respectively. The hardening rate in our models is found to be essentially N-independent within the given error margins. This finding is clearly in accord with the results of Berczik et al. (2006). The mean value of dE/dt for the combined set of models thus results in being roughly $\dot{E}\equiv dE/dt \approx -0.001$.

One can use the results of the previous paragraph to make some simple estimates regarding the binary evolution timescale in the purely Newtonian regime described above. For a Keplerian orbit, we get

$$\begin{equation} \frac{da}{dt} = \frac{2\,a^2}{G\,m_{1}\,m_{2}} \frac{dE}{dt} \,. \end{equation} \tag{ 17 }$$

We define the binary hardening time in the usual way as

$$\begin{eqnarray} \tau _{\mathrm{hard}} & \equiv & \left|{\frac{1}{a}}{\frac{da}{dt}}\right|^{-1} = \left|{\frac{1}{E}}{\frac{dE}{dt}}\right|^{-1} \nonumber \\ & = & {\frac{G\,m_1\,m_2}{2a}} \left|{\frac{dE}{dt}}\right|^{-1}. \end{eqnarray} \tag{ 18 }$$

Writing the rate of change of the binary energy in N-body units as |dE/dt| = 10⁻³ × K, where K ≈ 1, this becomes in physical units

$$\begin{eqnarray} \tau _{\mathrm{hard}}\! &=&\! 9.80 \!\times \! 10^3 K^{-1} {\frac{m_1\,m_2\,r_{c}^{5/2}}{ {G}^{1/2} M_{\mathrm{gal}}^{5/2} a}} \approx 4.62\times 10^4 {\mathrm{yr}} {\frac{m_1\,m_2}{\left(10^8 \,{M}_\odot \right)^2a}}\nonumber \\ &&\times \left({\frac{r_{c}}{10^2\, {\rm pc}}}\right)^{5/2}\left({\frac{M_{\mathrm{gal}}}{10^{11} M_\odot }}\right)^{-5/2} \!\!\left({\frac{a}{1\; {\rm pc}}}\right)^{-1}. \end{eqnarray} \tag{ 19 }$$

This expression shows clearly that a constant rate of energy loss corresponds to a gradually increasing timescale for binary hardening. If we assume that the expression holds for arbitrarily small a, we can predict a time to coalescence of

$$\begin{eqnarray} &&t_{\mathrm{coal}} - t_0 = \tau _{\mathrm{hard}}(a_0)\left({\frac{a_0}{a_{\mathrm{coal}}}} - 1\right), \end{eqnarray} \tag{ 20 }$$

where a₀ = a(t₀) and a_coal = a(t_coal) (Ferrarese & Ford 2005). Setting a₀ = 1 pc and a_coal = (10⁻² ⋅ ⋅ ⋅ 10⁻⁴) a₀ gives coalescence times ranging from ∼10 Myr to 1 Gyr for a galaxy with M_gal = 10¹¹M_☉ and r_c = 100 pc. Thus, we see clearly that the binaries in our galaxy models would have no difficulty reaching very small separations in a Hubble time, even in the absence of relativistic energy loss.

In Figure 3, we plot both the calculated timescale τ_hard for the stellar dynamical hardening, as well as the relativistic coalescence timescale t_gr (see Equation (10)) for orbits of different eccentricity as a function of their semimajor axis. As one can see, the two involved timescales are of the order of only some hundred time units or Myr, respectively, for the range of eccentricities found in our simulations. These timescales are short enough to allow, on average, for the SMBH binaries to coalesce by the time the next galaxy merger would occur.

Zoom In Zoom Out Reset image size

In order to test whether the purely Newtonian stellar dynamics in our models is sufficient enough to bring the SMBH binary into the regime where relativistic effects start to become important, we have done the following test. We first calculate the discrete time derivatives, i.e., finite differences, of the orbital elements and the total energy as directly computed from our simulations. An example is shown in Figure 4 (gray curves). In a second step, we can then use Peters & Mathews (1963) formalism (Equations (2)–(4)) in order to calculate the corresponding orbit-averaged rates, which encode the binary's secular drift due to the emission of GWs at the lowest, quadrupolar, order. The resulting curves are also plotted in Figure 4 (black curve). Finally, we use the $\cal {P\!N\!}$ expression for the (instantaneous) energy loss due to GW emission (Equation (16)).

Zoom In Zoom Out Reset image size

One can see that by the time t ≈ 100 in model units (or some 150 Myr), the two BHs have reached a separation that is small enough for the relativistic effects to set in and reach the same order of magnitude as the (Newtonian) perturbations from the field particles. This analysis demonstrates that with these rotating King models, the SMBH binary quickly (within a cosmological context) hardens to a point where it reaches the relativistic regime and the final parsec problem can be overcome. In the following section, we study the self-consistent evolution of the (nonspinning) SMBH binaries in the field of stars of the galactic nuclei remnant, including the $\cal {P\!N\!}$ corrections up to $2.5 \,\cal {P\!N\!}$ order.

To test if the final parsec problem can really be overcome self-consistently just by stellar-dynamical effects and if the binary BHs in our models eventually reach the relativistic inspiral regime, we take the following approach.

We start with models having exactly the same initial conditions for the set of 25k and 50k models described in the previous section, but this time we take the $\cal {P\!N\!}$ corrections for the two BH particles into account. Most numerical simulations of similar kind which are known to us from the literature are restricted—for simplicity—to the first dissipative $\cal {P\!N\!}$ term in the expansion, i.e., the $2.5 \cal {P\!N\!}$ correction. The effects of the $1 \,\cal {P\!N\!}$ and $2 \,\cal {P\!N\!}$, which both are conservative and are known to result in a pericenter shift, have often been neglected. However, these corrections are of lower order in (v/c) as compared to the $2.5 \,\cal {P\!N\!}$ and thus should clearly affect the binary evolution. Therefore, we decided to run all sets of simulations presented below with (1) only including the $2.5 \,\cal {P\!N\!}$ correction and (2) including the full $\cal {P\!N\!}$ corrections up to $2.5 \,\cal {P\!N\!}$. Since higher order corrections such as the $3 \,\cal {P\!N\!}$ term in the equations of motion are of order c⁻⁶ or smaller, their effect on the dynamical evolution of the BHs are expected to be small as well, and will not affect the main results presented in this work. Hence, we will loosely refer to the models in our Set (2) as the ones with full $\cal {P\!N\!}$ corrections. We should note here that the $\cal {P\!N\!}$ corrections in our runs are applied permanently, independent of the BH separation or their relative velocities (see for comparison Aarseth 2007). Finally, in order to enhance the relativistic effects we use different values of the speed of light c (in models units). We note that by fixing the values for G and c, the ratio between the units of mass and length is also fixed. Therefore, one has only one free parameter for setting the corresponding physical units in contrast to classical N-body simulations, which are scale-free in that respect.

In Figures 5 and 6 (upper panels), we show some examples of the evolution of the semimajor axis a, the eccentricity e, and the total energy E of the binary.⁹ The corresponding time derivatives, i.e., the discrete, Peters & Mathews (1963), and $\cal {P\!N\!}$ as defined in the previous section are shown in the lower panels of these figures, respectively. Both simulations shown in these figures start with identical initial conditions and only differ by the orders of $\cal {P\!N\!}$ corrections taken into account.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

At early times, the dynamical evolution of the binary is qualitatively very similar as compared to the evolution in our purely Newtonian simulations (see Figures 1 and 4): the two BH particles are rapidly driven toward the nuclei center by dynamical friction and the triaxial bar structure in the stellar nucleus. The BHs form a pair at roughly t = 25 (or some 37 Myr, correspondingly). Following the results and interpretation of Berczik et al. (2006), the initially full loss cone around the BH binary is quickly depleted and the evolution of the hard binary thereafter takes place in the empty loss-cone regime. Due to the internal rotation and the triaxial structure of the nucleus, the loss cone is repopulated constantly by centrophilic orbits with a rate which is found to be independent of the number of field particles in the model.

During this phase, the binary evolution is mainly dominated by just the Newtonian stellar-dynamical effects. The average rate of energy loss during this phase is again of the order of −0.001 as has been also found in the Newtonian simulations (Figure 2). Furthermore, this value is independent of the eccentricity with which the binary has formed and remains constant during the Newtonian-dominated phase. That relativistic effects are not yet important during this phase of evolution is supported in the lower panels in Figures 5 and 6 in which we compare the results from the simulations with the ones as expected by applying Peters & Mathews (1963) equations.

The relativistic changes of the orbital elements become of the same order of the Newtonian ones after roughly t ≳ 120 (in Figure 5) and t ≳ 60 (in Figure 6). More generally, we find that the time of transition from the Newtonian-dominated regime to the one in which relativistic effects are of about the same order, depends strongly on the (mean) eccentricity $\bar{e}$ of the binary's orbit. This result is not surprising since as we have already stated before, the dissipation rate due to GW emission becomes already important at earlier times for more eccentric binaries.

Finally, at about times t ≳ 140 and 70 for Figures 5 and 6, respectively, the rate of energy loss increases significantly due to the emission of GWs and becomes very nonlinear. The loss of orbital energy and angular momentum due to the $2.5 \,\cal {P\!N\!}$ corrections leads to the inspiral and circularization of the orbit, and the binary's dynamics in this phase is dominated by the relativistic effects, i.e., the binary almost fully completely decouples from the surrounding nucleus.

It is noteworthy that these are the first systematic astrophysical N-body simulations which actually follow the evolution of the SMBHs from their unbound state, i.e. with an order of kiloparsec separation, down to the subparsec scale close to relativistic coalescence and thus overcome the final parsec problem self-consistently. It seems that the origin of the latter problem partly has been the use of oversimplified models for the galactic nuclei, i.e., assuming spherically symmetric distributions without (net) angular momentum.

We have performed of order 150 different simulations, varying the initial data (statistical realization, particle numbers) and different physical scalings (speed of light, different $\cal {P\!N\!}$ orders). In all models that include the $\cal {P\!N\!}$ equations of motion for the BHs, we are able to follow the binaries' evolution down to the relativistic inspiral and virtually to coalescence. Based on our simulations we find that inclusion of the $1 \cal {P\!N\!}$ and $2 \,\cal {P\!N\!}$ corrections, which have been neglected in most earlier works (again with the notable exception of Aarseth 2003a), strongly affects the dynamical evolution of the binary. Figure 7 shows this, where we directly compare the evolution of 1/a and eccentricity as a function of time for two otherwise equal models with full $\cal {P\!N\!}$ and $2.5 \,\cal {P\!N\!}$ corrections only, respectively. The time of coalescence differs by almost a factor of two between these two cases. We find that part of the reason for this diverging behavior is the different eccentricity with which the binary forms. Using full $\cal {P\!N\!}$ corrections rather than only $2.5 \,\cal {P\!N\!}$ changes the details of the binary's trajectory already during their first encounters. Consequently later on, the binary forms with different orbital elements. Since all close encounters in the system—up to the one leading to the binary formation—are stochastic (due to the dynamical instability of close few-body encounters) we do not find any preferred systematic trend in our models toward higher or lower eccentricities between the two types of simulations as a function of the type of $\cal {P\!N\!}$ corrections used.

Zoom In Zoom Out Reset image size

In Figure 8 we show the distribution of eccentricities $\bar{e}$ obtained from the sample of all our simulations. It clearly shows that the binaries form preferentially with very high eccentricities. This is favored by the transient triaxial feature in our models (see also, Berczik et al. 2006) which brings the two BHs together with relatively small impact parameter. Since galactic mergers lead dominantly to the formation of stellar bars or other triaxial configurations (see, e.g., Naab et al. 2006; Jesseit et al. 2007; Johansson et al. 2009), we expect that high eccentricities in SMBH binaries after galactic mergers are a robust phenomenon occurring in many real cases. Furthermore, the merger of spherical models of galactic nuclei along nearly parabolic trajectories also leads to the formation of SMBH binaries with a similar distribution of high eccentricities (Preto et al. 2009).

Zoom In Zoom Out Reset image size

It is known that during the relativistic inspiral, emission of GWs leads to the circularization of the binary (Peters 1964). In the same Figure 8, we therefore also show the (shaded) histogram for the binary's eccentricities when their semimajor axis has fallen down to 100 Schwarzschild radii R_BH, for the runs with the realistic values of c = 447. We can clearly see that, although circularization has substantially reduced the binaries' eccentricity, the remaining distribution at such small separation is still significantly different from zero. Since the orbital eccentricity of the binaries is very important to predict gravitational waveforms from these binaries, and the expected distribution of their orbital parameters strongly determines the complexity of the data analysis for GW instruments such as the planned LISA satellites (e.g., Babak et al. 2008), we are currently generalizing the present study to a wider range of more realistic initial conditions. It may be interesting to note that this type of statistics provides an astrophysically motivated set of initial conditions for numerical relativity simulations of SMBH binary mergers. During the last couple of years several groups have made significant progress in modeling BH mergers by the solution of the full Einstein equations and of numerical-relativity simulations (see Rezzolla et al. 2008 and references 8–13 therein).

As we have described in Section 3.2, we can distinguish three dynamical phases of SMBH binary evolution: (1) Newtonian phase, (2) mixed phase, and (3) relativistic phase. It is thus interesting to see how much time a binary will spend the corresponding separations. In Figure 9, we plot a histogram showing the time the binary spends at a given semimajor axis normalized to the total binary lifetime, measured from its formation till coalescence. The binary orbit of the model selected for Figure 9 spends some 10% during its evolution at semimajor axis of 10⁻³, or 1 pc in physical units. We fit a Gaussian distribution to the binned data to provide a rough measure of the average value 〈a〉 and of the dispersion around it.

Zoom In Zoom Out Reset image size

In Figure 10, we then plot 〈a〉 as a function of the orbital eccentricity. While the SMBH binaries with high eccentricity tend to have larger 〈a〉, the fitted slope of the double logarithmic plot of 〈a〉 versus (1 − e) is −0.87. We conclude therefore that the average pericenter r_p = 〈a〉(1 − e) of the SMBH orbit at any given average 〈a〉 has a remarkably small variation, i.e., is nearly constant. Therefore, to first order, the orbits of SMBH binary in our models form a one-parameter family. This result is relevant and interesting for the determination of GW backgrounds in the ultralow and low-frequency areas from SMBH binaries in the universe: eccentric SMBH binaries emit GWs with a spectrum of frequencies g(n, e), where n denotes the higher harmonics over the basic mode of the circular orbit (see the Appendix). The harmonic which leads to the maximal emission of GWs n_max(e) is for high eccentricities (e ⩾ 0.8) completely determined by the frequency of motion at pericenter n_peri ∝ (1 − e)^−3/2 (Pierro et al. 2001; Amaro-Seoane et al. 2008). Thus, for every value of average semimajor axis 〈a〉, there exists a typical pericenter value and thus a typical GW spectrum. The expected signals of these objects (in particular, highly eccentric ones) lie in the pulsar timing and lower LISA bands (see for more details, including triple BHs, also Enoki & Nagashima 2007; Hoffman & Loeb 2007).

Zoom In Zoom Out Reset image size

Figure 11 illustrates the cumulative probability to find the SMBH binary with a semimajor axis smaller than a. For each of our runs with full $\cal {P\!N\!}$ terms and realistic value of c = 447, one curve is plotted in the figure. From this information, one can deduce that SMBH binaries with high eccentricity stay longer at larger a, which is an information consistent with the previous figure. The probabilities are normalized for each run separately, therefore these data cannot be used to determine probability distributions of a for samples of SMBH binaries in the universe.

Zoom In Zoom Out Reset image size

In Section 3.1, we estimated the time required for the coalescence of SMBH binaries in our models driven by stellar–dynamical effects only. On the other hand, using the pure relativistic estimate of the coalescence timescale t_gr (Equation (10)), we find for typical semimajor axis values of a_h ≈ 10⁻² ⋅ ⋅ ⋅ 10⁻⁴ and eccentricities of e ≈ 0.5 ⋅ ⋅ ⋅ 0.99 that the resulting timescale covers a range of t_gr ≈ 10⁹ ⋅ ⋅ ⋅ 10⁻⁴ in our model units. The interpretation of these numbers is that for small semimajor axes a and high eccentricities the two SMBHs basically directly plunge, while for larger a and moderate e the relativistic timescale becomes unreasonable large. The evolution and hardening of the binary in such case would therefore mainly be driven by the Newtonian stellar–dynamics over some time period. Eventually, the semimajor axis may reach small enough values for the emission of GWs to become efficient. The binary hardening would then equally be driven by both the classical Newtonian and the relativistic effects, before the latter eventually takes over, i.e., when the binary dynamically decouples from the stellar background.

To get a more precise estimate of the time until relativistic coalescence, we numerically integrate the following differential equation:

$$\begin{equation} T_{\mathrm{merg}} = T_{\mathrm{form}} + \int \left[ \left(\frac{\partial a}{\partial t}\right)_{\mathrm{hard}} + \left(\frac{\partial a}{\partial t} \right)_{\mathrm{P\&M}} \right] dt \,, \end{equation} \tag{ 21 }$$

where the corresponding da/dt are given by Equations (6) and (17), respectively. For the integration, we start with the semimajor axis separation of some 10⁻³ and different eccentricities. For the formation time of the binary, we assume T_form = 15 and 35 as a lower and upper limit, respectively.

In Figure 12, we plot the coalescence time as measured from our simulations for different values of c as a function of eccentricity. The results are then compared to the results obtained from integration of Equation (21). Considering the different times it takes the binaries to form a bound pair in the simulations, our results agree (within the error limits) with the theoretical value.

Zoom In Zoom Out Reset image size

Not surprisingly, we find that the merging time increases with increasing values of c, since the relativistic corrections decrease for simulations starting with the same configuration. However, since the hardening rate due to superelastic scattering with the binary is found to be constant in our simulations, we expect an upper limit of T_merg of roughly 500 time units. It is important to note that in all our simulations the two BHs coalesce in less than 0.5 Gyr after they have formed a bound pair! Therefore, we conclude that our results provide a dynamical substantiation to the picture of prompt SMBH coalescence advocated by Volonteri et al. (2003) and Sesana et al. (2005, 2007). Moreover, we should also note that these times for BH coalescences are of the same order as the average interval between consecutive mergers (see, e.g., Figure 2 in Khochfar & Burkert 2001 and Khochfar & Burkert 2006). Such timescales make it unlikely that there is enough time for a third BH to interact with the binary before it has coalesced. If three-body interactions between SMBHs were the norm, then it would be quite difficult to understand the notably small amount of scatter in the M–σ relation, as well as the scarcity of observational evidence for off-centered nuclei. Note, however, that some fraction of triple interactions of SMBHs may be present in the universe and could cause extremely large eccentricities (e.g., through Kozai oscillations) and become visible through pulsar timing or even in the lower LISA band at comparatively large separations (orbital timescales; Iwasawa et al. 2006, 2008; Hoffman & Loeb 2007; Amaro-Seoane et al. 2008).

Merger events of SMBH binaries are expected to be one of the brightest possible sources of GW emission to be detected by LISA (Danzmann 1997; Phinney 2005), a proposed space-borne GW observatory, scheduled for launch in 2018 or later, is most sensitive to GWs in the low-frequency regime (10⁻⁴–0.1 Hz); for SMBHs of more than about 10⁷ M_☉ LISA will preferentially detect the final merger and ringdown phases (Babak et al. 2008), while for smaller masses earlier evolutionary phases of SMBH evolution will become detectable in the LISA band, in particular for eccentric SMBHs (see discussion and references in Section 4.3).

In most of the present literature, the strain amplitude of the GW is usually estimated under the assumption of circular orbits. This has been motivated by the idea that massive binary systems are expected to have been completely circularized due to the emission of GWs (see, e.g., Peters 1964; Sesana et al. 2005) by the time they enter the LISA frequency band. In this section, we will use our results to test this assumption. The details on how to extract the information on the GW strain amplitude from our numerical simulations calculation are given in the Appendix.

The space-based GW detector LISA will be sensitive to signals of inspiralling compact objects with up to the maximum total mass M₁₂ ≈ 7.5 × 10⁶ M_☉, where M₁₂ is the total mass of the binary. Due to the involved low frequencies of the signal, such sources are out of reach of current ground-based detectors. The physical units adopted in Section 2.1 translate into BH masses of about 10⁹ M_☉ in our simulations. BH binaries in this mass regime, however, will inspiral with frequencies even lower than those of the LISA band and thus would not be detectable by LISA before their actual coalescence.

In order to follow the BH binary inspiral up to frequencies in range of interest for LISA, i.e., f ⩾ 10⁻⁴ Hz, we first need to rescale our models to binary masses of the order of 10⁵ or 10⁶ M_☉. As mentioned earlier, the speed of light in model units scales as c ∝ (R/M)^1/2, where R and M are the units for length-scale and mass, respectively. Since c is already given with a fixed value from our simulations, by changing our unit of mass by factors of 10⁻³ to 10⁻⁴, the length-scale automatically changes to 0.1 or 1 pc, respectively. Finally, following the results of Sesana et al. (2007), we place the system at the redshift z = 4, where LISA's detection rate for objects of that mass range is expected to be significant.

Before actually calculating the strain amplitude of the GW signal, we increase the time resolution of our output data for the BH's trajectory during the late phases of inspiral. This is done by evolving the binary in isolation, starting with initial conditions (in the binaries center of mass frame) as given from the large-scale simulations at a stage when it is already decoupled from the surrounding stellar system, i.e., when the binary evolution is dominated by relativistic dynamics (compare Figure 6). We convince ourselves that the binary is really decoupled from the stellar system by comparing the evolution of its orbital elements to the one of the "low" time resolution. In fact, during the late stages of inspiral, the perturbations from stars vary over timescales much longer than both the orbital period and the time for coalescence (at that moment) which allows us to safely ignore them.

Using the information of the combined time-series of the binary evolution, we then calculate the orbital elements in the $\cal {P\!N\!}$ generalization (Memmesheimer et al. 2004). These allow us to calculate the different modes of the characteristic strain amplitude h_c,n using Equation (A6) from the Appendix.

In Figure 13, we show h_c,n for the first six harmonics of the signal generated by the inspiral in one of our 50k models for (rescaled) binary masses of 2 × 10⁵ M_☉ (left panel) and 2 × 10⁶ M_☉ (right panel). The LISA sensitivity curve S_h(f) displayed in the figures was generated by the Online Sensitivity Curve Generator with default settings, corresponding to a three-year observation period (Larson 2007). We confirm that the BH binaries in the mass range 10⁵–10⁶M_☉ indeed enter the LISA frequency band during the relativistic inspiral in our simulations. The height of h_c,n above the sensitivity curve in Figure 13 provides an approximate indication of the signal-to-noise ratio (S/N; see also the Appendix) and is found to be high enough for detection by LISA.

Zoom In Zoom Out Reset image size

As a consequence of the very high initial eccentricity with which the MBH binaries form in our simulations (compare Figure 8), they reach the LISA band with an eccentricity which may be significantly different from zero—the exact distribution depending on the redshift z of the sources (M. Preto et al. 2009, in preparation). As a result, there will be a significant contribution to the measured power from higher harmonics f_n = n f_orb/(1 + z), with n ⩾ 2 (note that n = 2 for a circular orbit). This can be seen from the plots shown here, where several harmonics contribute significantly to the GW signal and are well above the threshold for detection. Furthermore, we find that, as the binary chirps in the LISA band, it circularizes quite rapidly with the consequence that the n = 2 harmonic always becomes dominant at the last stages of the inspiral.

These findings may have an important impact on the accurate computation of the inspiral waveforms as well as potentially leading to an increased upper limit of the range of detectable masses by LISA. Based on our results, we suggest that orbits with nonvanishing eccentricities should indeed seriously be considered for the LISA data analysis. Further details and discussion of the consequences for the LISA detection of SMBH binary inspirals are beyond the scope of this work and will be published elsewhere (M. Preto et al. 2009, in preparation).