SELF CONSISTENT MODEL FOR THE EVOLUTION OF ECCENTRIC MASSIVE BLACK HOLE BINARIES IN STELLAR ENVIRONMENTS: IMPLICATIONS FOR GRAVITATIONAL WAVE OBSERVATIONS

We consider an MBHB with mass M = M₁ + M₂ (M₁ > M₂) described by its semimajor axis a and eccentricity e. The system is embedded in a purely stellar background with a density profile described by a double power law, as follows:

$$\begin{equation} \begin{array} {l} \rho (r)=\rho _{\rm inf}\left(\frac{r}{r_{\rm inf}}\right)^{-\gamma }\,\,\,\,\,\,r<r_{\rm inf}\\ \rho (r)=\frac{\sigma ^2}{2\pi G r^2}\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,r>r_{\rm inf}. \end{array} \end{equation} \tag{ 1 }$$

Here, r_inf is the influence radius of the binary, identifying the region where the gravitational potential is dominated by the two MBHs, and formally defined as the radius containing a stellar mass equal to M, and ρ_inf is the stellar density at r_inf. The density distribution is normalized to an isothermal sphere for r > r_inf, such a condition sets

$$\begin{equation} \rho _{\rm inf}=\frac{\sigma ^2}{2\pi G r_{\rm inf}^2} \end{equation} \tag{ 2 }$$

and

$$\begin{equation} r_{\rm inf}=(3-\gamma)\frac{GM}{\sigma ^2}\approx 0.8\,{\rm pc}\,(3-\gamma)M_6^{1/2}, \end{equation} \tag{ 3 }$$

where M₆ is the total mass of the binary in units of 10⁶ M_☉, and we made use of the well established M–σ relation in the form (Tremaine et al. 2002)

$$\begin{equation} M_6=0.84\sigma _{70}^4 \end{equation} \tag{ 4 }$$

(σ₇₀ is the velocity dispersion in units of 70 km s⁻¹) to get rid of σ in the last approximation. We identify the region r < r_inf as the inner cusp, and we use γ = 1, 1.5, and2 corresponding to nuclei characterized by cores/weak cusps, mild cusps, and steep cusps, respectively.

N-body simulations of unequal mass binaries in the mass ratio range 0.01–10⁻³ (Baumgardt et al. 2006; Matsubayashi et al. 2007) have shown that dynamical friction is efficient in driving the secondary MBH much deeper than r_inf in the potential well of the primary-cusp system. The two MBH pairs together form an MBHB and continue to harden down to a separation at which the enclosed mass in the binary is of the order of M₂, without significantly affecting the stellar density profile. This is an indication that the hardening is still driven by the dynamical friction exerted by the overall distribution of stars, rather than by close individual encounters with stars intersecting the binary orbit. In our model, we assume that M₂ is driven by dynamical friction down to a separation a₀ where the enclosed stellar mass in the binary is twice the mass of the secondary MBH,

$$\begin{equation} a_0=(3-\gamma)\frac{GM}{\sigma ^2}\left(\frac{q}{1+q}\right)^{1/(3-\gamma)}= r_{\rm inf}\left(\frac{q}{1+q}\right)^{1/(3-\gamma)}, \end{equation} \tag{ 5 }$$

(q = M₂/M₁ is the mass ration of the binary system) without affecting the stellar distribution in the cusp significantly. At that point, three-body interactions take over, and the binary evolution is dictated by individual encounters with stars intersecting its orbit.

On its way to final coalescence starting from a₀, the binary is subject to three main dynamical mechanisms driving its evolution, namely: (1) the erosion of the cusp bound to the primary MBH, (2) the scattering of unbound stars supplied into the binary loss cone by relaxation processes once the stellar distribution is significantly modified by the MBHB, and (3) the emission of GWs. The detailed description of each mechanism has been presented elsewhere, and the reader will be referred throughout this section to the appropriate references for further insights. The focus of the present work is to add them together coherently, to produce sensible, although admittedly very simplistic, evolutionary tracks for MBHBs hardening in stellar environments.

2.2.1. Bound Cusp Erosion

The hardening of an unequal mass MBHB, with an initial semimajor axis a₀ and eccentricity e₀, in a bound cusp was extensively studied by Sesana et al. (2008), hereinafter SHM08. Using their formalism, two differential equations determine the rate of change of the orbital separation and eccentricity:

$$\begin{equation} \frac{da}{dt}\big |_b=-\frac{2a^2}{GM_1M_2}\int _0^\infty \Delta {{\cal E}} \frac{d^2N_{\rm ej}}{da_* dt}da_*, \end{equation} \tag{ 6 }$$

$$\begin{equation} \frac{de}{dt}\big |_b=\int _0^\infty \Delta {e}\frac{d^2N_{\rm ej}}{da_* dt}da_*. \end{equation} \tag{ 7 }$$

Here, a_* is the semimajor axis of a star bound to M₁ and d²N_ej/da_*dt is the number of stars subject to slingshot ejection in the semimajor axis and time intervals [a_*, a_* + da_*], [t, t + dt]. The terms Δe, $\Delta {\cal E}$ are measured from scattering experiments. The ejection rate d²N_ej/da_*dt is instead computed by coupling the numerical results of the experiments to an analytic framework for the binary evolution, which is embedded in a cusp of the form described by Equation (1), as detailed in SHM08. The major finding of SHM08 is that the binary hardens by a factor of ∼10 by extracting the binding energy of the stars in the cusp. During this process, e usually increases by a large factor, depending on the binary mass ratio and on the cusp slope. Results are tabulated in Table 1 of SHM08.

2.2.2. Slingshot of Unbound Stars

The theory of MBHB hardening in a distribution of unbound field stars characterized by a density ρ and a velocity dispersion σ was singled out by Quinlan (1996) and extensively revisited by Sesana et al. (2006), hereinafter SHM06. The binary evolution can be expressed as a function of the dimensionless hardening rate H and eccentricity growth rate K as

$$\begin{equation} \frac{da}{dt}\big |_u=-\frac{a^2G\rho }{\sigma }H, \end{equation} \tag{ 8 }$$

$$\begin{equation} \frac{de}{dt}\big |_u=\frac{aG\rho }{\sigma }HK. \end{equation} \tag{ 9 }$$

The quantities H and K are related to the average energy and angular momentum exchange between the stars and the binary in a single encounter and are computed via extensive three-body scattering experiments, as described, e.g., in SHM06. In general, hardening by scattering of unbound stars becomes effective when the binary reaches the so-called hardening radius, defined as (Quinlan 1996)

$$\begin{equation} a_h\approx \frac{GM_2}{4\sigma ^2}. \end{equation} \tag{ 10 }$$

This is the separation at which the specific binding energy of the binary is of the order of the specific kinetic energy of the field stars. For a > a_h, stars are basically too fast to effectively exchange energy and angular momentum with the binary (soft binary regime); when a < a_h, the binary tends to capture stars in short living weakly bound orbits, kicking them to infinity with v > σ (hard binary regime). The transition soft/hard binary is rather smooth and happens at about a_h. Once the binary is hard, its hardening proceeds at about constant rate, as shown by the H tracks plotted in Figure 3 of SHM06. Perfectly circular binaries tend to stay circular (because of the conservation of the Jacobian integral of motion in the three-body problem), while even slightly eccentric binaries tend to increase their eccentricity, as shown by the K rates plotted in Figure 4 of SHM06.

2.2.3. Gravitational Wave Emission

For our purposes, the effect of GW emission can be modeled in the quadrupole approximation. Under this assumption, the evolution equations for the system are given by Peters & Mathews (1963):

$$\begin{eqnarray} \frac{da}{dt}\big |_{\rm GW}&=&-\!\frac{64}{5}\frac{G^3}{c^5}\frac{M_1M_2M}{a^3(1-e^2)^{7/2}}\left(1+\frac{73}{24}e^2+\frac{37}{96}e^4\right)\nonumber \\ &=&-\!\frac{64}{5}\frac{G^3}{c^5}\frac{M_1M_2M}{a^3}F(e) \end{eqnarray} \tag{ 11 }$$

$$\begin{equation} \frac{de}{dt}\big |_{\rm GW}=-\frac{304}{15}\frac{G^3}{c^5}\frac{M_1M_2M}{a^4(1-e^2)^{5/2}}e\left(1+\frac{121}{304}e^2\right). \end{equation} \tag{ 12 }$$

The function F(e) is defined by the last equality in Equation (11). The shrinking rate is a strong factor of a, meaning that GW-driven hardening is effective only at small separations. The eccentricity evolution rate is also a strong function of a and e itself, and it is always negative. GW emission, therefore, is very effective in circularizing MBHBs, which, in turn, is the reason why little attention has been paid so far to eccentric systems in the context of GW detection.

Having identified the relevant mechanisms at play, we can put the pieces together by writing the evolution of the binary as

$$\begin{equation} \frac{da}{dt}=\sum _i\frac{da}{dt}\big |_i \end{equation} \tag{ 13 }$$

$$\begin{equation} \frac{de}{dt}=\sum _i\frac{de}{dt}\big |_i, \end{equation} \tag{ 14 }$$

where i = b, u, GW labels the three mechanisms considered. Here, we handle MBHBs in stellar environments, and the relevant scale of the stellar distribution is defined by the sphere of influence r_inf of the MBHBs. It is then natural to consider as relevant parameters the stellar density and velocity dispersion at the influence radius, ρ_inf and σ_inf = σ (for an isothermal sphere, the velocity dispersion is independent on radius). It is then instructive to recast the MBHB dynamics in the dimensionless H and K formalism proposed by Quinlan, to compare the dimensionless rates given by each mechanism.

The global evolution of the system can be then written as a generalization of Equations (8) and (9) in the form

$$\begin{equation} \frac{da}{dt}=-\frac{a^2G\rho _{\rm inf}}{\sigma }\sum _iH_i \end{equation} \tag{ 15 }$$

$$\begin{equation} \frac{de}{dt}=\frac{aG\rho _{\rm inf}}{\sigma }\sum _iH_iK_i, \end{equation} \tag{ 16 }$$

where H_i is trivially defined as

$$\begin{equation} H_i=\frac{\sigma }{a^2G\rho _{\rm inf}}\frac{da}{dt}\big |_i \end{equation} \tag{ 17 }$$

and

$$\begin{equation} K_i=a\frac{de}{dt}\big |_i\left(\sum _j\frac{da}{dt}\big |_j\right)^{-1}. \end{equation} \tag{ 18 }$$

We integrate the coupled differential Equations (15) and (16) starting from a₀. As described in Section 2.1(Equation (5)), the value of a₀ is set by the total mass of the binary M, the mass ratio q = M₂/M₁, the cusp slope γ, and the stellar velocity dispersion σ. In our default models, we force σ to obey the M–σ relation given by Equation (4). In this manner, there is a one-to-one correspondence between the mass of the binary and σ, i.e., equally massive binaries are embedded in identical isothermal spheres. Having set the normalization of the stellar distribution and the initial separation a₀, the evolution of the binary depends on the four parameters M₁, q, e₀, and γ. We extensively sample this parameter space as following:

• 1.  
  log(M₁/M_☉) = 2, 3, 4, 5, 6, 7, 8, 9, 10, 11
• 2.  
  q = 1, 1/3, 1/9, ..., 1/729
• 3.  
  e₀ = 0.01, 0.1, 0.3, 0.6, 0.9
• 4.  
  γ = 1, 1.5, 2

for a grand total of 10 × 7 × 5 × 3 = 1050 simulations. Even though the binary evolution under the effect of stellar encounters is basically scale free, the simulation of systems with different absolute masses was necessary to match together the scattering-driven phase to the GW-driven phase, which instead is highly mass dependent. The assumption of different eccentricities at the moment of pairing takes the environmental effects affecting the dynamical friction stage into account. In general, during the merging process, galaxies capture each other on a very eccentric orbit, which is reflected in the initial trajectories of the two MBHs (still at kpc separations at this stage). Dynamical friction against massive, rotationally supported, circumbinary disks has been proven to circularize the orbit (Dotti et al. 2006). However, this is not in general true in gas poor environments, where the process is driven by interaction with the stellar distribution (Colpi et al. 1999), and the eccentricity of the MBHB at the moment of pairing may retain memory of its initial value, or may, in general, be different from zero.

We also consider four alternative models to address the impact of the assumed M–σ relation and the choice of normalizing the efficiency of unbound scatterings to ρ_inf. Let us denote with $\hat{\sigma }$ the value of the velocity dispersion predicted by the M–σ relation for a given M. We consider models with velocity dispersions equal to $0.7\hat{\sigma }$ and $1.3\hat{\sigma }$, which is approximately the range of variance of σ for a given MBH mass measured in the M–σ relation (Haring & Rix 2004). We also run models with two different normalizations for the unbound scattering process: a fast model normalized to 10ρ_inf and a slow model normalized to 0.1ρ_inf. The motivation for this set of runs will be clarified in Section 3.2.

To practically evolve the binary, we make use of the results of the scattering experiments with unbound and bound stars performed in SHM06 and SHM08. In those papers, the quantities Δe, $\Delta {\cal E}$, d²N_ej/da_*dt (for the bound scatterings), H and K (for the unbound scatterings) were recorded on a grid of a and e, covering the relevant dynamical range. The evolution of the coupled differential Equations (15) and (16) is performed by interpolation over the grid as the binary evolves. In SHM06, unbound scatterings were carried out for all the considered mass ratios down to q = 1/243. To complete the sample, we carried out additional experiments for the case q = 1/729. SHM08, instead, focused on unequal MBHBs, with q ≲ 0.1. To complete the mass ratio sample, we ran experiments for the cases q = 1, 1/3. We then have all the bound and unbound scattering experiments' results spanning the q and e₀ range of interest.

Our evolutionary tracks are computed in a self-consistent way, summing together the effects of different mechanisms. However, we should be aware of the several limitations and simplifications we have adopted. One major caveat is the extension of the bound scattering experiments to mass ratios of q > 1/9. It is in fact unlikely that, in such cases, M₂ would reach a₀ without affecting the stellar cusp at all. Cusp disruption would start earlier (especially in the equal mass case), and the cusp erosion phase described here may not be a trustworthy description of reality. However, we find that the impact of the bound cusp erosion on the binary evolution is smaller for larger mass ratios. This is because for q → 1, a₀ ∼ r_inf (see Figure 1), and the impact of the binding energy extraction in the total energy budget of the system becomes less significant. The eccentricity evolution in the cusp erosion phase is only mild when q = 1, 1/3, implying that our approximate treatment would not significantly affect the overall results. Another caveat to bear in mind is that the three-body scattering is a scale-free problem as long as m_* ≪ M₂. Even though we compute "a posteriori" evolutionary tracks for systems with M₁ = 100 M_☉ and q = 1/729, we consider our results meaningful only when M₂ > 100 M_☉. This is also reasonable, since we are interested in the evolution of MBHBs. Moreover, if M₂ is small (<10⁴ M_☉), the amount of stars interacting with the binary is also quite small, i.e., the granularity of the problem increases. Our smooth evolutionary tracks should then be interpreted more as "trends" or "mean evolutions" rather than paths followed by each individual binary. We have also not included the possibility of stellar tidal disruption, which has been shown to be an efficient process in the cusp erosion phase, especially in the mass ratio range 0.01 < q < 0.1 (Chen et al. 2009). In general, the inclusion of tidal disruptions mildly enhances the increase of eccentricity (Chen et al. 2010), because disrupted stars preferentially have a_* < a, i.e., they would drive the binary toward circularization if ejected (see SHM08 for a detailed discussion of this effect). Lastly, there is a somewhat net distinction between bound and unbound scatterings in our formalism, which is certainly oversimplistic, since in reality relaxation processes will mix up the different stellar populations. The appearance of distinctive features in the transition between the bound and the unbound regime can be therefore considered somewhat artificial; the reality would probably be more gentle. We do not believe that this has a major impact on our main results.

Zoom In Zoom Out Reset image size

It is instructive at this point to frame the hardening rate given by Equation (8) in the context of loss-cone refilling theory (Frank & Rees 1976; Lightman & Shapiro 1977; Cohn & Kulsrud 1978). After the binary has depleted all the stars intersecting its orbit (formally defining what is called the "binary loss cone" of the stellar distribution function, see Milosavljevic & Merritt 2003 for a comprehensive review), its further hardening depends on the rate Γ at which such orbits (i.e., the loss cone) are refilled. Loss-cone refilling can proceed by diffusion of stars either in energy (_*) or in angular momentum (j_*). In general, stellar encounters are much more efficient in changing the star angular momentum, and diffusion in the j_* space is the relevant process. The loss-cone refilling depends on the average Δj_* experienced by a star on an almost radial orbit during one orbital period. If this change is larger than the size of the binary loss cone in the angular momentum space, $j_{\rm lc}\sim \sqrt{2GMa}$, then stars are easily scattered back and forth into the loss cone and the loss cone is filled. In this regime, the supply rate of stars from a given distance to the binary r is given by Lightman & Shapiro (1977) and Perets & Alexander (2008):

$$\begin{equation} \frac{d\Gamma _f}{d\,{\rm log}\,r}\approx \frac{a}{r}\frac{N_*(<r)}{P(r)}. \end{equation} \tag{ 22 }$$

Here, P(r) is the typical period of a star on an almost radial orbit coming from a distance r and N_*(<r) is the number of stars enclosed in a sphere of radius r around the MBHB. Let us focus on stars coming from r > r_inf. Assuming that an isothermal sphere N_*(<r) = (M/m_*)(r/r_inf) and that the typical period of a star on a radial orbit is P(r) ∼ r/σ, integrating Equation (22) from r_inf to ∞ and using Equation (2), we get

$$\begin{equation} \Gamma _f\approx 2\pi \frac{M}{m_*}\frac{G\rho _{\rm inf}a}{\sigma }. \end{equation} \tag{ 23 }$$

Let us contrast this result with Equation (8). In the Quinlan formulation, the binary is embedded in a homogeneous stellar field with density ρ. The hardening rate H is then derived writing the interaction rate as a flux of stars through the binary cross section, namely

$$\begin{equation} \Gamma _Q=\frac{\rho }{m_*}\Sigma v, \end{equation} \tag{ 24 }$$

where ρ/m_* is the number density of field stars, v their velocity at infinity (i.e., far from the binary), and Σ the binary cross section. If b is the star impact parameter at infinity and if we assume the encounter to be relevant only for b < b_max, then Σ = πb²_max. Relating b to the maximum approach x to the binary via gravitational focusing (b² = 2GMx/v²) and replacing x_max = a (stars have to cross the binary semimajor axis to exchange energy and angular momentum efficiently), we get

$$\begin{equation} \Gamma _Q=2\pi \frac{M}{m_*}\frac{G\rho a}{v}. \end{equation} \tag{ 25 }$$

By comparing Equations (23) and (25), if we identify the intruder velocity at infinity v with the dispersion velocity in the isothermal sphere σ, we see that our H_u and K_u prescriptions for the scattering of unbound stars correspond to considering the loss cone always full at r_inf.

By normalizing H_u and K_u to σ and ρ_inf, our model implicitly assumes that the loss cone is always full at r > r_inf (i.e., in the so-called pinhole regime) and empty for r < r_inf (i.e., in the so-called diffusive regime, Cohn & Kulsrud 1978). The status of the loss cone then enters as a parameter in our formulation, set by our choice of normalizing the system to ρ_inf. The issue of what is the physical mechanism that keeps the loss cone full at r > r_inf is not addressed in this paper. We will only briefly discuss here the plausibility of such a scenario. After the loss cone is depleted, in the absence of any other physical mechanism, two-body relaxation (Binney & Tremaine 2008) sets the timescale for loss-cone refilling. This is usually longer than the Hubble times in real galaxies (Merritt & Szell 2006), and under this assumption, the loss cone is, in general, in the diffusive regime way beyond r_inf. However, in more realistic situations, a myriad of other physical factors play a substantial role, shortening the loss-cone refilling timescale. It has been shown that axisymmetry and in particular triaxiality (Yu 2002; Merritt & Poon 2004; Berczik et al. 2006) are very effective in repopulating the loss cone, driving stars on chaotic and centrophilic orbits toward the black hole. Moreover, the matter distribution in stellar bulges is far from being smooth. Inhomogeneous concentrations of matter, such as massive star clusters or giant molecular clouds (the so-called massive perturbers), have been proven efficient in perturbing stellar orbits, significantly shortening the loss-cone refilling timescale (Perets & Alexander 2008). We should note that these mechanisms are likely to operate efficiently for r > r_inf. Inside the MBH influence radius, indeed, the presence of a central massive object dominating the gravitational potential tends to force the stellar system toward a more spherical symmetry (and the triaxial structure may be erased in the central region, see, e.g., Merritt & Quinlan 1998). For the same reason, massive perturbers tend to be stripped by the MBH tidal field and would hardly survive in this region. Also, the two-body relaxation timescale increases for decreasing r if the potential is dominated by a central object, because the velocity dispersion of the system, which in this region is the Keplerian velocity around the MBH, scales as r^−1/2, and the relaxation timescale has a strong dependence on the velocity dispersion (Binney & Tremaine 2008). Moreover, inside r_inf, the ejection of stars bound in the cusp depletes those orbits with energy _* < GM/(2r_inf), and since _* diffusion is much less efficient than j_* diffusion the contribution to the loss-cone refilling coming from r < r_inf should be, in general, negligible. It is then reasonable to assume a full loss cone for r > r_inf, and otherwise empty. We shall test the implication of such a hypothesis by varying the normalization used in Equations (15) and (16). We therefore test models normalized to 10ρ_inf and 0.1ρ_inf, which correspond to consider the loss cone full only for r > 10^1/2r_inf and down to r = 10^{−(3−γ)/2}r_inf (where γ is the slope of the inner cusp as defined by Equation (1)), respectively. The consideration of smaller values of the relevant density assumed for the diffusion process (0.1ρ_inf) also serves as a test for the robustness of our results against different outer density profiles. Although our refilling rate is derived for an isothermal distribution, which reasonably fits the measured density profiles of several nearby galaxies and of the Milky Way (Lauer et al. 1995; Dehnen & Binney 1998), many other galaxies show shallower outer density profiles (Ferrarese et al. 1994; Gebhardt et al. 1996; Rest et al. 2001) and the bulk of stars participating to the refilling mechanism may come from slightly larger radii, where the density is lower. We will show (Section 5.2 and Figure 7) that this has a minor impact on our results, changing the eccentricity of the system in the observable GW bands by a factor ≲2.

The evolution of the system depends on the chosen values for the parameters M₁, q, e₀, and γ. Such dependencies are extensively illustrated in Figures 2, 3, and 4. Let us consider each parameter separately, by starting with those defining the masses of the system: M₁ and q.

The evolution of the MBHB as a function of M₁ and for different q is plotted in Figure 2, assuming γ = 1.5 and e_i = 0.1. Since our treatment of stellar scattering is scale free, the value of M₁ affects the system evolution only, by setting the relative gap between a_h and a_GW, which has a mild dependence on M₁. By substituting the M–σ relation in Equation (20), we have in fact a_h/a_GW ∝ M^−1/4, i.e., the gap is larger for lighter systems. This means that, in general, lighter binaries become more eccentric, because, after the bound scattering phase (which is basically scale free by construction) they evolve under the effect of unbound scattering for a larger portion of their dynamical range, as it becomes clear by looking at the K panels of Figure 2. For equal mass binaries with e₀ = 0.1, the eccentricity grows only to 0.15 for M₁ = 10⁹ M_☉ and up to 0.3 for M₁ = 10⁵ M_☉. The mass dependence of the eccentricity growth is also evident for binaries with q = 1/9, while it tends to disappear for lower mass ratios where the eccentricity evolution is dominated by the bound phase. The q dependence of the eccentricity evolution is emphasized by the four different quadrants of Figures 2, 3, and 4. Let us consider again Figure 2. The main result here is that the eccentricity growth in the bound phase, at least when e₀ is small, is in general much larger for lower values of q, as explained in detail in Section 4.1 of SHM08. This is nicely shown by the K panels: as q decreases from 1 to 1/729, the peak in the K rate increases from ∼0.05 to ∼0.6. The value of q also sets the relative weight of the bound and of the unbound phases in the hardening process. Although a₀/a_h is only mildly dependent on q (depending on the cusp slope γ, see Equation (19)), a_h/a_GW is a strong function of q (see Equation (20)) and it can be even less than one for high M₁ and small q (as shown in Figure 1). Therefore, as q decreases, the eccentricity growth is dominated by the bound phase. For e₀ = 0.1, the maximum eccentricity reached by the binary at the end of the stellar driven phase increases from ∼0.2 for q = 1 to ∼0.9 for q = 1/729. The trends with M₁ and q presented in Figure 2 are preserved when changing γ and e₀.

The dependence of the MBHB evolution on e₀ is studied in Figure 3, where we fixed M₁ = 10⁶ M_☉ and γ = 1.5. The binary eccentricity, in general, tends to increase in the scattering phase regardless of the value of e₀. When e₀ = 0.01, binaries with q > 0.1 experience only a mild increase in their eccentricity, up to a value ∼0.2, while binaries with smaller q can reach e > 0.8. Binaries with e₀ > 0.3 tend to reach eccentricities larger than 0.9 regardless of q, M₁, and γ. The evolution of H is basically unaffected by the eccentricity of the system in the bound and unbound phases (in general the average energy subtracted to the binary by the star is not affected by the binary eccentricity, see e.g., SHM06); while K is interestingly larger for higher e₀ when q is large, and vice versa, it decreases with increasing e₀ for small q.

The impact of the assumed slope γ of the density profile is highlighted in Figure 4 for binaries with M₁ = 10⁶ M_☉ and e₀ = 0.1. In general, MBHBs in steeper cusps evolve faster, but to a lower maximum eccentricity, during the scattering phase. As explained in SHM08 this is because, in the scattering process, stars with a_* > a tend to increase e while stars with a_* < a tend to decrease it, and the relative weight of the former is larger in shallower cusps. Moreover, by increasing γ, the dynamical range covered by the scattering process is much shorter (especially for low q), because the a₀–a_GW gap is smaller, and there is less room for significant eccentricity growth. The K rates are mildly affected by γ, with a higher value of γ resulting in smaller K, as explained before. Also note that the absolute value of H in the bound phase increases a lot with γ. This is because the value of a₀ is much smaller for high γ, and the shrinkage in the bound phase is accordingly much faster. In general, for any value of γ, the eccentricity reached at the end of the star scattering phase is >0.7 for q < 0.1, while, again, equal mass binaries experience a less pronounced increase in the eccentricity.

Figures 2, 3, and 4 allow also a detailed study of the evolutionary timescale as a function of the system parameters. Since the bound phase is usually much faster than the unbound one, the evolution timescale of the system is set by $a/\dot{a}=\sigma /(G\rho H a)\propto 1/a$. The coalescence timescale is then set by the bound-GW transition occurring at a_GW. By substituting a_GW given by Equation (20) in the timescale definition above, we get for the coalescence timescale

$$\begin{equation} \tau _c\propto F(e)^{-1/5}(3-\gamma)^{9/5}q^{-1/5}(1+q)^{2/5}. \end{equation} \tag{ 26 }$$

Firstly, we note that τ_c is independent of the absolute value of the MBHB mass, in agreement with Figure 2. This is a consequence of normalizing the stellar distribution outside r_inf to an isothermal sphere obeying the M–σ relation. τ_c also has a very mild dependence on q, increasing by a factor of ∼3 when the mass ratio drops from q = 1 to q = 1/729, as shown by the correspondent panels in Figure 2. High values of the maximum eccentricity accelerate the coalescence by a factor F(e)^1/5, which is ∼5 for e = 0.9; this effect is clear in the (a) and (e) panels of Figure 3. The impact of γ is also quite mild, in spite of the 9/5 exponent, and it modifies τ_c by a factor of ∼3–4 (see (a) and (e) panels in Figure 4). In general, we find 10⁷ yr < τ_c < few × 10⁸ yr.

The evolution of MBHBs in stellar environments has been tackled by several authors by means of full N-body simulations (Milosavljevic & Merritt 2001; Hemsendorf et al. 2002; Aarseth 2003; Makino & Funato 2004; Baumgardt et al. 2006; Matsubayashi et al. 2007; Merritt et al. 2007; Berentzen et al. 2009; Amaro-Seoane et al. 2009, 2010). However, the limited number of particles (N < 10⁶) in such simulations results in very noisy behavior for the binary eccentricity, and it is difficult to draw conclusions about the general trends behind the numerical noise. We can compare our results with N-body simulations carried out in two regimes: q = 1 (equal mass inspirals) and q = 1/1000 (intermediate MBH–MBH inspiral). Milosavljevic & Merritt (2001) carried out numerical integration of equal MBHBs embedded in two merging isothermal cusps (γ = 2). Starting with circular orbits, they find a mild eccentricity increase to a value of ≲0.2 during the stellar-driven hardening phase, consistent with our findings. Merritt et al. (2007) considered equal MBHBs embedded in Dehnen density profiles (Dehnen 1993) with γ = 1.2 with different initial eccentricities. Again, they find that circular binaries tend to stay circular, while eccentric binaries tend to increase their eccentricities in reasonable agreement with the prediction of scattering experiments, and, consequently, with the tracks we presented in Figure 3 for the q = 1 case. Simulations carried out by Aarseth (2003) and Hemsendorf et al. (2002) produce MBHBs with e₀ ≈ 0.8 at the moment of pairing, with e subsequently increasing up to ≳0.95, again consistent with our findings. Berentzen et al. (2009) studied the evolution of equal MBHBs in rotating systems described by a King stellar distribution. They also find quite eccentric binaries at the moment of pairing (e > 0.4), and the subsequent evolution leads to eccentricities larger than 0.95 at which point GW emission takes over, again consistent with our findings. Amaro-Seoane et al. (2009) focused on intermediate MBHBs (M ∼ 10³ M_☉) in massive star clusters. They employ a machinery similar to ours, coupling full N-body simulations to three-body scattering experiments. Their binaries have significant eccentricity (∼0.5–0.6) at the moment of pairing, and the predicted range in the LISA band is 0.1 < e < 0.3, in good agreement with the results shown by the eccentricity maps in Figure 5 that we will describe in the next section. Similar conclusions are reported by Amaro-Seoane et al. (2010). On the small q side, simulations were performed by Baumgardt et al. (2006) and Matsubayashi et al. (2007), assuming a stellar density profile γ = 1.75. When properly rescaled, the eccentricity increase found in both papers agrees surprisingly well with our predictions based on the hybrid cusp erosion model. Unfortunately, we did not find any mention in the N-body literature about the eccentricity evolution for intermediate values of q (0.01 < q < 0.1), and it would be useful to compare our results with N-body simulations in this intermediate range. We find, however, the overall good agreement, at least in the trends, shown at the two extremes of the mass ratio range comforting.

Zoom In Zoom Out Reset image size

To convert our evolutionary track into predictions for GW observations, we proceed as follows. For each MBHB (uniquely defined by M₁, q, γ, and e₀), we convert the a/a₀(t) tracks into a(t) tracks, and then we compute f_k(t), the orbital frequency of the binary, simply by assuming Kepler's law. Having e(t) and f_k(t), we then construct the e(f_k) evolution from the moment of pairing to the final coalescence. We then select two frequencies appropriate for LISA and PTA campaigns, f_LISA and f_PTA, respectively, and evaluate $e|_{f_{LISA}}$ and $e|_{f_{\rm PTA}}$. Remember that for circular binaries, gravitational radiation is emitted at f_GW = 2f_k. We pick f_LISA = 5 × 10⁻⁵ Hz (corresponding to f_GW = 10⁻⁴ Hz, which is approximately the lower bound of the LISA band); on the other hand, we use f_PTA = 5 × 10⁻⁹ Hz (corresponding to f_GW = 10⁻⁸ Hz, which is approximately the frequency at which a 5-to-10 yr PTA campaign will be most sensitive to). The results are shown in Figures 5 and 6 as contour plots $e|_{f_{LISA}}(M_1,q)$ and $e|_{f_{\rm PTA}}(M_1,q)$, for selected values of γ and e₀, as labeled in the figures. Let us start discussing the LISA case. First, we limited M₁ to an upper value of 10⁷ M_☉, since the inspiral of binaries with higher masses will fall outside the LISA band.¹ We also excluded from our contour plots systems with M₂ < 100 M_☉, for the reasons discussed in Section 2.4. The general trend is that lighter unequal mass binaries tend to have larger e when they enter the relevant frequency range. The mass trend is easily explained by the fact that our treatment is largely mass invariant, but the absolute frequency of the system is not. Same stages of the MBHB evolution correspond to progressively lower frequencies as M₁ increases, and since we are in the GW-dominated phase (and thus de/df < 0), more massive systems have lower eccentricities at a given frequency. Binaries with smaller q are in general more eccentric because the eccentricity growth they experience in the stellar scattering phase is larger. If binaries are approximately circular at the moment of pairing (e₀ = 0.01), then the maximum eccentricity in the LISA band is ∼0.2 when M₁ ∼ 10⁴ M_☉ and q ≲ 0.1, and the general trend is largely independent on γ. This is the result of two competitive effects: milder cusps lead to larger values of e, but in this case GW takes over earlier (because the timescale of the three-body scattering evolution is set by the density at r_inf which is lower for milder central cusps, being r_inf itself larger), and it has more time to circularize the orbit before the binary get to f_LISA. If binaries are significantly eccentric at the moment of pairing (e₀ = 0.6 as a study case), then the q dependence in the contour plots almost disappear, because the eccentricity growth in the scattering phase is very efficient irrespective of q when e₀ is large. In this case, light MBHBs (M₁ < 10⁴) may reach f_LISA with eccentricities up to ∼0.5.

Zoom In Zoom Out Reset image size

The situation is even more "dramatic" for PTA observations. PTAs are sensitive to much larger masses (M₁ > 10⁷ M_☉) emitting at much lower frequencies (f ≈ 10⁻⁸ Hz). Systems are, in general, caught far from coalescence (the typical time to coalescence is ∼10⁴ yr, Sesana & Vecchio 2010) and they did not have much time to circularize under the effect of GW emission. Because of this, even if binaries were circular at the moment of pairing, eccentricities can be as high as 0.7 in the PTA band, with the same trend observed for the LISA case (i.e., lighter binaries with smaller q are more eccentric). If e₀ = 0.6, then all the systems with M₁ < 10⁹ M_☉ are expected to have e > 0.5. Again, the results are only mildly dependent on γ. Note that frequencies are computed in the reference frame of the source, the actual observed frequency has then to be appropriately redshifted by a factor (1 + z) according to the redshift of the emitting system. This means that for high redshift sources, an observed frequency of 10⁻⁴ Hz corresponds to a higher intrinsic frequency. High redshift sources may then have milder eccentricities in the observable band, which may be relevant for LISA sources (whereas typical PTA sources are at z < 1).

We want to check at this point how the M–σ normalization and the tuning of the loss-cone refilling efficiency to r_inf impact on our results. Selected cases of alternative models are presented in Figure 7. The left-hand plot is representative of LISA sources. We see that changing σ merely shifts the timescale of the binary evolution. For larger σ, the system is more compact, the timescale for three-body scattering is shorter, GW emission takes over later, and consequently the residual e at f_LISA is larger. However, as shown by the e(f_k) tracks, this is at most a factor-of-2 effect. Changing the normalization ρ in the loss-cone refilling process has instead a major impact on the evolutionary timescale and on the eccentricity evolution of the system, but still the residual eccentricity at f_LISA is basically unaffected when e₀ = 0.01, and it changes by at most a factor of 3 in the e₀ = 0.6 case. In the right-hand plot, we consider instead the typical PTA source. All the considerations made for the LISA case still hold, and in the relevant frequency range 3 × 10⁻⁹ Hz < f < 3 × 10⁻⁸ Hz, the expected e changes by at most a factor of 2. We therefore consider our results quite robust irrespective of the assumed normalizations.

Zoom In Zoom Out Reset image size

As a final step, we quantify the eccentricity distribution of GW sources in the relevant frequency range resulting by applying our eccentricity evolution scheme to standard MBHB population models.

For the LISA case, we use two of the models utilized by the LISA parameter estimation task force (Arun et al. 2009): in the first case seeds are light (M ≳ 100 M_☉, VHM model; Volonteri et al. 2003), being the remnant of the first POPIII star explosions (Madau & Rees 2001); in the second case already quite heavy (M ≳ 10⁴ M_☉) seed BHs form by direct collapse of massive protogalactic discs (BVR model; Begelman et al. 2006). We ran 50 Monte Carlo realizations of each model, producing 50 catalogs of coalescing binaries over a period of three years. We then estimate the signal-to-noise ratio (S/N) of each binary in the LISA detector by assuming circular inspiral and computing the waveform to the second post-Newtonian order. We then consider only those events resulting in an S/N >8 in the detector, and we compute the expected eccentricity distribution at f_LISA. Results are shown in Figure 8. In the upper and in the middle panels, we plot the contour plots of the differential distribution of GW sources as a function of M₁ and q, d²N/dM₁dq, averaged over the 50 Monte Carlo realizations, superposed to the contour plots for $e|_{f_{{LISA}}}(M_1,q)$. The two observed MBHB populations are extremely different: in the VHM model, the bulk of sources have M ∼ 10⁴ M_☉ with q ∼ 0.1; while in the BVR model, most of the sources have M > 10⁴ M_☉ and q ≈ 1. The resulting eccentricity distributions are plotted in the lower panel, where we plot the probability density function p(e) against e for the observed population. When e₀ = 0.01 (binaries approximately circular), eccentricity is expected to be <10⁻² in the BVR model, with a peak at about 2 × 10⁻³, but a broad eccentricity spectrum covering the range 10⁻³–0.2 is expected in the VHM case. In this latter scenario, in fact, sources are on average less massive and with low q, a condition that maximizes the eccentricity increase during the stellar scattering phase. If e₀ is already large (0.6 in our study case), then the observed eccentricity at  f_LISA is peaked at ∼0.1 for the BVR case and at ∼0.4 in the VHM case.

Zoom In Zoom Out Reset image size

In exploring the consequences for PTA observations, we adopt the standard Tu-SA population model employed by Sesana et al. (2009), where merging galaxies are populated by MBHs according to the M − M_bulge in the form given by Tundo et al. (2007) and accretion is triggered onto the more MBH before the final coalescence. The reader is referred to Sesana et al. (2009) for details. We ran 50 Monte Carlo realizations of the model (assuming binaries in circular orbit) and we pick only the individually resolvable sources generating a timing residual larger than 1 ns. Again, the obtained differential distribution of the individually resolvable sources, d²N/dM₁dq, averaged over the 50 Monte Carlo realizations, is superposed to the contour plots for e|_PTA(M₁, q) in Figure 9. The source distribution is strongly peaked around M₁ = 10⁹ M_☉ and q = 0.1, with a long tail extending to q = 10⁻³. If binaries are approximately circular at the moment of pairing (e₀ = 0.01), then the expected p(e) is basically flat in the range [0.03, 0.3], while for e₀ = 0.6, p(e) has a sharp peak in the range [0.5, 0.7], highlighting the possible significant impact of eccentricity for PTA observations.

Zoom In Zoom Out Reset image size